Development of samir xarray

code/modspa_samir.py

@@ -10,7 +10,7 @@ Usage of the SAMIR model in the Modspa framework.
 import os  # for path exploration
 import csv  # open csv files
 from fnmatch import fnmatch  # for character string comparison
-from typing import List, Tuple  # to declare variables
+from typing import List, Tuple, Union  # to declare variables
 import xarray as xr  # to manage dataset
 import pandas as pd  # to manage dataframes
 import numpy as np  # for math and array operations
@@ -18,7 +18,42 @@ import rasterio as rio  # to open geotiff files
 import geopandas as gpd  # to manage shapefile crs projections
 from parameters.params_samir_class import samir_parameters
 
-def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Dataset, land_cover_raster: str) -> xr.Dataset:
+
+def xr_maximum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:
+    """
+    Equivalent of `numpy.maximum(ds, value)` for xarray DataArrays
+
+    ## Arguments
+    1. ds: `xr.DataArray`
+        datarray to compare
+    2. value: `Union[xr.DataArray, float, int]`
+        value (scalar or dataarray) to compare
+
+    ## Returns
+    1. output: `xr.DataArray`
+        resulting dataarray with maximum value element-wise
+    """
+    return xr.where(ds <= value, value, ds)
+
+
+def xr_minimum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:
+    """
+    Equivalent of `numpy.minimum(ds, value)` for xarray DataArrays
+
+    ## Arguments
+    1. ds: `xr.DataArray`
+        datarray to compare
+    2. value: `Union[xr.DataArray, float, int]`
+        value (scalar or dataarray) to compare
+
+    ## Returns
+    1. output: `xr.DataArray`
+        resulting dataarray with minimum value element-wise
+    """
+    return xr.where(ds >= value, value, ds)
+
+
+def rasterize_samir_parameters(csv_param_file: str, empty_dataset: xr.Dataset, land_cover_raster: str, chunk_size: dict) -> Tuple[xr.Dataset, dict]:
     """
     Creates a raster `xarray` dataset from the csv parameter file, the land cover raster and an empty dataset
     that contains the right structure (emptied ndvi dataset for example). For each parameter, the function loops
@@ -27,14 +62,18 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
     ## Arguments
     1. csv_param_file: `str`
         path to csv paramter file
-    2. parameter_dataset: `xr.Dataset`
-        empty dataset that contains the right structure (emptied ndvi dataset for example)
+    2. empty_dataset: `xr.Dataset`
+        empty dataset that contains the right structure (emptied ndvi dataset for example).
     3. land_cover_raster: `str`
         path to land cover netcdf raster
+    4. chunk_size: `dict`
+        chunk_size for dask computation
 
     ## Returns
     1. parameter_dataset: `xr.Dataset`
         the dataset containing all the rasterized Parameters
+    2. scale_factor: `dict`
+        dictionnary containing the scale factors for each parameter
     """
     
     # Load samir params into an object
@@ -44,7 +83,13 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
     class_count = table_param.table.shape[1] - 2  # remove dtype and default columns
     
     # Open land cover raster
-    land_cover = xr.open_dataarray(land_cover_raster)
+    land_cover = xr.open_dataarray(land_cover_raster, chunks = chunk_size)
+    
+    # Create dataset
+    parameter_dataset = empty_dataset.copy(deep = True)
+    
+    # Create dictionnary containing the scale factors
+    scale_factor = {}
     
     # Loop on samir parameters and create 
     for parameter in table_param.table.index[1:]:
@@ -53,6 +98,10 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
         parameter_dataset[parameter] = land_cover.astype('i2')
         parameter_dataset[parameter].attrs['name'] = parameter
         parameter_dataset[parameter].attrs['description'] = 'cf SAMIR Doc for detail'
+        parameter_dataset[parameter].attrs['scale factor'] = str(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])
+        
+        # Assigne value in dictionnary
+        scale_factor[parameter] = 1/int(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])
         
         # Loop on classes to set parameter values for each class
         for class_val, class_name in zip(range(1, class_count + 1), table_param.table.columns[2:]):
