diff --git a/DEV_inputs_test/land_cover.nc b/DEV_inputs_test/land_cover.nc
new file mode 100644
index 0000000000000000000000000000000000000000..381acc9d9e01422e5725274e38e090a7357e6b3c
Binary files /dev/null and b/DEV_inputs_test/land_cover.nc differ
diff --git a/DEV_inputs_test/ndvi.nc b/DEV_inputs_test/ndvi.nc
new file mode 100644
index 0000000000000000000000000000000000000000..28e2e5c50cfebff2b904645c8b5fc476bc0d1e50
Binary files /dev/null and b/DEV_inputs_test/ndvi.nc differ
diff --git a/DEV_inputs_test/soil.nc b/DEV_inputs_test/soil.nc
new file mode 100644
index 0000000000000000000000000000000000000000..005a8d4325786221abf7ea9b1890480682e2d1e8
Binary files /dev/null and b/DEV_inputs_test/soil.nc differ
diff --git a/DEV_inputs_test/weather.nc b/DEV_inputs_test/weather.nc
new file mode 100644
index 0000000000000000000000000000000000000000..0d7f57fe048eeab4f18a23da0b727ed72d6e8ceb
Binary files /dev/null and b/DEV_inputs_test/weather.nc differ
diff --git a/code/modspa_samir.py b/code/modspa_samir.py
index cce94e486f6319e62251b799f5a1b7337cea6e82..39575984bcd9099a595d15176d8dc1f3d513b21e 100644
--- a/code/modspa_samir.py
+++ b/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():
diff --git a/dev_samir_xarray.ipynb b/dev_samir_xarray.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..99173abc4a5c9dfdb74056f58e23f66f8d95b61a
--- /dev/null
+++ b/dev_samir_xarray.ipynb
@@ -0,0 +1,1257 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import xarray as xr\n",
+ "from dask.distributed import Client\n",
+ "import os, sys\n",
+ "import pandas as pd\n",
+ "import numpy as np\n",
+ "from typing import List, Tuple, Union\n",
+ "import warnings\n",
+ "import gc\n",
+ "\n",
+ "from parameters.params_samir_class import samir_parameters\n",
+ "from config.config import config\n",
+ "ndvi_path = '/mnt/e/DATA/DEV_inputs_test/ndvi.nc'\n",
+ "weather_path = '/mnt/e/DATA/DEV_inputs_test/weather.nc'\n",
+ "land_cover_path = '/mnt/e/DATA/DEV_inputs_test/land_cover.nc'\n",
+ "\n",
+ "param_file = '/home/auclairj/GIT/modspa-pixel/parameters/csv_files/params_samir_test.csv'"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "\n",
+ "def rasterize_samir_parameters(csv_param_file: str, empty_dataset: xr.Dataset, land_cover_raster: str, chunk_size: dict) -> Tuple[xr.Dataset, dict]:\n",
+ " \"\"\"\n",
+ " Creates a raster `xarray` dataset from the csv parameter file, the land cover raster and an empty dataset\n",
+ " that contains the right structure (emptied ndvi dataset for example). For each parameter, the function loops\n",
+ " on land cover classes to fill the raster.\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. csv_param_file: `str`\n",
+ " path to csv paramter file\n",
+ " 2. empty_dataset: `xr.Dataset`\n",
+ " empty dataset that contains the right structure (emptied ndvi dataset for example).\n",
+ " 3. land_cover_raster: `str`\n",
+ " path to land cover netcdf raster\n",
+ " 4. chunk_size: `dict`\n",
+ " chunk_size for dask computation\n",
+ "\n",
+ " ## Returns\n",
+ " 1. parameter_dataset: `xr.Dataset`\n",
+ " the dataset containing all the rasterized Parameters\n",
+ " 2. scale_factor: `dict`\n",
+ " dictionnary containing the scale factors for each parameter\n",
+ " \"\"\"\n",
+ " \n",
+ " # Load samir params into an object\n",
+ " table_param = samir_parameters(csv_param_file)\n",
+ " \n",
+ " # Set general variables\n",
+ " class_count = table_param.table.shape[1] - 2 # remove dtype and default columns\n",
+ " \n",
+ " # Open land cover raster\n",
+ " land_cover = xr.open_dataarray(land_cover_raster, chunks = chunk_size)\n",
+ " \n",
+ " # Create dataset\n",
+ " parameter_dataset = empty_dataset.copy(deep = True)\n",
+ " \n",
+ " # Create dictionnary containing the scale factors\n",
+ " scale_factor = {}\n",
+ " \n",
+ " # Loop on samir parameters and create \n",
+ " for parameter in table_param.table.index[1:]:\n",
+ " \n",
+ " # Create new variable and set attributes\n",
+ " parameter_dataset[parameter] = land_cover.copy(deep = True).astype('f4')\n",
+ " parameter_dataset[parameter].attrs['name'] = parameter\n",
+ " parameter_dataset[parameter].attrs['description'] = 'cf SAMIR Doc for detail'\n",
+ " parameter_dataset[parameter].attrs['scale factor'] = str(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])\n",
+ " \n",
+ " # Assigne value in dictionnary\n",
+ " scale_factor[parameter] = 1/int(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])\n",
+ " \n",
+ " # Loop on classes to set parameter values for each class\n",
+ " for class_val, class_name in zip(range(1, class_count + 1), table_param.table.columns[2:]):\n",
+ " \n",
+ " # Parameter values are multiplied by the scale factor in order to store all values as int16 types\n",
+ " # These values are then rounded to make sure there isn't any decimal point issues when casting the values to int16\n",
+ " 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('f4')\n",
+ " \n",
+ " # Return dataset converted to 'int16' data type to reduce memory usage\n",
+ " # and scale_factor dictionnary for later conversion\n",
+ " return parameter_dataset, scale_factor\n",
+ "\n",
+ "\n",
+ "def setup_time_loop(calculation_variables_t1: List[str], calculation_variables_t2: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:\n",
+ " \"\"\"\n",
+ " Creates two temporary `xarray Datasets` that will be used in the SAMIR time loop.\n",
+ " `variables_t1` corresponds to the variables for the previous day and `variables_t2`\n",
+ " corresponds to the variables for the current day. After each loop, `variables_t1`\n",
+ " takes the value of `variables_t2` for the corresponding variables.\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. calculation_variables_t1: `List[str]`\n",
+ " list of strings containing the variable names\n",
+ " for the previous day dataset\n",
+ " 2. calculation_variables_t2: `List[str]`\n",
+ " list of strings containing the variable names\n",
+ " for the current day dataset\n",
+ " 3. empty_dataset: `xr.Dataset`\n",
+ " empty dataset that contains the right structure\n",
+ "\n",
+ " ## Returns\n",
+ " 1. variables_t1: `xr.Dataset`\n",
+ " output dataset for previous day\n",
+ " 2. variables_t2: `xr.Dataset`\n",
+ " output dataset for current day\n",
+ " \"\"\"\n",
+ " \n",
+ " # Create new dataset\n",
+ " variables_t1 = empty_dataset.copy(deep = True)\n",
+ " \n",
+ " # Create empty DataArray for each variable\n",
+ " for variable in calculation_variables_t1:\n",
+ " \n",
+ " # Assign new empty DataArray\n",
+ " variables_t1[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))\n",
+ " variables_t1[variable].attrs['name'] = variable # set name in attributes\n",
+ " \n",
+ " # Create new dataset\n",
+ " variables_t2 = empty_dataset.copy(deep = True)\n",
+ " \n",
+ " # Create empty DataArray for each variable\n",
+ " for variable in calculation_variables_t2:\n",
+ " \n",
+ " # Assign new empty DataArray\n",
+ " variables_t2[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))\n",
+ " variables_t2[variable].attrs['name'] = variable # set name in attributes\n",
+ " \n",
+ " return variables_t1, variables_t2\n",
+ "\n",
+ "\n",
+ "def prepare_outputs(empty_dataset: xr.Dataset, additional_outputs: List[str] = None) -> xr.Dataset:\n",
+ " \"\"\"\n",
+ " Creates the `xarray Dataset` containing the outputs of the SAMIR model that will be saved.\n",
+ " Additional variables can be saved by adding their names to the `additional_outputs` list.\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. empty_dataset: `xr.Dataset`\n",
+ " empty dataset that contains the right structure\n",
+ " 2. additional_outputs: `List[str]`\n",
+ " list of additional variable names to be saved\n",
+ "\n",
+ " ## Returns\n",
+ " 1. model_outputs: `xr.Dataset`\n",
+ " model outputs to be saved\n",
+ " \"\"\"\n",
+ " \n",
+ " # Evaporation and Transpiraion\n",
+ " model_outputs = empty_dataset.copy(deep = True)\n",
+ " \n",
+ " model_outputs['E'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['E'].attrs['units'] = 'mm'\n",
+ " model_outputs['E'].attrs['standard_name'] = 'Evaporation'\n",
+ " model_outputs['E'].attrs['description'] = 'Accumulated daily evaporation in milimeters'\n",
+ " model_outputs['E'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " model_outputs['Tr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['Tr'].attrs['units'] = 'mm'\n",
+ " model_outputs['Tr'].attrs['standard_name'] = 'Transpiration'\n",
+ " model_outputs['Tr'].attrs['description'] = 'Accumulated daily plant transpiration in milimeters'\n",
+ " model_outputs['Tr'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " # Soil Water Content\n",
+ " model_outputs['SWCe'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['SWCe'].attrs['units'] = 'mm'\n",
+ " model_outputs['SWCe'].attrs['standard_name'] = 'Soil Water Content of the evaporative zone'\n",
+ " model_outputs['SWCe'].attrs['description'] = 'Soil water content of the evaporative zone in milimeters'\n",
+ " model_outputs['SWCe'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " model_outputs['SWCr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['SWCr'].attrs['units'] = 'mm'\n",
+ " model_outputs['SWCr'].attrs['standard_name'] = 'Soil Water Content of the root zone'\n",
+ " model_outputs['SWCr'].attrs['description'] = 'Soil water content of the root zone in milimeters'\n",
+ " model_outputs['SWCr'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " # Irrigation\n",
+ " model_outputs['Irr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['Irr'].attrs['units'] = 'mm'\n",
+ " model_outputs['Irr'].attrs['standard_name'] = 'Irrigation'\n",
+ " model_outputs['Irr'].attrs['description'] = 'Simulated daily irrigation in milimeters'\n",
+ " model_outputs['Irr'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " # Deep Percolation\n",
+ " model_outputs['DP'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " model_outputs['DP'].attrs['units'] = 'mm'\n",
+ " model_outputs['DP'].attrs['standard_name'] = 'Deep Percolation'\n",
+ " model_outputs['DP'].attrs['description'] = 'Simulated daily Deep Percolation in milimeters'\n",
+ " model_outputs['DP'].attrs['scale factor'] = '1'\n",
+ " \n",
+ " if additional_outputs:\n",
+ " for var in additional_outputs:\n",
+ " model_outputs[var] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+ " \n",
+ " return model_outputs\n",
+ "\n",
+ "\n",
+ "def xr_maximum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Equivalent of `numpy.maximum(ds, value)` for xarray DataArrays\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. ds: `xr.DataArray`\n",
+ " datarray to compare\n",
+ " 2. value: `Union[xr.DataArray, float, int]`\n",
+ " value (scalar or dataarray) to compare\n",
+ "\n",
+ " ## Returns\n",
+ " 1. output: `xr.DataArray`\n",
+ " resulting dataarray with maximum value element-wise\n",
+ " \"\"\"\n",
+ " return xr.where(ds <= value, value, ds, keep_attrs = True)\n",
+ "\n",
+ "\n",
+ "def xr_minimum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Equivalent of `numpy.minimum(ds, value)` for xarray DataArrays\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. ds: `xr.DataArray`\n",
+ " datarray to compare\n",
+ " 2. value: `Union[xr.DataArray, float, int]`\n",
+ " value (scalar or dataarray) to compare\n",
+ "\n",
+ " ## Returns\n",
+ " 1. output: `xr.DataArray`\n",
+ " resulting dataarray with minimum value element-wise\n",
+ " \"\"\"\n",
+ " return xr.where(ds >= value, value, ds, keep_attrs = True)\n",
+ "\n",
+ "\n",
+ "def calculate_diff_re(TAW: xr.DataArray, Dr: xr.DataArray, Zr: xr.DataArray, RUE: xr.DataArray, De: xr.DataArray, FCov: xr.DataArray, Ze_: xr.DataArray, DiffE_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Calculates the diffusion between the top soil layer and the root layer.\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TAW: `xr.DataArray`\n",
+ " water capacity of root layer\n",
+ " 2. Dr: `xr.DataArray`\n",
+ " depletion of root layer\n",
+ " 3. Zr: `xr.DataArray`\n",
+ " height of root layer\n",
+ " 4. RUE: `xr.DataArray`\n",
+ " total available surface water\n",
+ " 5. De: `xr.DataArray`\n",
+ " depletion of the evaporative layer\n",
+ " 6. FCov: `xr.DataArray`\n",
+ " fraction cover of plants\n",
+ " 7. Ze_: `xr.DataArray`\n",
+ " height of evaporative layer (paramter)\n",
+ " 8. DiffE_: `xr.DataArray`\n",
+ " diffusion coefficient between evaporative\n",
+ " and root layers (unitless, parameter)\n",
+ " 9. scale_dict: `dict`\n",
+ " dictionnary containing the scale factors for\n",
+ " the rasterized parameters\n",
+ "\n",
+ " ## Returns\n",
+ " 1. diff_re: `xr.Dataset`\n",
+ " the diffusion between the top soil layer and\n",
+ " the root layer\n",
+ " \"\"\"\n",
+ " \n",
+ " # Temporary variables to make calculation easier to read\n",
+ " tmp1 = (((TAW - Dr) / Zr - (RUE - De) / (scale_dict['Ze'] * Ze_)) / FCov) * (scale_dict['DiffE'] * DiffE_)\n",
+ " tmp2 = ((TAW * scale_dict['Ze'] * Ze_) - (RUE - De - Dr) * Zr) / (Zr + scale_dict['Ze'] * Ze_) - Dr\n",
+ " \n",
+ " # Calculate diffusion according to SAMIR equation\n",
+ " diff_re = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))\n",
+ "\n",
+ " # Return zero values where the 'DiffE' parameter is equal to 0\n",
+ " return xr.where(DiffE_ == 0, 0, diff_re)\n",
+ "\n",
+ "\n",
+ "def calculate_diff_dr(TAW: xr.DataArray, TDW: xr.DataArray, Dr: xr.DataArray, Zr: xr.DataArray, Dd: xr.DataArray, FCov: xr.DataArray, Zsoil_: xr.DataArray, DiffR_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Calculates the diffusion between the root layer and the deep layer.\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TAW: `xr.DataArray`\n",
+ " water capacity of root layer\n",
+ " 2. TDW: `xr.DataArray`\n",
+ " water capacity of deep layer\n",
+ " 3. Dr: `xr.DataArray`\n",
+ " depletion of root layer\n",
+ " 4. Zr: `xr.DataArray`\n",
+ " height of root layer\n",
+ " 5. Dd: `xr.DataArray`\n",
+ " depletion of deep layer\n",
+ " 6. FCov: `xr.DataArray`\n",
+ " fraction cover of plants\n",
+ " 7. Zsoil_: `xr.DataArray`\n",
+ " total height of soil (paramter)\n",
+ " 8. DiffR_: `xr.DataArray`\n",
+ " Diffusion coefficient between root\n",
+ " and deep layers (unitless, parameter)\n",
+ " 9. scale_dict: `dict`\n",
+ " dictionnary containing the scale factors for\n",
+ " the rasterized parameters\n",
+ "\n",
+ " ## Returns\n",
+ " 1. diff_dr: `xr.Dataset`\n",
+ " the diffusion between the root layer and the\n",
+ " deep layer\n",
+ " \"\"\"\n",
+ " \n",
+ " # Temporary variables to make calculation easier to read\n",
+ " tmp1 = (((TDW - Dd) / (scale_dict['Zsoil'] * Zsoil_ - Zr) - (TAW - Dr) / Zr) / FCov) * scale_dict['DiffR'] * DiffR_\n",
+ " tmp2 = (TDW *Zr - (TAW - Dr - Dd) * (scale_dict['Zsoil'] * Zsoil_ - Zr)) / (scale_dict['Zsoil'] * Zsoil_) - Dd\n",
+ " \n",
+ " # Calculate diffusion according to SAMIR equation\n",
+ " diff_dr = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))\n",
+ " \n",
+ " # Return zero values where the 'DiffR' parameter is equal to 0\n",
+ " return xr.where(DiffR_ == 0, 0, diff_dr)\n",
+ "\n",
+ "\n",
+ "def calculate_W(TEW: xr.DataArray, Dei: xr.DataArray, Dep: xr.DataArray, fewi: xr.DataArray, fewp: xr.DataArray) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Calculate W, the weighting factor to split the energy available\n",
+ " for evaporation depending on the difference in water availability\n",
+ " in the two evaporation components, ensuring that the larger and\n",
+ " the wetter, the more the evaporation occurs from that component\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TEW: `xr.DataArray`\n",
+ " water capacity of evaporative layer\n",
+ " 2. Dei: `xr.DataArray`\n",
+ " depletion of the evaporative layer\n",
+ " (irrigation part)\n",
+ " 3. Dep: `xr.DataArray`\n",
+ " depletion of the evaporative layer\n",
+ " (precipitation part)\n",
+ " 4. fewi: `xr.DataArray`\n",
+ " soil fraction which is wetted by irrigation\n",
+ " and exposed to evaporation\n",
+ " 5. fewp: `xr.DataArray`\n",
+ " soil fraction which is wetted by precipitation\n",
+ " and exposed to evaporation\n",
+ "\n",
+ " ## Returns\n",
+ " 1. W: `xr.DataArray`\n",
+ " weighting factor W\n",
+ " \"\"\"\n",
+ " \n",
+ " # Temporary variables to make calculation easier to read\n",
+ " tmp = fewi * (TEW - Dei)\n",
+ " \n",
+ " # Calculate the weighting factor to split the energy available for evaporation\n",
+ " W = 1 / (1 + (fewp * (TEW - Dep) / tmp ))\n",
+ "\n",
+ " # Return W \n",
+ " return xr.where(tmp > 0, W, 0)\n",
+ "\n",
+ "\n",
+ "def calculate_Kr(TEW: xr.DataArray, De: xr.DataArray, REW_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " calculates of the reduction coefficient for evaporation dependent \n",
+ " on the amount of water in the soil using the FAO-56 method\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TEW: `xr.DataArray`\n",
+ " water capacity of evaporative layer\n",
+ " 2. De: `xr.DataArray`\n",
+ " depletion of evaporative layer\n",
+ " 3. REW_: `xr.DataArray`\n",
+ " readily evaporable water\n",
+ " 4. scale_dict: `dict`\n",
+ " dictionnary containing the scale factors for\n",
+ " the rasterized parameters\n",
+ "\n",
+ " ## Returns\n",
+ " 1. Kr: `xr.DataArray`\n",
+ " Kr coefficient\n",
+ " \"\"\"\n",
+ " \n",
+ " # Formula for calculating Kr\n",
+ " Kr = (TEW - De) / (TEW - scale_dict['REW'] * REW_)\n",
+ " \n",
+ " # Return Kr\n",
+ " return xr_maximum(0, xr_minimum(Kr, 1))\n",
+ "\n",
+ "\n",
+ "def update_Dr(TAW: xr.DataArray, TDW: xr.DataArray, Zr: xr.DataArray, TAW0: xr.DataArray, TDW0: xr.DataArray, Dr0: xr.DataArray, Dd0: xr.DataArray, Zr0: xr.DataArray) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Return the updated depletion for the root layer\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TAW: `xr.DataArray`\n",
+ " water capacity of root layer for current day\n",
+ " 2. TDW: `xr.DataArray`\n",
+ " water capacity of deep layer for current day\n",
+ " 3. Zr: `xr.DataArray`\n",
+ " root layer height for current day\n",
+ " 4. TAW0: `xr.DataArray`\n",
+ " water capacity of root layer for previous day\n",
+ " 5. TDW0: `xr.DataArray`\n",
+ " water capacity of deep layer for previous day\n",
+ " 6. Dr0: `xr.DataArray`\n",
+ " depletion of the root layer for previous day\n",
+ " 7. Dd0: `xr.DataArray`\n",
+ " depletion of the deep laye for previous day\n",
+ " 8. Zr0: `xr.DataArray`\n",
+ " root layer height for previous day\n",
+ "\n",
+ " ## Returns\n",
+ " 1. output: `xr.DataArray`\n",
+ " updated depletion for the root layer\n",
+ " \"\"\"\n",
+ " \n",
+ " # Temporary variables to make calculation easier to read\n",
+ " tmp1 = xr_maximum(Dr0 + Dd0 * (TAW - TAW0) / TDW0, 0)\n",
+ " tmp2 = xr_minimum(Dr0 + Dd0 * (TAW - TAW0) / TDW0, TDW)\n",
+ "\n",
+ " # Return updated Dr\n",
+ " return xr.where(Zr > Zr0, tmp1, tmp2)\n",
+ "\n",
+ "\n",
+ "def update_Dd(TAW: xr.DataArray, TDW: xr.DataArray, Zr: xr.DataArray, TAW0: xr.DataArray, TDW0: xr.DataArray, Dd0: xr.DataArray, Zr0: xr.DataArray) -> xr.DataArray:\n",
+ " \"\"\"\n",
+ " Return the updated depletion for the deep layer\n",
+ "\n",
+ " ## Arguments\n",
+ " 1. TAW: `xr.DataArray`\n",
+ " water capacity of root layer for current day\n",
+ " 2. TDW: `xr.DataArray`\n",
+ " water capacity of deep layer for current day\n",
+ " 3. TAW0: `xr.DataArray`\n",
+ " water capacity of root layer for previous day\n",
+ " 5. TDW0: `xr.DataArray`\n",
+ " water capacity of deep layer for previous day\n",
+ " 6. Dd0: `xr.DataArray`\n",
+ " depletion of the deep laye for previous day\n",
+ " 7. Zr0: `xr.DataArray`\n",
+ " root layer height for previous day\n",
+ "\n",
+ " ## Returns\n",
+ " 1. output: `xr.DataArray`\n",
+ " updated depletion for the deep layer\n",
+ " \"\"\"\n",
+ " \n",
+ " # Temporary variables to make calculation easier to read\n",
+ " tmp1 = xr_maximum(Dd0 - Dd0 * (TAW - TAW0) / TDW0, 0)\n",
+ " tmp2 = xr_minimum(Dd0 - Dd0 * (TAW - TAW0) / TDW0, TDW)\n",
+ " \n",
+ " # Return updated Dd\n",
+ " return xr.where(Zr > Zr0, tmp1, tmp2)\n",
+ "\n",
+ "\n",
+ "def run_samir(json_config_file: str, csv_param_file: str, ndvi_cube_path: str, weather_cube_path: str, soil_params_path: str, land_cover_path: str, chunk_size: dict, save_path: str) -> None:\n",
+ " \n",
+ " # warnings.simplefilter(\"error\", category = RuntimeWarning())\n",
+ " warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in cast\")\n",
+ " warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in divide\")\n",
+ " np.errstate(all = 'raise')\n",
+ " gc.disable()\n",
+ " \n",
+ " #============ General parameters ============#\n",
+ " config_params = config(json_config_file)\n",
+ " calculation_variables_t2 = ['diff_rei', 'diff_rep', 'diff_dr' , 'Dd', 'De', 'Dei', 'Dep', 'Dr', 'FCov', 'Irrig', 'Kcb', 'Kei', 'Kep', 'Ks', 'Kti', 'Ktp', 'RUE', 'TAW', 'TDW', 'TEW', 'Tei', 'Tep', 'W', 'Zr', 'fewi', 'fewp']\n",
+ " calculation_variables_t1 = ['Dr', 'Dd', 'TAW', 'TDW', 'Zr']\n",
+ " \n",
+ " #============ Manage inputs ============#\n",
+ " # NDVI\n",
+ " ndvi_cube = xr.open_dataset(ndvi_cube_path, chunks = chunk_size).astype('f4')\n",
+ " \n",
+ " # # Create a daily DateTimeIndex for the desired date range\n",
+ " # daily_index = pd.date_range(start = config_params.start_date, end = config_params.end_date, freq = 'D')\n",
+ "\n",
+ " # # Resample the dataset to a daily frequency and reindex with the new DateTimeIndex\n",
+ " # ndvi_cube = ndvi_cube.resample(time = '1D').asfreq().reindex(time = daily_index)\n",
+ "\n",
+ " # # Interpolate the dataset along the time dimension to fill nan values\n",
+ " # ndvi_cube = ndvi_cube.interpolate_na(dim = 'time', method = 'linear', fill_value = 'extrapolate').astype('u1')\n",
+ " \n",
+ " # Weather\n",
+ " weather_cube = xr.open_dataset(weather_cube_path, chunks = chunk_size).astype('f4')\n",
+ " \n",
+ " # Soil\n",
+ " soil_params = xr.open_dataset(soil_params_path, chunks = chunk_size).astype('f4')\n",
+ " \n",
+ " # SAMIR Parameters\n",
+ " param_dataset, scale_factor = rasterize_samir_parameters(csv_param_file, ndvi_cube.drop_vars(['ndvi', 'time']), land_cover_path, chunk_size = chunk_size)\n",
+ " \n",
+ " # SAMIR Variables\n",
+ " variables_t1, variables_t2 = setup_time_loop(calculation_variables_t1, calculation_variables_t2, ndvi_cube.drop_vars(['ndvi', 'time']))\n",
+ "\n",
+ " #============ Prepare outputs ============#\n",
+ " model_outputs = prepare_outputs(ndvi_cube.drop_vars(['ndvi']))\n",
+ " \n",
+ " #============ Prepare time iterations ============#\n",
+ " dates = ndvi_cube.time.values\n",
+ " \n",
+ " #============ Create aliases for better readability ============#\n",
+ " \n",
+ " # Variables for current day\n",
+ " diff_rei = variables_t2['diff_rei']\n",
+ " diff_rep = variables_t2['diff_rep']\n",
+ " diff_dr = variables_t2['diff_dr']\n",
+ " Dd = variables_t2['Dd']\n",
+ " De = variables_t2['De']\n",
+ " Dei = variables_t2['Dei']\n",
+ " Dep = variables_t2['Dep']\n",
+ " Dr = variables_t2['Dr']\n",
+ " FCov = variables_t2['FCov']\n",
+ " Irrig = variables_t2['Irrig']\n",
+ " Kcb = variables_t2['Kcb']\n",
+ " Kei = variables_t2['Kei']\n",
+ " Kep = variables_t2['Kep']\n",
+ " Ks = variables_t2['Ks']\n",
+ " Kti = variables_t2['Kti']\n",
+ " Ktp = variables_t2['Ktp']\n",
+ " RUE = variables_t2['RUE']\n",
+ " TAW = variables_t2['TAW']\n",
+ " TDW = variables_t2['TDW']\n",
+ " TEW = variables_t2['TEW']\n",
+ " Tei = variables_t2['Tei']\n",
+ " Tep = variables_t2['Tep']\n",
+ " Zr = variables_t2['Zr']\n",
+ " W = variables_t2['W']\n",
+ " fewi = variables_t2['fewi']\n",
+ " fewp = variables_t2['fewp']\n",
+ " \n",
+ " # Variables for previous day\n",
+ " TAW0 = variables_t1['TAW']\n",
+ " TDW0 = variables_t1['TDW']\n",
+ " Dr0 = variables_t1['Dr']\n",
+ " Dd0 = variables_t1['Dd']\n",
+ " Zr0 = variables_t1['Zr']\n",
+ " \n",
+ " # Parameters\n",
+ " # Parameters have an underscore (_) behind their name for recognition \n",
+ " DiffE_ = param_dataset['DiffE']\n",
+ " DiffR_ = param_dataset['DiffR']\n",
+ " FW_ = param_dataset['FW']\n",
+ " Fc_stop_ = param_dataset['Fc_stop']\n",
+ " FmaxFC_ = param_dataset['FmaxFC']\n",
+ " Foffset_ = param_dataset['Foffset']\n",
+ " Fslope_ = param_dataset['Fslope']\n",
+ " Init_RU_ = param_dataset['Init_RU']\n",
+ " Irrig_auto_ = param_dataset['Irrig_auto']\n",
+ " Kcmax_ = param_dataset['Kcmax']\n",
+ " KmaxKcb_ = param_dataset['KmaxKcb']\n",
+ " Koffset_ = param_dataset['Koffset']\n",
+ " Kslope_ = param_dataset['Kslope']\n",
+ " Lame_max_ = param_dataset['Lame_max']\n",
+ " REW_ = param_dataset['REW']\n",
+ " Ze_ = param_dataset['Ze']\n",
+ " Zsoil_ = param_dataset['Zsoil']\n",
+ " maxZr_ = param_dataset['maxZr']\n",
+ " minZr_ = param_dataset['minZr']\n",
+ " p_ = param_dataset['p']\n",
+ " \n",
+ " # scale factors\n",
+ " # Scale factors have the following name scheme : s_ + parameter_name\n",
+ " s_DiffE = scale_factor['DiffE']\n",
+ " s_DiffR = scale_factor['DiffR']\n",
+ " s_FW = scale_factor['FW']\n",
+ " s_Fc_stop = scale_factor['Fc_stop']\n",
+ " s_FmaxFC = scale_factor['FmaxFC']\n",
+ " s_Foffset = scale_factor['Foffset']\n",
+ " s_Fslope = scale_factor['Fslope']\n",
+ " s_Init_RU = scale_factor['Init_RU']\n",
+ " # s_Irrig_auto = scale_factor['Irrig_auto']\n",
+ " s_Kcmax = scale_factor['Kcmax']\n",
+ " s_KmaxKcb = scale_factor['KmaxKcb']\n",
+ " s_Koffset = scale_factor['Koffset']\n",
+ " s_Kslope = scale_factor['Kslope']\n",
+ " s_Lame_max = scale_factor['Lame_max']\n",
+ " s_REW = scale_factor['REW']\n",
+ " s_Ze = scale_factor['Ze']\n",
+ " s_Zsoil = scale_factor['Zsoil']\n",
+ " s_maxZr = scale_factor['maxZr']\n",
+ " s_minZr = scale_factor['minZr']\n",
+ " s_p = scale_factor['p']\n",
+ " \n",
+ " #============ First day initialization ============#\n",
+ " # Fraction cover\n",
+ " FCov = s_Fslope * Fslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]})/255) + s_Foffset * Foffset_\n",
+ " FCov = xr_minimum(xr_maximum(FCov, 0), s_Fc_stop * Fc_stop_)\n",
+ " \n",
+ " # Root depth upate\n",
+ " Zr = s_minZr * minZr_ + (FCov / (s_FmaxFC * FmaxFC_)) * s_maxZr * (maxZr_ - minZr_)\n",
+ " \n",
+ " # Water capacities\n",
+ " TEW = (soil_params['FC'] - soil_params['WP']/2) * s_Ze * Ze_\n",
+ " RUE = (soil_params['FC'] - soil_params['WP']) * s_Ze * Ze_\n",
+ " TAW = (soil_params['FC'] - soil_params['WP']) * Zr\n",
+ " TDW = (soil_params['FC'] - soil_params['WP']) * (s_Zsoil * Zsoil_ - Zr) # Zd = Zsoil - Zr\n",
+ " \n",
+ " # Depletions\n",
+ " Dei = RUE * (1 - s_Init_RU * Init_RU_)\n",
+ " Dep = RUE * (1 - s_Init_RU * Init_RU_)\n",
+ " Dr = TAW * (1 - s_Init_RU * Init_RU_)\n",
+ " Dd = TDW * (1 - s_Init_RU * Init_RU_)\n",
+ " \n",
+ " # Irrigation ==============!!!!!\n",
+ " Irrig = xr_minimum(xr_maximum(Dr - weather_cube['tp'].sel({'time': dates[0]}) / 1000, 0), s_Lame_max * Lame_max_) * Irrig_auto_\n",
+ " Irrig = xr.where(Dr > TAW * s_p * p_, Irrig, 0)\n",
+ " \n",
+ " # Kcb\n",
+ " Kcb = xr_minimum(s_Kslope * Kslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]}) / 255) + s_Koffset * Koffset_, s_KmaxKcb * KmaxKcb_)\n",
+ " \n",
+ " # Update depletions with rainfall and/or irrigation \n",
+ " ## DP\n",
+ " model_outputs['DP'].loc[{'time': dates[0]}] = -xr_minimum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), 0)\n",
+ " \n",
+ " ## De\n",
+ " Dei = xr_minimum(xr_maximum(Dei - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig / (s_FW * FW_ / 100), 0), TEW)\n",
+ " Dep = xr_minimum(xr_maximum(Dep - (weather_cube['tp'].sel({'time': dates[0]}) / 1000), 0), TEW)\n",
+ " \n",
+ " fewi = xr_minimum(1 - FCov, (s_FW * FW_ / 100))\n",
+ " fewp = 1 - FCov - fewi\n",
+ " \n",
+ " De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+ " # variables_t1['De'] = xr.where(variables_t1['De'].isfinite(), variables_t1['De'], variables_t1['Dei'] * (s_FW * FW_ / 100) + variables_t1['Dep'] * (1 - (s_FW * FW_ / 100))) #================= find replacement for .isfinite() method !!!!!!!!!\n",
+ "\n",
+ " ## Dr\n",
+ " Dr = xr_minimum(xr_maximum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), TAW)\n",
+ " \n",
+ " ## Dd\n",
+ " Dd = xr_minimum(xr_maximum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), 0), TDW)\n",
+ " \n",
+ " # Diffusion coefficients\n",
+ " diff_rei = calculate_diff_re(TAW, Dr, Zr, RUE, Dei, FCov, Ze_, DiffE_, scale_factor)\n",
+ " diff_rep = calculate_diff_re(TAW, Dr, Zr, RUE, Dep, FCov, Ze_, DiffE_, scale_factor)\n",
+ " diff_dr = calculate_diff_dr(TAW, TDW, Dr, Zr, Dd, FCov, Zsoil_, DiffR_, scale_factor) \n",
+ " \n",
+ " # Weighing factor W\n",
+ " W = calculate_W(TEW, Dei, Dep, fewi, fewp)\n",
+ " \n",
+ " # Soil water contents\n",
+ " model_outputs['SWCe'].loc[{'time': dates[0]}] = 1 - De/TEW\n",
+ " model_outputs['SWCr'].loc[{'time': dates[0]}] = 1 - Dr/TAW\n",
+ " \n",
+ " # Water Stress coefficient\n",
+ " Ks = xr_minimum((TAW - Dr) / (TAW * (1 - s_p * p_)), 1)\n",
+ " \n",
+ " # Reduction coefficient for evaporation\n",
+ " Kei = xr_minimum(W * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewi * s_Kcmax * Kcmax_)\n",
+ " Kep = xr_minimum((1 - W) * calculate_Kr(TEW, Dep, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewp * s_Kcmax * Kcmax_)\n",
+ " \n",
+ " # Prepare coefficients for evapotranspiration\n",
+ " Kti = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dei / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+ " Ktp = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dep / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+ " Tei = Kti * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+ " Tep = Ktp * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+ " \n",
+ " # Update depletions\n",
+ " Dei = xr.where(fewi > 0, xr_minimum(xr_maximum(Dei + (weather_cube['ET0'].sel({'time': dates[0]}) / 1000) * Kei / fewi + Tei - diff_rei, 0), TEW), xr_minimum(xr_maximum(Dei + Tei - diff_rei, 0), TEW))\n",
+ " Dep = xr.where(fewp > 0, xr_minimum(xr_maximum(Dep + (weather_cube['ET0'].sel({'time': dates[0]}) / 1000) * Kep / fewp + Tep - diff_rep, 0), TEW), xr_minimum(xr_maximum(Dep + Tep - diff_rep, 0), TEW))\n",
+ " \n",
+ " De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+ " # De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100))) #================= find replacement for .isfinite() method !!!!!!!!!\n",
+ " \n",
+ " # Evaporation\n",
+ " model_outputs['E'].loc[{'time': dates[0]}] = xr_maximum((Kei + Kep) * weather_cube['ET0'].sel({'time': dates[0]}) / 1000, 0)\n",
+ " \n",
+ " # Transpiration\n",
+ " model_outputs['Tr'].loc[{'time': dates[0]}] = Kcb * Ks * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+ " \n",
+ " # Potential evapotranspiration and evaporative fraction ??\n",
+ " \n",
+ " # Update depletions (root and deep zones) at the end of the day\n",
+ " Dr = xr_minimum(xr_maximum(Dr + model_outputs['E'].loc[{'time': dates[0]}] + model_outputs['Tr'].loc[{'time': dates[0]}] - diff_dr, 0), TAW)\n",
+ " Dd = xr_minimum(xr_maximum(Dd + diff_dr, 0), TDW)\n",
+ " \n",
+ " # Write outputs\n",
+ " model_outputs['Irr'].loc[{'time': dates[0]}] = Irrig\n",
+ " \n",
+ " # Update variable_t1 values\n",
+ " for variable in calculation_variables_t1:\n",
+ " variables_t1[variable] = variables_t2[variable].copy(deep = True)\n",
+ " \n",
+ " #============ Time loop ============#\n",
+ " for i in range(1, len(dates)):\n",
+ " \n",
+ " # Update variables\n",
+ " ## Fraction cover\n",
+ " FCov = s_Fslope * Fslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]})/255) + s_Foffset * Foffset_\n",
+ " FCov = xr_minimum(xr_maximum(FCov, 0), s_Fc_stop * Fc_stop_)\n",
+ " \n",
+ " ## Root depth upate\n",
+ " Zr = s_minZr * minZr_ + (FCov / (s_FmaxFC * FmaxFC_)) * s_maxZr * (maxZr_ - minZr_)\n",
+ " \n",
+ " # Water capacities\n",
+ " TAW = (soil_params['FC'] - soil_params['WP']) * Zr\n",
+ " TDW = (soil_params['FC'] - soil_params['WP']) * (s_Zsoil * Zsoil_ - Zr)\n",
+ " \n",
+ " # Update depletions\n",
+ " Dr = update_Dr(TAW, TDW, Zr, TAW0, TDW0, Dr0, Dd0, Zr0)\n",
+ " Dd = update_Dd(TAW, TDW, Zr, TAW0, TDW0, Dd0, Zr0)\n",
+ " \n",
+ " # Update param p\n",
+ " p_ = (xr_minimum(xr_maximum(s_p * p_ + 0.04 * (5 - weather_cube['ET0'].sel({'time': dates[i-1]}) / 1000), 0.1), 0.8) * (1 / s_p)).round(0).astype('i2')\n",
+ " \n",
+ " # Irrigation ==============!!!!!\n",
+ " Irrig = xr_minimum(xr_maximum(Dr - weather_cube['tp'].sel({'time': dates[i]}) / 1000, 0), s_Lame_max * Lame_max_) * Irrig_auto_\n",
+ " Irrig = xr.where(Dr > TAW * s_p * p_, Irrig, 0)\n",
+ " \n",
+ " # Kcb\n",
+ " Kcb = xr_minimum(s_Kslope * Kslope_ * (ndvi_cube['ndvi'].sel({'time': dates[i]}) / 255) + s_Koffset * Koffset_, s_KmaxKcb * KmaxKcb_)\n",
+ " \n",
+ " # Update depletions with rainfall and/or irrigation \n",
+ " ## DP\n",
+ " model_outputs['DP'].loc[{'time': dates[i]}] = -xr_minimum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), 0)\n",
+ " \n",
+ " ## De\n",
+ " Dei = xr_minimum(xr_maximum(Dei - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig / (s_FW * FW_ / 100), 0), TEW)\n",
+ " Dep = xr_minimum(xr_maximum(Dep - (weather_cube['tp'].sel({'time': dates[i]}) / 1000), 0), TEW)\n",
+ " \n",
+ " fewi = xr_minimum(1 - FCov, (s_FW * FW_ / 100))\n",
+ " fewp = 1 - FCov - fewi\n",
+ " \n",
+ " De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+ " # variables_t1['De'] = xr.where(variables_t1['De'].isfinite(), variables_t1['De'], variables_t1['Dei'] * (s_FW * FW_ / 100) + variables_t1['Dep'] * (1 - (s_FW * FW_ / 100))) #================= find replacement for .isfinite() method !!!!!!!!!\n",
+ "\n",
+ " ## Dr\n",
+ " Dr = xr_minimum(xr_maximum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), TAW)\n",
+ " \n",
+ " ## Dd\n",
+ " Dd = xr_minimum(xr_maximum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), 0), TDW)\n",
+ " \n",
+ " # Diffusion coefficients\n",
+ " diff_rei = calculate_diff_re(TAW, Dr, Zr, RUE, Dei, FCov, Ze_, DiffE_, scale_factor)\n",
+ " diff_rep = calculate_diff_re(TAW, Dr, Zr, RUE, Dep, FCov, Ze_, DiffE_, scale_factor)\n",
+ " diff_dr = calculate_diff_dr(TAW, TDW, Dr, Zr, Dd, FCov, Zsoil_, DiffR_, scale_factor) \n",
+ " \n",
+ " # Weighing factor W\n",
+ " W = calculate_W(TEW, Dei, Dep, fewi, fewp)\n",
+ " \n",
+ " # Soil water contents\n",
+ " model_outputs['SWCe'].loc[{'time': dates[i]}] = 1 - De/TEW\n",
+ " model_outputs['SWCr'].loc[{'time': dates[i]}] = 1 - Dr/TAW\n",
+ " \n",
+ " # Water Stress coefficient\n",
+ " Ks = xr_minimum((TAW - Dr) / (TAW * (1 - s_p * p_)), 1)\n",
+ " \n",
+ " # Reduction coefficient for evaporation\n",
+ " Kei = xr_minimum(W * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewi * s_Kcmax * Kcmax_)\n",
+ " Kep = xr_minimum((1 - W) * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewp * s_Kcmax * Kcmax_)\n",
+ " \n",
+ " # Prepare coefficients for evapotranspiration\n",
+ " Kti = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dei / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+ " Ktp = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dep / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+ " Tei = Kti * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+ " Tep = Ktp * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+ " \n",
+ " # Update depletions\n",
+ " Dei = xr.where(fewi > 0, xr_minimum(xr_maximum(Dei + (weather_cube['ET0'].sel({'time': dates[i]}) / 1000) * Kei / fewi + Tei - diff_rei, 0), TEW), xr_minimum(xr_maximum(Dei + Tei - diff_rei, 0), TEW))\n",
+ " Dep = xr.where(fewp > 0, xr_minimum(xr_maximum(Dep + (weather_cube['ET0'].sel({'time': dates[i]}) / 1000) * Kep / fewp + Tep - diff_rep, 0), TEW), xr_minimum(xr_maximum(Dep + Tep - diff_rep, 0), TEW))\n",
+ " \n",
+ " De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+ " # De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100))) #================= find replacement for .isfinite() method !!!!!!!!!\n",
+ " \n",
+ " # Evaporation\n",
+ " model_outputs['E'].loc[{'time': dates[i]}] = xr_maximum((Kei + Kep) * weather_cube['ET0'].sel({'time': dates[i]}) / 1000, 0)\n",
+ " \n",
+ " # Transpiration\n",
+ " model_outputs['Tr'].loc[{'time': dates[i]}] = Kcb * Ks * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+ " \n",
+ " # Potential evapotranspiration and evaporative fraction ??\n",
+ " \n",
+ " # Update depletions (root and deep zones) at the end of the day\n",
+ " Dr = xr_minimum(xr_maximum(Dr + model_outputs['E'].loc[{'time': dates[i]}] + model_outputs['Tr'].loc[{'time': dates[i]}] - diff_dr, 0), TAW)\n",
+ " Dd = xr_minimum(xr_maximum(Dd + diff_dr, 0), TDW)\n",
+ " \n",
+ " # Write outputs\n",
+ " model_outputs['Irr'].loc[{'time': dates[i]}] = Irrig\n",
+ " \n",
+ " # Update variable_t1 values\n",
+ " for variable in calculation_variables_t1:\n",
+ " variables_t1[variable] = variables_t2[variable].copy(deep = True)\n",
+ " \n",
+ " print('day ', i+1, '/', len(dates), ' ', end = '\\r')\n",
+ " \n",
+ " # Scale the model_outputs variable to save in int16 format\n",
+ " model_outputs = model_outputs * 1000\n",
+ " \n",
+ " # Write encoding dict\n",
+ " encoding_dict = {}\n",
+ " for variable in list(model_outputs.keys()):\n",
+ " encod = {}\n",
+ " encod['dtype'] = 'i2'\n",
+ " encoding_dict[variable] = encod\n",
+ " \n",
+ " # Save model outputs to netcdf\n",
+ " model_outputs.to_netcdf(save_path, encoding = encoding_dict)\n",
+ " \n",
+ " return None"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.\n",
+ "Perhaps you already have a cluster running?\n",
+ "Hosting the HTTP server on port 38039 instead\n",
+ " warnings.warn(\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "<div>\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <h3 style=\"margin-bottom: 0px;\">Client</h3>\n",
+ " <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Client-03e36644-264a-11ee-9c09-00155d33b451</p>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ "\n",
+ " <tr>\n",
+ " \n",
+ " <td style=\"text-align: left;\"><strong>Connection method:</strong> Cluster object</td>\n",
+ " <td style=\"text-align: left;\"><strong>Cluster type:</strong> distributed.LocalCluster</td>\n",
+ " \n",
+ " </tr>\n",
+ "\n",
+ " \n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\"></td>\n",
+ " </tr>\n",
+ " \n",
+ "\n",
+ " </table>\n",
+ "\n",
+ " \n",
+ "\n",
+ " \n",
+ " <details>\n",
+ " <summary style=\"margin-bottom: 20px;\"><h3 style=\"display: inline;\">Cluster Info</h3></summary>\n",
+ " <div class=\"jp-RenderedHTMLCommon jp-RenderedHTML jp-mod-trusted jp-OutputArea-output\">\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\">\n",
+ " </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <h3 style=\"margin-bottom: 0px; margin-top: 0px;\">LocalCluster</h3>\n",
+ " <p style=\"color: #9D9D9D; margin-bottom: 0px;\">a5d645bb</p>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Workers:</strong> 4\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads:</strong> 8\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total memory:</strong> 23.47 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " \n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\"><strong>Status:</strong> running</td>\n",
+ " <td style=\"text-align: left;\"><strong>Using processes:</strong> True</td>\n",
+ "</tr>\n",
+ "\n",
+ " \n",
+ " </table>\n",
+ "\n",
+ " <details>\n",
+ " <summary style=\"margin-bottom: 20px;\">\n",
+ " <h3 style=\"display: inline;\">Scheduler Info</h3>\n",
+ " </summary>\n",
+ "\n",
+ " <div style=\"\">\n",
+ " <div>\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #FFF7E5; border: 3px solid #FF6132; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <h3 style=\"margin-bottom: 0px;\">Scheduler</h3>\n",
+ " <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Scheduler-2d9db510-69be-4ff5-86b0-f61af0227954</p>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Comm:</strong> tcp://127.0.0.1:37671\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Workers:</strong> 4\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads:</strong> 8\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Started:</strong> Just now\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total memory:</strong> 23.47 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " </table>\n",
+ " </div>\n",
+ " </div>\n",
+ "\n",
+ " <details style=\"margin-left: 48px;\">\n",
+ " <summary style=\"margin-bottom: 20px;\">\n",
+ " <h3 style=\"display: inline;\">Workers</h3>\n",
+ " </summary>\n",
+ "\n",
+ " \n",
+ " <div style=\"margin-bottom: 20px;\">\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <details>\n",
+ " <summary>\n",
+ " <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 0</h4>\n",
+ " </summary>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Comm: </strong> tcp://127.0.0.1:45285\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads: </strong> 2\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:41621/status\" target=\"_blank\">http://127.0.0.1:41621/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Memory: </strong> 5.87 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Nanny: </strong> tcp://127.0.0.1:36255\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\"></td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td colspan=\"2\" style=\"text-align: left;\">\n",
+ " <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-jebkacuo\n",
+ " </td>\n",
+ " </tr>\n",
+ "\n",
+ " \n",
+ "\n",
+ " \n",
+ "\n",
+ " </table>\n",
+ " </details>\n",
+ " </div>\n",
+ " </div>\n",
+ " \n",
+ " <div style=\"margin-bottom: 20px;\">\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <details>\n",
+ " <summary>\n",
+ " <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 1</h4>\n",
+ " </summary>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Comm: </strong> tcp://127.0.0.1:42337\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads: </strong> 2\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:37155/status\" target=\"_blank\">http://127.0.0.1:37155/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Memory: </strong> 5.87 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Nanny: </strong> tcp://127.0.0.1:36249\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\"></td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td colspan=\"2\" style=\"text-align: left;\">\n",
+ " <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-iprr4c5j\n",
+ " </td>\n",
+ " </tr>\n",
+ "\n",
+ " \n",
+ "\n",
+ " \n",
+ "\n",
+ " </table>\n",
+ " </details>\n",
+ " </div>\n",
+ " </div>\n",
+ " \n",
+ " <div style=\"margin-bottom: 20px;\">\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <details>\n",
+ " <summary>\n",
+ " <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 2</h4>\n",
+ " </summary>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Comm: </strong> tcp://127.0.0.1:38697\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads: </strong> 2\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:34703/status\" target=\"_blank\">http://127.0.0.1:34703/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Memory: </strong> 5.87 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Nanny: </strong> tcp://127.0.0.1:40667\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\"></td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td colspan=\"2\" style=\"text-align: left;\">\n",
+ " <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-92iu95t2\n",
+ " </td>\n",
+ " </tr>\n",
+ "\n",
+ " \n",
+ "\n",
+ " \n",
+ "\n",
+ " </table>\n",
+ " </details>\n",
+ " </div>\n",
+ " </div>\n",
+ " \n",
+ " <div style=\"margin-bottom: 20px;\">\n",
+ " <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+ " <div style=\"margin-left: 48px;\">\n",
+ " <details>\n",
+ " <summary>\n",
+ " <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 3</h4>\n",
+ " </summary>\n",
+ " <table style=\"width: 100%; text-align: left;\">\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Comm: </strong> tcp://127.0.0.1:32795\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Total threads: </strong> 2\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:41767/status\" target=\"_blank\">http://127.0.0.1:41767/status</a>\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Memory: </strong> 5.87 GiB\n",
+ " </td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td style=\"text-align: left;\">\n",
+ " <strong>Nanny: </strong> tcp://127.0.0.1:35569\n",
+ " </td>\n",
+ " <td style=\"text-align: left;\"></td>\n",
+ " </tr>\n",
+ " <tr>\n",
+ " <td colspan=\"2\" style=\"text-align: left;\">\n",
+ " <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-bx7cwqe3\n",
+ " </td>\n",
+ " </tr>\n",
+ "\n",
+ " \n",
+ "\n",
+ " \n",
+ "\n",
+ " </table>\n",
+ " </details>\n",
+ " </div>\n",
+ " </div>\n",
+ " \n",
+ "\n",
+ " </details>\n",
+ "</div>\n",
+ "\n",
+ " </details>\n",
+ " </div>\n",
+ "</div>\n",
+ " </details>\n",
+ " \n",
+ "\n",
+ " </div>\n",
+ "</div>"
+ ],
+ "text/plain": [
+ "<Client: 'tcp://127.0.0.1:37671' processes=4 threads=8, memory=23.47 GiB>"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "client = Client()\n",
+ "client"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "2023-07-19 17:36:42,921 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: ['waiting'] workers: []\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:42,922 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:42,924 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+ "2023-07-19 17:36:43,297 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: [None] workers: []\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:43,298 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:43,300 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+ "2023-07-19 17:36:43,454 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: [None] workers: []\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:43,455 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+ "NoneType: None\n",
+ "2023-07-19 17:36:43,456 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+ "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+ " return func(*(_execute_task(a, cache) for a in args))\n",
+ "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+ " return func(*(_execute_task(a, cache) for a in args))\n",
+ "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+ " return func(*(_execute_task(a, cache) for a in args))\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "day 2 / 366 \r"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+ " return func(*(_execute_task(a, cache) for a in args))\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "day 42 / 366 \r"
+ ]
+ },
+ {
+ "ename": "KeyboardInterrupt",
+ "evalue": "",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[4], line 14\u001b[0m\n\u001b[1;32m 9\u001b[0m save_path \u001b[39m=\u001b[39m data_path \u001b[39m+\u001b[39m os\u001b[39m.\u001b[39msep \u001b[39m+\u001b[39m \u001b[39m'\u001b[39m\u001b[39moutputs.nc\u001b[39m\u001b[39m'\u001b[39m\n\u001b[1;32m 11\u001b[0m chunk_size \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mx\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m'\u001b[39m\u001b[39my\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m}\n\u001b[0;32m---> 14\u001b[0m run_samir(json_config_file, param_file, ndvi_path, weather_path, soil_path, land_cover_path, chunk_size, save_path)\n",
+ "Cell \u001b[0;32mIn[2], line 734\u001b[0m, in \u001b[0;36mrun_samir\u001b[0;34m(json_config_file, csv_param_file, ndvi_cube_path, weather_cube_path, soil_params_path, land_cover_path, chunk_size, save_path)\u001b[0m\n\u001b[1;32m 730\u001b[0m De \u001b[39m=\u001b[39m (Dei \u001b[39m*\u001b[39m fewi \u001b[39m+\u001b[39m Dep \u001b[39m*\u001b[39m fewp) \u001b[39m/\u001b[39m (fewi \u001b[39m+\u001b[39m fewp)\n\u001b[1;32m 731\u001b[0m \u001b[39m# De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100)))\u001b[39;00m\n\u001b[1;32m 732\u001b[0m \n\u001b[1;32m 733\u001b[0m \u001b[39m# Evaporation\u001b[39;00m\n\u001b[0;32m--> 734\u001b[0m model_outputs[\u001b[39m'\u001b[39m\u001b[39mE\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39mloc[{\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}] \u001b[39m=\u001b[39m xr_maximum((Kei \u001b[39m+\u001b[39m Kep) \u001b[39m*\u001b[39m weather_cube[\u001b[39m'\u001b[39m\u001b[39mET0\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39msel({\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}) \u001b[39m/\u001b[39m \u001b[39m1000\u001b[39m, \u001b[39m0\u001b[39m)\n\u001b[1;32m 736\u001b[0m \u001b[39m# Transpiration\u001b[39;00m\n\u001b[1;32m 737\u001b[0m model_outputs[\u001b[39m'\u001b[39m\u001b[39mTr\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39mloc[{\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}] \u001b[39m=\u001b[39m Kcb \u001b[39m*\u001b[39m Ks \u001b[39m*\u001b[39m weather_cube[\u001b[39m'\u001b[39m\u001b[39mET0\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39msel({\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}) \u001b[39m/\u001b[39m \u001b[39m1000\u001b[39m\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/dataarray.py:223\u001b[0m, in \u001b[0;36m_LocIndexer.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m 220\u001b[0m key \u001b[39m=\u001b[39m \u001b[39mdict\u001b[39m(\u001b[39mzip\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array\u001b[39m.\u001b[39mdims, labels))\n\u001b[1;32m 222\u001b[0m dim_indexers \u001b[39m=\u001b[39m map_index_queries(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array, key)\u001b[39m.\u001b[39mdim_indexers\n\u001b[0;32m--> 223\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array[dim_indexers] \u001b[39m=\u001b[39m value\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/dataarray.py:840\u001b[0m, in \u001b[0;36mDataArray.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m 835\u001b[0m \u001b[39m# DataArray key -> Variable key\u001b[39;00m\n\u001b[1;32m 836\u001b[0m key \u001b[39m=\u001b[39m {\n\u001b[1;32m 837\u001b[0m k: v\u001b[39m.\u001b[39mvariable \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(v, DataArray) \u001b[39melse\u001b[39;00m v\n\u001b[1;32m 838\u001b[0m \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_item_key_to_dict(key)\u001b[39m.\u001b[39mitems()\n\u001b[1;32m 839\u001b[0m }\n\u001b[0;32m--> 840\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mvariable[key] \u001b[39m=\u001b[39m value\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/variable.py:977\u001b[0m, in \u001b[0;36mVariable.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m 974\u001b[0m value \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mmoveaxis(value, new_order, \u001b[39mrange\u001b[39m(\u001b[39mlen\u001b[39m(new_order)))\n\u001b[1;32m 976\u001b[0m indexable \u001b[39m=\u001b[39m as_indexable(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_data)\n\u001b[0;32m--> 977\u001b[0m indexable[index_tuple] \u001b[39m=\u001b[39m value\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/indexing.py:1338\u001b[0m, in \u001b[0;36mNumpyIndexingAdapter.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m 1336\u001b[0m array, key \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_indexing_array_and_key(key)\n\u001b[1;32m 1337\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m-> 1338\u001b[0m array[key] \u001b[39m=\u001b[39m value\n\u001b[1;32m 1339\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m:\n\u001b[1;32m 1340\u001b[0m \u001b[39m# More informative exception if read-only view\u001b[39;00m\n\u001b[1;32m 1341\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m array\u001b[39m.\u001b[39mflags\u001b[39m.\u001b[39mwriteable \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m array\u001b[39m.\u001b[39mflags\u001b[39m.\u001b[39mowndata:\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/array/core.py:1699\u001b[0m, in \u001b[0;36mArray.__array__\u001b[0;34m(self, dtype, **kwargs)\u001b[0m\n\u001b[1;32m 1698\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__array__\u001b[39m(\u001b[39mself\u001b[39m, dtype\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[0;32m-> 1699\u001b[0m x \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mcompute()\n\u001b[1;32m 1700\u001b[0m \u001b[39mif\u001b[39;00m dtype \u001b[39mand\u001b[39;00m x\u001b[39m.\u001b[39mdtype \u001b[39m!=\u001b[39m dtype:\n\u001b[1;32m 1701\u001b[0m x \u001b[39m=\u001b[39m x\u001b[39m.\u001b[39mastype(dtype)\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:381\u001b[0m, in \u001b[0;36mDaskMethodsMixin.compute\u001b[0;34m(self, **kwargs)\u001b[0m\n\u001b[1;32m 357\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mcompute\u001b[39m(\u001b[39mself\u001b[39m, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[1;32m 358\u001b[0m \u001b[39m\"\"\"Compute this dask collection\u001b[39;00m\n\u001b[1;32m 359\u001b[0m \n\u001b[1;32m 360\u001b[0m \u001b[39m This turns a lazy Dask collection into its in-memory equivalent.\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 379\u001b[0m \u001b[39m dask.compute\u001b[39;00m\n\u001b[1;32m 380\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 381\u001b[0m (result,) \u001b[39m=\u001b[39m compute(\u001b[39mself\u001b[39;49m, traverse\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 382\u001b[0m \u001b[39mreturn\u001b[39;00m result\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:660\u001b[0m, in \u001b[0;36mcompute\u001b[0;34m(traverse, optimize_graph, scheduler, get, *args, **kwargs)\u001b[0m\n\u001b[1;32m 652\u001b[0m \u001b[39mreturn\u001b[39;00m args\n\u001b[1;32m 654\u001b[0m schedule \u001b[39m=\u001b[39m get_scheduler(\n\u001b[1;32m 655\u001b[0m scheduler\u001b[39m=\u001b[39mscheduler,\n\u001b[1;32m 656\u001b[0m collections\u001b[39m=\u001b[39mcollections,\n\u001b[1;32m 657\u001b[0m get\u001b[39m=\u001b[39mget,\n\u001b[1;32m 658\u001b[0m )\n\u001b[0;32m--> 660\u001b[0m dsk \u001b[39m=\u001b[39m collections_to_dsk(collections, optimize_graph, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 661\u001b[0m keys, postcomputes \u001b[39m=\u001b[39m [], []\n\u001b[1;32m 662\u001b[0m \u001b[39mfor\u001b[39;00m x \u001b[39min\u001b[39;00m collections:\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:433\u001b[0m, in \u001b[0;36mcollections_to_dsk\u001b[0;34m(collections, optimize_graph, optimizations, **kwargs)\u001b[0m\n\u001b[1;32m 431\u001b[0m \u001b[39mfor\u001b[39;00m opt, val \u001b[39min\u001b[39;00m groups\u001b[39m.\u001b[39mitems():\n\u001b[1;32m 432\u001b[0m dsk, keys \u001b[39m=\u001b[39m _extract_graph_and_keys(val)\n\u001b[0;32m--> 433\u001b[0m dsk \u001b[39m=\u001b[39m opt(dsk, keys, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 435\u001b[0m \u001b[39mfor\u001b[39;00m opt_inner \u001b[39min\u001b[39;00m optimizations:\n\u001b[1;32m 436\u001b[0m dsk \u001b[39m=\u001b[39m opt_inner(dsk, keys, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/array/optimization.py:49\u001b[0m, in \u001b[0;36moptimize\u001b[0;34m(dsk, keys, fuse_keys, fast_functions, inline_functions_fast_functions, rename_fused_keys, **kwargs)\u001b[0m\n\u001b[1;32m 46\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39misinstance\u001b[39m(dsk, HighLevelGraph):\n\u001b[1;32m 47\u001b[0m dsk \u001b[39m=\u001b[39m HighLevelGraph\u001b[39m.\u001b[39mfrom_collections(\u001b[39mid\u001b[39m(dsk), dsk, dependencies\u001b[39m=\u001b[39m())\n\u001b[0;32m---> 49\u001b[0m dsk \u001b[39m=\u001b[39m optimize_blockwise(dsk, keys\u001b[39m=\u001b[39;49mkeys)\n\u001b[1;32m 50\u001b[0m dsk \u001b[39m=\u001b[39m fuse_roots(dsk, keys\u001b[39m=\u001b[39mkeys)\n\u001b[1;32m 51\u001b[0m dsk \u001b[39m=\u001b[39m dsk\u001b[39m.\u001b[39mcull(\u001b[39mset\u001b[39m(keys))\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1080\u001b[0m, in \u001b[0;36moptimize_blockwise\u001b[0;34m(graph, keys)\u001b[0m\n\u001b[1;32m 1078\u001b[0m \u001b[39mwhile\u001b[39;00m out\u001b[39m.\u001b[39mdependencies \u001b[39m!=\u001b[39m graph\u001b[39m.\u001b[39mdependencies:\n\u001b[1;32m 1079\u001b[0m graph \u001b[39m=\u001b[39m out\n\u001b[0;32m-> 1080\u001b[0m out \u001b[39m=\u001b[39m _optimize_blockwise(graph, keys\u001b[39m=\u001b[39;49mkeys)\n\u001b[1;32m 1081\u001b[0m \u001b[39mreturn\u001b[39;00m out\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1154\u001b[0m, in \u001b[0;36m_optimize_blockwise\u001b[0;34m(full_graph, keys)\u001b[0m\n\u001b[1;32m 1151\u001b[0m stack\u001b[39m.\u001b[39mappend(d)\n\u001b[1;32m 1153\u001b[0m \u001b[39m# Merge these Blockwise layers into one\u001b[39;00m\n\u001b[0;32m-> 1154\u001b[0m new_layer \u001b[39m=\u001b[39m rewrite_blockwise([layers[l] \u001b[39mfor\u001b[39;49;00m l \u001b[39min\u001b[39;49;00m blockwise_layers])\n\u001b[1;32m 1155\u001b[0m out[layer] \u001b[39m=\u001b[39m new_layer\n\u001b[1;32m 1157\u001b[0m \u001b[39m# Get the new (external) dependencies for this layer.\u001b[39;00m\n\u001b[1;32m 1158\u001b[0m \u001b[39m# This corresponds to the dependencies defined in\u001b[39;00m\n\u001b[1;32m 1159\u001b[0m \u001b[39m# full_graph.dependencies and are not in blockwise_layers\u001b[39;00m\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1341\u001b[0m, in \u001b[0;36mrewrite_blockwise\u001b[0;34m(inputs)\u001b[0m\n\u001b[1;32m 1339\u001b[0m sub \u001b[39m=\u001b[39m {}\n\u001b[1;32m 1340\u001b[0m \u001b[39m# Map from (id(key), inds or None) -> index in indices. Used to deduplicate indices.\u001b[39;00m\n\u001b[0;32m-> 1341\u001b[0m index_map \u001b[39m=\u001b[39m {(\u001b[39mid\u001b[39m(k), inds): n \u001b[39mfor\u001b[39;00m n, (k, inds) \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(indices)}\n\u001b[1;32m 1342\u001b[0m \u001b[39mfor\u001b[39;00m ii, index \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(new_indices):\n\u001b[1;32m 1343\u001b[0m id_key \u001b[39m=\u001b[39m (\u001b[39mid\u001b[39m(index[\u001b[39m0\u001b[39m]), index[\u001b[39m1\u001b[39m])\n",
+ "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1341\u001b[0m, in \u001b[0;36m<dictcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 1339\u001b[0m sub \u001b[39m=\u001b[39m {}\n\u001b[1;32m 1340\u001b[0m \u001b[39m# Map from (id(key), inds or None) -> index in indices. Used to deduplicate indices.\u001b[39;00m\n\u001b[0;32m-> 1341\u001b[0m index_map \u001b[39m=\u001b[39m {(\u001b[39mid\u001b[39;49m(k), inds): n \u001b[39mfor\u001b[39;00m n, (k, inds) \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(indices)}\n\u001b[1;32m 1342\u001b[0m \u001b[39mfor\u001b[39;00m ii, index \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(new_indices):\n\u001b[1;32m 1343\u001b[0m id_key \u001b[39m=\u001b[39m (\u001b[39mid\u001b[39m(index[\u001b[39m0\u001b[39m]), index[\u001b[39m1\u001b[39m])\n",
+ "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+ ]
+ }
+ ],
+ "source": [
+ "data_path = '/mnt/e/DATA/DEV_inputs_test'\n",
+ "\n",
+ "ndvi_path = data_path + os.sep + 'ndvi.nc'\n",
+ "weather_path = data_path + os.sep + 'weather.nc'\n",
+ "land_cover_path = data_path + os.sep + 'land_cover.nc'\n",
+ "json_config_file = '/home/auclairj/GIT/modspa-pixel/config/config_modspa.json'\n",
+ "param_file = '/home/auclairj/GIT/modspa-pixel/parameters/csv_files/params_samir_test.csv'\n",
+ "soil_path = data_path + os.sep + 'soil.nc'\n",
+ "save_path = data_path + os.sep + 'outputs.nc'\n",
+ "\n",
+ "chunk_size = {'x': 5, 'y': 5, 'time': -1}\n",
+ "\n",
+ "\n",
+ "run_samir(json_config_file, param_file, ndvi_path, weather_path, soil_path, land_cover_path, chunk_size, save_path)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "modspa_pixel",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.9"
+ },
+ "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/input/calculate_ndvi.py b/input/calculate_ndvi.py
index b24a581b6f62d3ff9e15c835f736bcbee01bc750..0796bec7c6539bbd4e0410bea0f474971d37a925 100644
--- a/input/calculate_ndvi.py
+++ b/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
diff --git a/input/lib_era5_land_pixel.py b/input/lib_era5_land_pixel.py
index 0a8f7d3560747d1ff4b531992bb5779dc24aa4d2..93b12c3f94de472ecf8ee17a3db57bb7099f9b7a 100644
--- a/input/lib_era5_land_pixel.py
+++ b/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
diff --git a/parameters/csv_files/params_samir_test.csv b/parameters/csv_files/params_samir_test.csv
new file mode 100644
index 0000000000000000000000000000000000000000..3ffa399662ff91acffa1e55a2bf47b19ee6f6245
--- /dev/null
+++ b/parameters/csv_files/params_samir_test.csv
@@ -0,0 +1,24 @@
+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