diff --git a/.gitignore b/.gitignore
index ee82cd4b4d569b91cc58b7028e8b4cea21ecc4a6..4b211e9bf052029d5b82e50e8a16a01524d92f7f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
tests.py
+test_samir.py
+test_numpy_xarray.py
*__pycache__*
*config_modspa.json
dl_S2.csv
diff --git a/dev_samir_xarray.ipynb b/dev_samir_xarray.ipynb
index ac5bdc27c022ac9ea739ed9f217028f77cdfa8aa..19265d3536a9c8d6a0c0dd4da227296ac9010903 100644
--- a/dev_samir_xarray.ipynb
+++ b/dev_samir_xarray.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -19,7 +19,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -451,7 +451,7 @@
" 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",
+ "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, max_GB: int = 2) -> None:\n",
" \n",
" # warnings.simplefilter(\"error\", category = RuntimeWarning())\n",
" warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in cast\")\n",
@@ -488,6 +488,13 @@
" \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",
+ " # Manage loading of data based on disk size of inputs\n",
+ " if ndvi_cube.nbytes < max_GB * (1024)**3:\n",
+ " ndvi_cube.load()\n",
+ " \n",
+ " if weather_cube.nbytes < max_GB * (1024)**3:\n",
+ " weather_cube.load()\n",
"\n",
" #============ Prepare outputs ============#\n",
" model_outputs = prepare_outputs(ndvi_cube.drop_vars(['ndvi']))\n",
@@ -498,62 +505,62 @@
" #============ 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",
+ " 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",
+ " 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",
+ " 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",
@@ -580,17 +587,17 @@
" \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 = 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",
+ " 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",
@@ -599,19 +606,19 @@
" 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_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",
+ " 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",
+ " 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",
+ " 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",
@@ -620,10 +627,10 @@
" # 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",
+ " 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",
+ " 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",
@@ -634,8 +641,8 @@
" 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",
+ " 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",
@@ -647,30 +654,30 @@
" # 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",
+ " 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",
+ " 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",
+ " 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",
+ " 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",
+ " 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",
+ " model_outputs.Irr.loc[{'time': dates[0]}] = Irrig\n",
" \n",
" # Update variable_t1 values\n",
" for variable in calculation_variables_t1:\n",
@@ -681,37 +688,37 @@
" \n",
" # Update variables\n",
" ## Fraction cover\n",
- " FCov = s_Fslope * Fslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]})/255) + s_Foffset * Foffset_\n",
+ " FCov = s_Fslope * Fslope_ * (ndvi_cube.ndvi.sel({'time': dates[i]})/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",
+ " 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",
+ " 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_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",
+ " 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",
+ " 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",
+ " 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",
@@ -720,10 +727,10 @@
" # 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",
+ " 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",
+ " 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",
@@ -734,8 +741,8 @@
" 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",
+ " 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",
@@ -747,30 +754,30 @@
" # 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",
+ " 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",
+ " 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",
+ " 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",
+ " 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",
+ " 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",
+ " model_outputs.Irr.loc[{'time': dates[i]}] = Irrig\n",
" \n",
" # Update variable_t1 values\n",
" for variable in calculation_variables_t1:\n",
@@ -796,332 +803,9 @@
},
{
"cell_type": "code",
- "execution_count": 3,
+ "execution_count": null,
"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 37667 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-29d78c7f-2653-11ee-9cc7-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:37667/status\" target=\"_blank\">http://127.0.0.1:37667/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;\">c5e1ba98</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:37667/status\" target=\"_blank\">http://127.0.0.1:37667/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-205d3daa-9675-4977-8d69-da12e45dc32c</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:42111\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:37667/status\" target=\"_blank\">http://127.0.0.1:37667/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:43845\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:44805/status\" target=\"_blank\">http://127.0.0.1:44805/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:46421\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-t_o8kxq0\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:34535\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:33207/status\" target=\"_blank\">http://127.0.0.1:33207/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:36817\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-ekfvxctp\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:38783\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:36777/status\" target=\"_blank\">http://127.0.0.1:36777/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:33311\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-pq9c33gu\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:43915\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:37865/status\" target=\"_blank\">http://127.0.0.1:37865/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:36963\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-b4dzhye5\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:42111' processes=4 threads=8, memory=23.47 GiB>"
- ]
- },
- "execution_count": 3,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
+ "outputs": [],
"source": [
"client = Client()\n",
"client"
@@ -1129,84 +813,9 @@
},
{
"cell_type": "code",
- "execution_count": 4,
+ "execution_count": null,
"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: "
- ]
- }
- ],
+ "outputs": [],
"source": [
"data_path = '/mnt/e/DATA/DEV_inputs_test'\n",
"\n",
@@ -1218,11 +827,27 @@
"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",
+ "chunk_size = {'x': -1, 'y': -1, '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)"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "client.close()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/input/calculate_ndvi.py b/input/calculate_ndvi.py
index 0796bec7c6539bbd4e0410bea0f474971d37a925..30a330f7ae6c717a860cb7afda792d23ae954d6b 100644
--- a/input/calculate_ndvi.py
+++ b/input/calculate_ndvi.py
@@ -15,11 +15,16 @@ import xarray as xr # to manage dataset
import pandas as pd # to manage dataframes
import rasterio as rio # to open geotiff files
import geopandas as gpd # to manage shapefile crs projections
+from numpy import nan # to use xr.interpolate_na()
from shapely.geometry import box # to create boundary box
+from config.config import config # to import config file
from input.input_toolbox import product_str_to_datetime
-def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, boundary_shapefile_path: str, resolution: int = 20, chunk_size: dict = {'x': 4000, 'y': 4000, 'time': 8}, acorvi_corr: int = 500) -> str:
+def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, boundary_shapefile_path: str, config_file: str, resolution: int = 20, chunk_size: dict = {'x': 512, 'y': 256, 'time': -1}, acorvi_corr: int = 500) -> str:
+
+ # Open config_file
+ config_params = config(config_file)
# Check resolution for Sentinel-2
if not resolution in [10, 20]:
@@ -72,13 +77,13 @@ def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, bounda
dates = [product_str_to_datetime(prod) for prod in red_paths]
# Open datasets with xarray
- red = xr.open_mfdataset(red_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'red'}).astype('f4')
- nir = xr.open_mfdataset(nir_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'nir'}).astype('f4')
- mask = xr.open_mfdataset(mask_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'mask'}).astype('f4')
+ red = xr.open_mfdataset(red_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'red'})
+ nir = xr.open_mfdataset(nir_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'nir'})
+ mask = xr.open_mfdataset(mask_paths, combine = 'nested', concat_dim = 'time', chunks = chunk_size, parallel = True).squeeze(dim = ['band'], drop = True).rename({'band_data': 'mask'})
if resolution == 10:
- mask = xr.where((mask == 4) | (mask == 5), 1, 0).interp(x = red.coords['x'], y = red.coords['y'], method = 'nearest')
+ mask = xr.where((mask == 4) | (mask == 5), 1, nan).interp(x = red.coords['x'], y = red.coords['y'], method = 'nearest')
else:
- mask = xr.where((mask == 4) | (mask == 5), 1, 0)
+ mask = xr.where((mask == 4) | (mask == 5), 1, nan)
# Set time coordinate
red['time'] = pd.to_datetime(dates)
@@ -94,8 +99,20 @@ def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, bounda
# Mask and scale ndvi
ndvi['ndvi'] = xr.where(ndvi.ndvi < 0, 0, ndvi.ndvi)
ndvi['ndvi'] = xr.where(ndvi.ndvi > 1, 1, ndvi.ndvi)
- ndvi['ndvi'] = (ndvi.ndvi*255).sortby('time')
+ ndvi['ndvi'] = (ndvi.ndvi*255).chunk(chunk_size)
+
+ # Sort images by time
+ ndvi = ndvi.sortby('time')
+ # Interpolates on a daily frequency
+ daily_index = pd.date_range(start = config_params.start_date, end = config_params.end_date, freq = 'D')
+
+ # Resample the dataset to a daily frequency and reindex with the new DateTimeIndex
+ ndvi = ndvi.resample(time = '1D').asfreq().reindex(time = daily_index)
+
+ # Interpolate the dataset along the time dimension to fill nan values
+ ndvi = ndvi.interpolate_na(dim = 'time', method = 'linear', fill_value = 'extrapolate').astype('u1')
+
# Write attributes
ndvi['ndvi'].attrs['units'] = 'None'
ndvi['ndvi'].attrs['standard_name'] = 'NDVI'
@@ -103,7 +120,7 @@ def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, bounda
ndvi['ndvi'].attrs['scale factor'] = '255'
# Create save path
- ndvi_cube_path = save_dir + os.sep + 'NDVI_precube_' + dates[0].strftime('%d-%m-%Y') + '_' + dates[-1].strftime('%d-%m-%Y') + '.nc'
+ ndvi_cube_path = save_dir + os.sep + 'NDVI_cube_' + dates[0].strftime('%d-%m-%Y') + '_' + dates[-1].strftime('%d-%m-%Y') + '.nc'
# Save NDVI cube to netcdf
ndvi.to_netcdf(ndvi_cube_path, encoding = {"ndvi": {"dtype": "u1", "_FillValue": 0}})
diff --git a/input/download_ERA5_weather.py b/input/download_ERA5_weather.py
index cc63de66ab0c12d91dbf0f4e6355007684c5aa0a..b82a9cc753195391ac168d968e5365b4ff31681b 100644
--- a/input/download_ERA5_weather.py
+++ b/input/download_ERA5_weather.py
@@ -18,7 +18,7 @@ import input.lib_era5_land_pixel as era5land # custom built functions for ERA5-
from config.config import config # to import config file
-def request_ER5_weather(input_file: str, ndvi_path: str) -> str:
+def request_ER5_weather(input_file: str, raw_S2_image_ref: str) -> str:
# Get config file
config_params = config(input_file)
@@ -120,7 +120,7 @@ def request_ER5_weather(input_file: str, ndvi_path: str) -> str:
print('----------')
# Save daily wheather data into ncfile
- weather_daily_ncFile = save_dir + os.sep + config_params.start_date + '_' + config_params.end_date + '_' + config_params.run_name + '_era5-land-daily-meteo.nc'
+ weather_daily_ncFile = save_dir + os.sep + config_params.start_date + '_' + config_params.end_date + '_' + config_params.run_name + '_era5-land-daily-meteo'
# Temporary save directory for daily file merge
variable_list = ['2m_dewpoint_temperature_daily_maximum', '2m_dewpoint_temperature_daily_minimum', '2m_temperature_daily_maximum', '2m_temperature_daily_minimum', 'total_precipitation_daily_mean', '10m_u_component_of_wind_daily_mean', '10m_v_component_of_wind_daily_mean', 'surface_solar_radiation_downwards_daily_mean']
@@ -129,7 +129,8 @@ def request_ER5_weather(input_file: str, ndvi_path: str) -> str:
aggregated_files = era5land.concat_monthly_nc_file(list_era5land_hourly_ncFiles, variable_list, save_dir)
# Calculate ET0 over the whole time period
- era5land.era5Land_nc_daily_to_ET0(aggregated_files, weather_daily_ncFile, ndvi_path, h = wind_height)
- print(weather_daily_ncFile)
+ era5land.era5Land_nc_daily_to_ET0(aggregated_files, weather_daily_ncFile, raw_S2_image_ref, config_params, h = wind_height)
+
+ print('\n', weather_daily_ncFile + '.nc', '\n')
- return weather_daily_ncFile
\ No newline at end of file
+ return weather_daily_ncFile + '.nc'
\ No newline at end of file
diff --git a/input/lib_era5_land_pixel.py b/input/lib_era5_land_pixel.py
index 93b12c3f94de472ecf8ee17a3db57bb7099f9b7a..63fca42d41e51b7ba1e52f9c1711b325313e10b5 100644
--- a/input/lib_era5_land_pixel.py
+++ b/input/lib_era5_land_pixel.py
@@ -10,7 +10,7 @@ Functions to call ECMWF Reanalysis with CDS-api
@author: rivalland
"""
-import os # for path exploration
+import os, shutil # for path exploration and file management
from typing import List # to declare variables
import numpy as np # for math on arrays
import xarray as xr # to manage nc files
@@ -18,6 +18,11 @@ from datetime import datetime # to manage dates
from p_tqdm import p_map # for multiprocessing with progress bars
from dateutil.rrule import rrule, MONTHLY
from fnmatch import fnmatch # for file name matching
+import rasterio # to manage geotiff images
+from pandas import date_range
+from rasterio.warp import reproject, Resampling # to reproject
+from dask.diagnostics import ProgressBar
+
import re # for string comparison
import warnings # to suppress pandas warning
@@ -429,19 +434,64 @@ 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, h: float = 10) -> None:
+def reproject_geotiff(source_image: str, destination_image: str, destination_crs: str):
+
+ # Open the original GeoTIFF file
+ with rasterio.open(source_image) as src:
+ # Get the source CRS and transform
+ src_crs = src.crs
+ src_transform = src.transform
+ # Read the data as a numpy array
+ source = src.read()
+
+ # Optionally, calculate the destination transform and shape based on the new CRS
+ dst_transform, dst_width, dst_height = rasterio.warp.calculate_default_transform(
+ src_crs, destination_crs, src.width, src.height, *src.bounds)
+
+ # Create an empty numpy array for the destination
+ destination = np.zeros((src.count, dst_height, dst_width))
+
+ # Reproject the source to the destination
+ reproject(
+ source,
+ destination,
+ src_transform=src_transform,
+ src_crs=src_crs,
+ dst_transform=dst_transform,
+ dst_crs=destination_crs,
+ resampling=Resampling.nearest)
+
+ # Save the reprojected data as a new GeoTIFF file
+ with rasterio.open(destination_image, "w", **src.meta) as dst:
+ # Update the metadata with the new CRS, transform and shape
+ dst.update(
+ crs=destination_crs,
+ transform=dst_transform,
+ width=dst_width,
+ height=dst_height)
+ # Write the reprojected data to the file
+ dst.write(destination)
+
+ return None
+
+
+def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_file: str, raw_S2_image_ref: str, config_params, h: float = 10, max_ram: int = 12288) -> None:
"""
Calculate ET0 values from the ERA5 netcdf weather variables.
Output netcdf contains the ET0 values for each day in the selected
- time period and for each ERA5 pixel covering the required area.
+ time period and reprojected on the same grid as the NDVI values.
## Arguments
1. list_era5land_files: `List[str]`
list of netcdf files containing the necessary variables
- 2. output_nc_file: `str`
- output netcdf file to save
- 3. h: `float` `default = 10`
+ 2. output_file: `str`
+ output file name without extension
+ 3. raw_S2_image_ref: `str`
+ raw Sentinel 2 image at right resolution for reprojection
+ 4. h: `float` `default = 10`
height of ERA5 wind measurements in meters
+ 5. max_ram: `int` `default = 12288`
+ max ram (in MiB) to give to OTB
## Returns
`None`
@@ -477,7 +527,10 @@ def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str
final_weather_ds['tp'] = final_weather_ds['tp'] * 1000 # conversion from m to mm
# Change datatype to reduce memory usage
- final_weather_ds = (final_weather_ds * 1000).astype('u2')
+ final_weather_ds = (final_weather_ds * 1000).astype('u2')
+
+ # Write projection
+ final_weather_ds = final_weather_ds.rio.write_crs('EPSG:4326')
# Set variable attributes
final_weather_ds['ET0'].attrs['units'] = 'mm'
@@ -487,9 +540,28 @@ def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str
final_weather_ds['tp'].attrs['units'] = 'mm'
final_weather_ds['tp'].attrs['standard_name'] = 'Precipitation'
final_weather_ds['tp'].attrs['comment'] = 'Volume of total daily precipitation expressed as water height in milimeters'
- final_weather_ds['tp'].attrs['scale factor'] = '1000'
+ final_weather_ds['tp'].attrs['scale factor'] = '1000'
- # 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"}})
+ # Save dataset to geotiff, still in wgs84 (lat, lon) coordinates
+ output_file_prec = output_file + '_prec.tif'
+ output_file_ET0 = output_file + '_ET0.tif'
+ final_weather_ds.tp.rio.to_raster(output_file_prec, dtype = 'uint16')
+ final_weather_ds.ET0.rio.to_raster(output_file_ET0, dtype = 'uint16')
+
+ # Reprojected image paths
+ output_file_prec_reproj = output_file + '_prec_reproj.tif'
+ output_file_ET0_reproj = output_file + '_ET0_reproj.tif'
+ # Run otbcli_SuperImpose
+ OTB_command = 'otbcli_Superimpose -inr ' + raw_S2_image_ref + ' -inm ' + output_file_prec + ' -out ' + output_file_prec_reproj + ' uint16 -ram ' + str(max_ram)
+ os.system(OTB_command)
+ OTB_command = 'otbcli_Superimpose -inr ' + raw_S2_image_ref + ' -inm ' + output_file_ET0 + ' -out ' + output_file_ET0_reproj + ' uint16 -ram ' + str(max_ram)
+ os.system(OTB_command)
+
+ # remove old files and rename outputs
+ os.remove(output_file_prec)
+ shutil.move(output_file_prec_reproj, output_file_prec)
+ os.remove(output_file_ET0)
+ shutil.move(output_file_ET0_reproj, output_file_ET0)
+
return None
\ No newline at end of file