@@ -62,23 +111,24 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
             parameter_dataset[parameter].values = np.where(parameter_dataset[parameter].values == class_val, round(table_param.table.loc[table_param.table.index == parameter][class_name].values[0]*table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0]), parameter_dataset[parameter].values).astype('i2')
     
     # Return dataset converted to 'int16' data type to reduce memory usage
-    # Scale factor for calculation is stored in the samir_parameters object
-    return parameter_dataset
+    # and scale_factor dictionnary for later conversion
+    return parameter_dataset, scale_factor
 
 
-def setup_time_loop(calculation_variables: List[str], calculation_constant_values: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset, xr.Dataset]:
+def setup_time_loop(calculation_variables_t1: List[str], calculation_variables_t2: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:
     """
-    Creates three temporary `xarray Datasets` that will be used in the SAMIR time loop.
+    Creates two temporary `xarray Datasets` that will be used in the SAMIR time loop.
     `variables_t1` corresponds to the variables for the previous day and `variables_t2`
     corresponds to the variables for the current day. After each loop, `variables_t1`
-    takes the value of `variables_t2`. `constant_values` corresponds to model values
-    that are not time dependant, and therefore stay constant during the SAMIR time loop.
+    takes the value of `variables_t2` for the corresponding variables.
 
     ## Arguments
-    1. calculation_variables: `List[str]`
+    1. calculation_variables_t1: `List[str]`
         list of strings containing the variable names
-    2. calculation_constant_values: `List[str]`
-        list of strings containing the constant value names
+        for the previous day dataset
+    2. calculation_variables_t2: `List[str]`
+        list of strings containing the variable names
+        for the current day dataset
     3. empty_dataset: `xr.Dataset`
         empty dataset that contains the right structure
 
@@ -87,35 +137,157 @@ def setup_time_loop(calculation_variables: List[str], calculation_constant_value
         output dataset for previous day
     2. variables_t2: `xr.Dataset`
         output dataset for current day
-    3. constant_values: `xr.Dataset`
-        output dataset for constant values
     """
     
     # Create new dataset
     variables_t1 = empty_dataset.copy(deep = True)
     
     # Create empty DataArray for each variable
-    for variable in calculation_variables:
+    for variable in calculation_variables_t1:
         
         # Assign new empty DataArray
         variables_t1[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))
         variables_t1[variable].attrs['name'] = variable  # set name in attributes
     
-    # Create new copy of the variables_t1 dataset
-    variables_t2 = variables_t1.copy(deep = True)
-    
-    # Create dataset for constant values
-    constant_values = empty_dataset.copy(deep = True)
+    # Create new dataset
+    variables_t2 = empty_dataset.copy(deep = True)
     
-    # Create empty DataArray for each value
-    for value in calculation_constant_values:
+    # Create empty DataArray for each variable
+    for variable in calculation_variables_t2:
         
         # Assign new empty DataArray
