initial ridging push

ridging/README.md

@@ -0,0 +1,23 @@
+Installation
+#############
+
+Software requirements:
+Python 3.9 with following packages:
+  - numpy
+  - xarray (including dask, netCDF4 and bottleneck)
+  - pyyaml
+
+Use the provided file conda-env_NOCOS-RIDGING.yml to install a suitable conda environment.
+
+Running  TODO redi
+############
+1) Adjust the configuration file config_RIOcalc_default.yml or make a new one and provide its name inside ridg_engine_NOCOS.py
+2) conda activate ridgNOCOS
+3) python3 ridg_engine_NOCOS.py
+
+Info
+#############
+
+
+
+

ridging/calculate_RID.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Andrea Gierisch, DMI
+@author: Ilja Maljutenko, MSI
+
+"""
+
+def calc_RID_multicat(icedata,configs):
+
+    import numpy as np
+    import warnings
+    import time # for sleep and cputime measurement
+    import sys # For flushing of buffer to stdout
+
+    # Simpler names for relevant variables
+    ############################################
+    timestr=configs['coordinates']['time_name']
+    nistr=configs['coordinates']['ncat_name']
+#    print(nistr)
+    nt = icedata[configs['coordinates']['time_name']].shape[0] # Number of time steps
+    ni = icedata.sit_c.shape[1] 
+    ny = icedata.sit_c.shape[2] # Number of grid points in y-direction
+    nx = icedata.sit_c.shape[3] # Number of grid points in x-direction
+    sit= icedata.sit_c.sel(ncatice=slice(ni-7,ni))
+    sic= icedata.sic_c.sel(ncatice=slice(ni-7,ni))
+ #   print(sic)
+    sits=sit.sum(dim='ncatice')
+    sics=sic.sum(dim='ncatice')
+# Set up an empty output array
+#    ridgfinal=np.array(np.nan*np.zeros([ny,nx]))
+#    print(  nt, ny,nx)
+    thr_thck=0.01
+    mask_t=sits.where(sits > thr_thck, 0) 
+    thr_conc=0.01
+    mask_c=sics.where(sics > thr_conc, 0) 
+    mask_ct=mask_c #*mask_t
+    mask_ct=mask_ct.sum(dim=timestr)/(nt)
+#    print(nt)    
+#    mask_ct=mask_ct.sum(dim=timestr)/(nt)
+    ridgfinal=mask_ct                             
+    return(ridgfinal)
+
+
+def calc_RID_intcat(icedata,configs):
+
+    import numpy as np
+    import warnings
+    import time # for sleep and cputime measurement
+    import sys # For flushing of buffer to stdout
+#    print(icedata)
+    # Simpler names for relevant variables
+    ############################################
+    timestr=configs['coordinates']['time_name']
+    nt = icedata[configs['coordinates']['time_name']].shape[0] # Number of time steps
+    ny = icedata.sit.shape[1] # Number of grid points in y-direction
+    nx = icedata.sit.shape[2] # Number of grid points in x-direction
+    sit= icedata.sit
+    sic= icedata.sit
+# Set up an empty output array
+#    ridgfinal=np.array(np.nan*np.zeros([ny,nx]))
+#    print(  nt, ny,nx)
+    thr_thck=0.7 #0.6
+    mask_t=sit.where(sit > thr_thck, 0) 
+    thr_conc=0.95
+    mask_c=sic.where(sic > thr_conc, 0) 
+#    mask_ct = (sit> thr_thck) & (sic > thr_conc)
+#    mask_ct=mask_ct.where(mask_ct > 0, 1)     
+    mask_ct=mask_c*mask_t
+    ridged_ice=mask_ct.sum(dim=timestr)/(nt)
+#    mask_ct=mask_c*mask_t
+#    mask_ct_p=mask_t.sum(dim=timestr) #/(nt)
+#    ridgfinal=mask_ct_p                             
+    return(ridged_ice)
+
+################################################
+################################################
+
+if __name__ == "__main__":
+    # If this script is called directly without RIOengine_NOCOS-py, do everything necessary using config from config_RIOcalc.yml
+    from read_config import read_configfile
+    from read_ice_data import read_ice_data
+    from output_RID import save_RID_toNetCDF
+
+    configs=read_configfile()
+    icedata=read_ice_data(configs)
+    riodata=calc_RID_multicat(icedata,configs)
+    save_RID_toNetCDF(riodata,icedata,configs)

ridging/conda-env_NOCOS-RIDGING.yml

@@ -0,0 +1,19 @@
+# To create the environment:
+# #########################
+# conda env create -f conda-env_NOCOS-RIDGING.yml
+#    To update if yml-file has changed:
+#    conda env update --file conda-env_NOCOS-RIDGING.yml --prune 
+# To use the environment:
+# ######################
+# conda activate ridgingNOCOS
+name: ridgingNOCOS
+channels:
+  - conda-forge
+dependencies:
+  - python>=3.9
+  - numpy
+  - xarray
+  - pyyaml
+  - dask        # for xarray
+  - netCDF4     # for xarray
+  - bottleneck  # for xarray

ridging/config_ridg_calc_REF.yml

@@ -0,0 +1,49 @@
+# Name of the input file(s) containing ice information
+# You can specify one filename or a string glob in the form "path/to/my/files/*.nc" to merge multiple files (e.g. along time dimension).
+ice_filename: "/home/ilja/hpcshare/nemo_est/out_ice/N_run7_ice.selgor.mergetime*.nc" 
+
+# Information about variable oordinates
+coordinates:
+   ncat: 15          # Number of ice thickness categories
+   ncatnl: 5          # Number of ice thickness categories   
+   ncat_name: ncatice     # Name of ncat dimension (if ncat > 1)
+   time_name: time_counter   # Name of time dimension
+   lat_name: x    # Name of latitude field (only for plotting)
+   lon_name: y    # Name of longitude field (only for plotting)
+
+# Which (relevant) variables are contained in your file?
+variables:
+   siconc_name:  icefrac         # Percentage of grid cell covered by sea ice
+   sithick_name: icethic         # Actual (floe) thickness of sea ice (not volume)
+   sivol_name:            # Total volume of sea ice divided by grid-cell area
+   siitdconc_name:  iceconc_cat # Percentage of grid cell covered by each ice-thickness categor
+   siitdthick_name: icethic_cat      # Actual (floe) thickness of sea ice in each category
+   siitdvol_name:    # Ice volume in each categroy divided by whole grid cell area
+   siage_name:        # Age of sea ice
+   sisali_name:       # Mean sea-ice salinity of all sea ice in grid cell
+
+# Method how RIDGING should be calculated
+RIDGE_method:
+   Lthickconc: True  # Calculate RIDGING purely based on ice thickness info
+   Ldthick:   False  # Calculate RIDGING based on ice thickness  difference
+   Ldproc:   False  # Calculate RIDGING based on RIDGING auxilary parameters
+
+# Select time or space range
+sel: # Not implemented yet
+   t1: 2011-01-01 
+   t2: 2011-04-01
+   x1: 1
+   x2: 151
+   y1: 1
+   y2: 181
+#   x1: 2
+#   x2: 80
+#   y1: 30
+#   y2: 120
+
+# Output
+output:
+   filename_automatic: True
+   filename:       test_output_rid     # *.nc (Only the name, not a full path)
+   output_folder: .   # If different from working directory
+   title:         Ridging Index Outcome (Ridging) calculated from  XX # Given as global attribute to the RIO netcdf file

ridging/output_RID.py

@@ -0,0 +1,68 @@
+
+def save_RID_toNetCDF(ridg_data,icedata,configs):
+
+    import numpy as np
+    import datetime
+    import os
+
+    # RIDmethod as string
+    rid_method=[i for i in configs['RIDGE_method'] if configs['RIDGE_method'][i]==True][0]
+    # Create the output filename
+    if configs['output']['filename_automatic']==True:
+        outfile=configs['output']['output_folder']+'/RID-'+rid_method+'_'+os.path.basename(configs['ice_filename']).replace('*','_XXX').replace('?','X')
+    else:
+        outfile=configs['output']['output_folder']+'/'+configs['output']['filename']
+
+    # Make an empty dataset with basic coordinates for the netcdf file, copied from the icedata-file
+    timespace_coords=[configs["coordinates"]["lon_name"], configs["coordinates"]["lat_name"], configs["coordinates"]["time_name"]]
+    ridds=icedata[timespace_coords] # Create an "empty" DataSet with same time and space dimentions as in icedata
+
+    # Add coordinate for shipclasses
+#    shipclasses = configs['output']['shipclasses']
+#    ridds.coords["shipclass"] =  np.array(configs['output']['shipclassNUMs'], dtype='int32') # Would be nice with ship classes as string values, but cdo cannot handle that. Therefore we use integer values.
+
+    # Add the RID variable with dimensions, attributes and data from riddata
+    rid_dims = (icedata.sic_c.dims[2],icedata.sic_c.dims[3]) # (time, shipclass, y, x)
+    ridds["RID"] = (rid_dims, ridg_data)
+    rid_attrs={'long_name':'ridging probability',
+                'units':'[]',
+                'comment':"Ridging probability based on ice thickness and concentration ",
+                }
+    ridds.RID.attrs=rid_attrs
+
+    ## Try to correct the attribute 'coordinates' of the RIO variable in netcdf. For DMI data, xarray's automatic does it wrongly. But this doesn't work:
+    #rio_coords=configs["coordinates"]["lon_name"] +' '+ configs["coordinates"]["lat_name"] +' shipclass '+ configs["coordinates"]["time_name"]
+    #riods["RIO"] = (rio_dims, riofinal, {'coordinates':rio_coords})
+    #riods.RIO.encoding["coordinates"]=rio_coords
+
+    # Set global attributes of the netcdf file
+    ridds.attrs={
+    'title':configs['output']['title'],
+    'histroy':'Created by ridg_engine_NOCOS.py on '+datetime.datetime.today().strftime("%Y-%m-%d"),
+    'description':'Used ridging calculation method: '+rid_method ,
+    'institution':"Marine Systems Deprtment at TalTech, MSI",
+    'references':'-'
+    }
+
+    # Write netcdf file
+    ridds.to_netcdf(outfile)
+#, unlimited_dims=configs["coordinates"]["time_name"], encoding={'RID':{'_FillValue':np.nan},'time':{"dtype": "double", 'units': "days since 1900-01-01 00:00:00"}})
+    print("RIDGING output written to: "+outfile)
+
+
+def plot_RID(ridg_data,icedata,configs):
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import datetime
+    import os
+
+        #print(ridg_data)
+    plt.imshow(ridg_data, cmap='Blues', origin='lower',vmin=0., vmax=1. )
+    plt.title('Ridged Ice Distribution')
+    plt.colorbar(label='Ridged Ice Probability')
+    plt.xticks([])  # Remove x-axis ticks
+    plt.yticks([])  # Remove y-axis ticks
+    plt.gcf().set_facecolor('#F0F0F0')  # Use a hex code for light gray
+    plt.savefig('figure_plot.png', dpi=300, bbox_inches='tight')
+    return

ridging/read_config.py

@@ -0,0 +1,34 @@
+def read_configfile(*argv):
+    import yaml
+
+# Open the config file
+######################
+    if len(argv)==1:
+        config_filename=argv[0]  # Given as input argument
+    elif len(argv)==0:
+        config_filename='conda-env_NOCOS-RIDGING.yml' # Default filename
+    else:
+        raise RuntimeError("READ_CONFIGFILE can only be called with zero or one argument, not with "+str(len(argv)))
+
+    with open(config_filename,"r") as ymlfile:
+         configs = yaml.load(ymlfile,  Loader=yaml.Loader)
+
+
+    if configs['coordinates']['ncat'] == 1:
+        configs['multiCAT']=False                 
+    else:
+        configs['multiCAT']=True     
+
+
+# f) Check output settings
+    if configs['output']['filename_automatic']!=True and configs['output']['filename']==None:
+        raise RuntimeError("You need to set FILENAME_AUTOMATIC to True if you don't provide an output filename.")
+    if configs['output']['output_folder']==None:
+        configs['output']['output_folder']="."  # If no output folder is given, set it to the current working directory
+
+    print("These are the configurations you have chosen:")
+    print(configs)
+    print()
+
+    return configs
+

ridging/read_ice_data.py

@@ -0,0 +1,51 @@
+
+#from read_config import read_configfile
+#configs=read_configfile()
+
+######################################################
+
+def read_ice_data(configs):
+    import xarray as xr
+    print(configs['multiCAT']) 
+    ds = xr.open_mfdataset(configs['ice_filename']) # Will try to infer automatically in which dimension (e.g. time) the files should be concatenated.
+#    icedataraw = xr.open_mfdataset(configs['ice_filename'], concat_dim=configs['coordinates']['time_name'])
+
+
+    time_slice = slice(configs['sel']['t1'], configs['sel']['t2'])  
+    lon_slice=slice(configs['sel']['x1'], configs['sel']['x2'])  
+    lat_slice=slice(configs['sel']['y1'], configs['sel']['y2'])  
+    icedataraw=ds.sel(time_counter=time_slice,y=lat_slice, x=lon_slice)
+
+
+    # Sanity checks
+    ######################
+    # ToDO: Check if all provided variable names actually exist in icedataraw
+    # to be implemented
+
+    # Check if NCAT matches to dimensions in icedataraw:
+    if configs['multiCAT']:
+        # Check siitdconc
+        if icedataraw[configs['variables']['siitdconc_name']][configs['coordinates']['ncat_name']].shape[0] != configs['coordinates']['ncat']:
+            raise RuntimeError("The given number for NCAT doesn't match the dimension 'ncat_name' for siitdconc_name.")
+
+
+    ## Select relevant variables from input file(s) and assign them standard names
+    ##########################################
+    if configs['multiCAT']:
+        relevantvars_list=[configs['variables']['siitdconc_name'], configs['variables']['siitdthick_name']]
+        rename_vardict={configs['variables']['siitdconc_name']:'sic_c', configs['variables']['siitdthick_name']:'sit_c'}
+    elif configs['multiCAT']==False:
+        relevantvars_list=[configs['variables']['siconc_name'], configs['variables']['sithick_name']]
+        rename_vardict={configs['variables']['siconc_name']:'sic', configs['variables']['sithick_name']:'sit'}
+
+    if configs['RIDGE_method']['Lthickconc']:
+        relevantvars_list.append(configs['variables']['siconc_name'])
+        relevantvars_list.append(configs['variables']['sithick_name'])
+
+    if configs['RIDGE_method']['Ldthick']:
+        relevantvars_list.append(configs['variables']['sithick_name'])
+
+    icedata = icedataraw[relevantvars_list].rename(rename_vardict)
+#    icedata = icedataraw[relevantvars_list]
+
+    return(icedata)

ridging/ridg_engine_NOCOS.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+
+@author: Andrea Gierisch, DMI
+@author: Ilja Maljutenko, MSI
+
+"""
+# Load RID functions
+from read_config import read_configfile
+from read_ice_data import read_ice_data
+from calculate_RID import calc_RID_multicat
+from calculate_RID import calc_RID_intcat
+from output_RID import save_RID_toNetCDF
+from output_RID import plot_RID
+
+
+# Read config file
+configs=read_configfile('config_ridg_calc_REF.yml')
+
+# Read ice data
+icedata=read_ice_data(configs)
+print(icedata)
+# Calculate ridging
+if configs['multiCAT']==True:
+    ridg_data=calc_RID_multicat(icedata,configs)
+elif configs['multiCAT']==False:
+    ridg_data=calc_RID_intcat(icedata,configs)
+
+
+save_RID_toNetCDF(ridg_data,icedata,configs)
+plot_RID(ridg_data,icedata,configs)
+
+
+
+