-        constant_values[value] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
-        constant_values[value].attrs['name'] = value  # set name in attributes
-        constant_values[value].attrs['description'] = 'Values which stays constant during the SAMIR time loop'  # set description in attributes
+        variables_t2[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))
+        variables_t2[variable].attrs['name'] = variable  # set name in attributes
+    
+    return variables_t1, variables_t2
+
+
+def prepare_outputs(empty_dataset: xr.Dataset, additional_outputs: List[str] = None) -> xr.Dataset:
+    """
+    Creates the `xarray Dataset` containing the outputs of the SAMIR model that will be saved.
+    Additional variables can be saved by adding their names to the `additional_outputs` list.
+
+    ## Arguments
+    1. empty_dataset: `xr.Dataset`
+        empty dataset that contains the right structure
+    2. additional_outputs: `List[str]`
+        list of additional variable names to be saved
+
+    ## Returns
+    1. model_outputs: `xr.Dataset`
+        model outputs to be saved
+    """
+    
+    # Evaporation and Transpiraion
+    model_outputs = empty_dataset.copy(deep = True)
+    
+    model_outputs['E'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['E'].attrs['units'] = 'mm'
+    model_outputs['E'].attrs['standard_name'] = 'Evaporation'
+    model_outputs['E'].attrs['description'] = 'Accumulated daily evaporation in milimeters'
+    model_outputs['E'].attrs['scale factor'] = '1'
+    
+    model_outputs['Tr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['Tr'].attrs['units'] = 'mm'
+    model_outputs['Tr'].attrs['standard_name'] = 'Transpiration'
+    model_outputs['Tr'].attrs['description'] = 'Accumulated daily plant transpiration in milimeters'
+    model_outputs['Tr'].attrs['scale factor'] = '1'
+    
+    # Soil Water Content
+    model_outputs['SWCe'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['SWCe'].attrs['units'] = 'mm'
+    model_outputs['SWCe'].attrs['standard_name'] = 'Soil Water Content of the evaporative zone'
+    model_outputs['SWCe'].attrs['description'] = 'Soil water content of the evaporative zone in milimeters'
+    model_outputs['SWCe'].attrs['scale factor'] = '1'
+    
+    model_outputs['SWCr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['SWCr'].attrs['units'] = 'mm'
+    model_outputs['SWCr'].attrs['standard_name'] = 'Soil Water Content of the root zone'
+    model_outputs['SWCr'].attrs['description'] = 'Soil water content of the root zone in milimeters'
+    model_outputs['SWCr'].attrs['scale factor'] = '1'
+    
+    # Irrigation
+    model_outputs['Irr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['Irr'].attrs['units'] = 'mm'
+    model_outputs['Irr'].attrs['standard_name'] = 'Irrigation'
+    model_outputs['Irr'].attrs['description'] = 'Simulated daily irrigation in milimeters'
+    model_outputs['Irr'].attrs['scale factor'] = '1'
+    
+    # Deep Percolation
+    model_outputs['DP'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['DP'].attrs['units'] = 'mm'
+    model_outputs['DP'].attrs['standard_name'] = 'Deep Percolation'
+    model_outputs['DP'].attrs['description'] = 'Simulated daily Deep Percolation in milimeters'
+    model_outputs['DP'].attrs['scale factor'] = '1'
+    
+    if additional_outputs:
+        for var in additional_outputs:
+            model_outputs[var] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    
+    return model_outputs
+
+
+def calculate_diff_re(variable_ds: xr.Dataset, param_ds: xr.Dataset, scale_dict: dict, var: str) -> xr.DataArray:
+    """
+    Calculates the diffusion between the top soil layer and the root layer.
+
+    ## Arguments
+    1. variable_ds: `xr.Dataset`
+        dataset containing calculation variables
+    2. param_ds: `xr.Dataset`
+        dataset containing the rasterized parameters
+    3. scale_dict: `dict`
+        dictionnary containing the scale factors for
+        the rasterized parameters
+    4. var: `str`
+        name of depletion variable to use (Dei or Dep or De)
+
+    ## Returns
+    1. diff_re: `xr.Dataset`
+        the diffusion between the top soil layer and
+        the root layer
+    """
+    
+    # Temporary variables to make calculation easier to read
+    # tmp1 = (((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE'])
+    # tmp2 = ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr']
+    
+    # Calculate diffusion according to SAMIR equation
+    # diff_re = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))
+
+    # Return zero values where the 'DiffE' parameter is equal to 0
+    return xr.where(param_ds['DiffE'] == 0, 0, xr.where((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']) < 0, xr_maximum((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']), ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr']), xr_minimum((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']), ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr'])))
+
+
+def calculate_diff_dr(variable_ds: xr.Dataset, param_ds: xr.Dataset, scale_dict: dict) -> xr.DataArray:
+    """
+    Calculates the diffusion between the root layer and the deep layer.
+
+    ## Arguments
+    1. variable_ds: `xr.Dataset`
+        dataset containing calculation variables
+    2. param_ds: `xr.Dataset`
+        dataset containing the rasterized parameters
+    3. scale_dict: `dict`
+        dictionnary containing the scale factors for
+        the rasterized parameters
+
+    ## Returns
+    1. diff_dr: `xr.Dataset`
+        the diffusion between the root layer and the
+        deep layer
+    """
+    
+    # Temporary variables to make calculation easier to read
+    # tmp1 = (((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR']
+    # tmp2 = (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd']
+    
+    # Calculate diffusion according to SAMIR equation
+    # diff_dr = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))
     
-    return variables_t1, variables_t2, constant_values
+    # Return zero values where the 'DiffR' parameter is equal to 0
+    # return xr.where(param_ds['DiffR'] == 0, 0, diff_dr)
+    return xr.where(param_ds['DiffR'] == 0, 0, xr.where((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'] < 0, xr_maximum((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'], (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd']), xr_minimum((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'], (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd'])))
 
 
 def run_samir():

input/calculate_ndvi.py

@@ -97,10 +97,9 @@ def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, bounda
     ndvi['ndvi'] = (ndvi.ndvi*255).sortby('time')
     
     # Write attributes
-    ndvi['ndvi'].attrs['long_name'] = 'Normalized Difference Vegetation Index'
     ndvi['ndvi'].attrs['units'] = 'None'
     ndvi['ndvi'].attrs['standard_name'] = 'NDVI'
-    ndvi['ndvi'].attrs['comment'] = 'Normalized difference of the near infrared and red band. A value of one is a high vegetation presence.'
+    ndvi['ndvi'].attrs['description'] = 'Normalized difference Vegetation Index (of the near infrared and red band). A value of one is a high vegetation presence.'
     ndvi['ndvi'].attrs['scale factor'] = '255'
     
     # Create save path

input/lib_era5_land_pixel.py

@@ -429,7 +429,7 @@ def calculate_ET0_pixel(pixel_dataset: xr.Dataset, lat: float, lon: float, h: fl
     return ET0_values
 
 
-def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str, ndvi_path: str, h: float = 10) -> None:
+def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str, h: float = 10) -> None:
     """
     Calculate ET0 values from the ERA5 netcdf weather variables.
     Output netcdf contains the ET0 values for each day in the selected
@@ -491,4 +491,5 @@ def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str
 
     # Save dataset to netcdf, still in wgs84 (lat, lon) coordinates
     final_weather_ds.to_netcdf(path = output_nc_file, encoding = {"ET0": {"dtype": "u2"}, "tp": {"dtype": "u2"}})
+    
     return None
\ No newline at end of file

parameters/csv_files/params_samir_test.csv

+ClassName,scale_factor,Default,class1,class2,class3,class4,class5,class6,class7
+ClassNumber,1,0,1,2,3,4,5,6,7
+FmaxFC,1000,1,0.9,1,1,1,1,,1
+Fslope,1000,1.4,1.5,1.5,1.5,1.2,1.5,0,1.5
+Foffset,1000,-0.1,-0.1,-0.1,-0.1,-0.17,-0.1,0,-0.1
+Plateau,1,70,,,,,,,
+KmaxKcb,1000,0.98,1,1.1,1.1,0.65,1.13,,1.13
+Kslope,1000,1.6,1.6,1.6,1.6,1.2,1.6,0,1.6
+Koffset,1000,-0.1,-0.1,-0.1,-0.1,-0.18,-0.1,0,-0.1
+Zsoil,1,2000,1600,1550,1550,2000,1550,2000,1550
+Ze,1,125,100,125,125,100,100,125,100
+Init_RU,1,0,,,,,,,
+DiffE,1,5,5,5,5,5,,,
+DiffR,1,10,10,10,10,10,,,
+REW,1,0,,,,,,,
+minZr,1,150,1550,125,125,1550,125,125,125
+maxZr,1,1000,1550,800,800,1550,1000,126,1000
+p,1000,0,0.5,0.55,0.55,0.65,0.7,0.45,0.7
+FW,1,100,30,100,100,30,100,30,100
+Irrig_auto,1,0,,,,,,0,
+Irrig_man,1,0,,,,,,,
+Lame_max,1,0,,,,,,,
+Kcmax,1000,1.15,,,,,,,
+Fc_stop,1000,0.95,,,,,,,
\ No newline at end of file