From 162dfda32c5b4a02c9d3a67ca31af4fa23da576e Mon Sep 17 00:00:00 2001 From: schwemro Date: Fri, 20 Oct 2023 17:48:42 +0200 Subject: [PATCH] modified examples --- .../StressRes/Moehlin/moehlin_setup.py | 6 +- .../eberbaechle/svat_distributed/make_gif.py | 134 ++++++++++++++++++ .../svat_distributed/plot_average.py | 81 +++++++++++ .../eberbaechle/svat_distributed/plot_sum.py | 25 +++- .../eberbaechle/svat_distributed/svat.py | 27 +++- .../ruetlibach/svat_distributed/svat.py | 9 +- .../oneD_distributed_routing_tutorial/oneD.py | 21 ++- .../oneD_distributed_tutorial/oneD.py | 12 +- .../oneD_event.py | 12 +- .../oneD.py | 14 +- .../svat_distributed/svat.py | 24 ++-- roger/variables.py | 10 +- 12 files changed, 313 insertions(+), 62 deletions(-) create mode 100644 examples/catchment_scale/eberbaechle/svat_distributed/make_gif.py create mode 100644 examples/catchment_scale/eberbaechle/svat_distributed/plot_average.py diff --git a/examples/catchment_scale/StressRes/Moehlin/moehlin_setup.py b/examples/catchment_scale/StressRes/Moehlin/moehlin_setup.py index c3c4f6331..4703cbaa3 100644 --- a/examples/catchment_scale/StressRes/Moehlin/moehlin_setup.py +++ b/examples/catchment_scale/StressRes/Moehlin/moehlin_setup.py @@ -216,8 +216,8 @@ def set_parameters_setup(self, state): self._read_var_from_nc("F_n_h_y", input_dir, para_fn)/100 ) # weight factor of air temperature (-) - vs.ta_weight = update( - vs.ta_weight, + vs.ta_offset = update( + vs.ta_offset, at[2:-2, 2:-2], self._read_var_from_nc("F_t", input_dir, para_fn) ) @@ -313,7 +313,7 @@ def set_forcing(self, state): vs.prec_day, at[:, :, :], vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.prec_weight[:, :, npx.newaxis] ) vs.ta_day = update( - vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.ta_weight[:, :, npx.newaxis] + vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] + vs.ta_offset[:, :, npx.newaxis] ) vs.pet_day = update( vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.pet_weight[:, :, npx.newaxis] diff --git a/examples/catchment_scale/eberbaechle/svat_distributed/make_gif.py b/examples/catchment_scale/eberbaechle/svat_distributed/make_gif.py new file mode 100644 index 000000000..dbba7b1aa --- /dev/null +++ b/examples/catchment_scale/eberbaechle/svat_distributed/make_gif.py @@ -0,0 +1,134 @@ +import os +from pathlib import Path +import xarray as xr +import numpy as onp +import matplotlib.pyplot as plt +import yaml +import click +from PIL import Image +import imageio +import copy +import matplotlib as mpl +import seaborn as sns +mpl.use("agg") +import matplotlib.pyplot as plt # noqa: E402 +mpl.rcParams['font.size'] = 8 +mpl.rcParams['axes.titlesize'] = 8 +mpl.rcParams['axes.labelsize'] = 9 +mpl.rcParams['xtick.labelsize'] = 8 +mpl.rcParams['ytick.labelsize'] = 8 +mpl.rcParams['legend.fontsize'] = 8 +mpl.rcParams['legend.title_fontsize'] = 9 +sns.set_style("ticks") +sns.plotting_context("paper", font_scale=1, rc={'font.size': 8.0, + 'axes.labelsize': 9.0, + 'axes.titlesize': 8.0, + 'xtick.labelsize': 8.0, + 'ytick.labelsize': 8.0, + 'legend.fontsize': 8.0, + 'legend.title_fontsize': 9.0}) + +_lab_unit2 = { + "theta": r"$\theta$ [-]", +} + + +@click.command("main") +def main(): + base_path = Path(__file__).parent + # directory of output + base_path_output = base_path / "output" + if not os.path.exists(base_path_output): + os.mkdir(base_path_output) + # directory of figures + base_path_figs = base_path / "figures" + if not os.path.exists(base_path_figs): + os.mkdir(base_path_figs) + + # load configuration file + file_config = base_path / "config.yml" + with open(file_config, "r") as file: + config = yaml.safe_load(file) + + # load hydrologic simulations + states_hm_file = base_path_output / f"{config['identifier']}.nc" + ds_sim = xr.open_dataset(states_hm_file, engine="h5netcdf") + # load catchment + mask_file = base_path / "mask.nc" + ds_mask = xr.open_dataset(mask_file, engine="h5netcdf") + mask = ds_mask["MASK"].values + + # assign date + days_sim = (ds_sim['Time'].values / onp.timedelta64(24 * 60 * 60, "s")) + ds_sim = ds_sim.assign_coords(date=("Time", days_sim)) + + # make GIF for Online-Documentation + frames = [] + for t in range(1, 3*364): + fig, axes = plt.subplots(3, 1, figsize=(4, 6)) + axes[0].imshow(onp.where(mask, ds_sim['prec'].isel(Time=t).values.T, onp.nan), origin="lower", cmap='viridis_r', vmin=0, vmax=20) + axes[0].set_xticks(onp.arange(-12.5, 80, 25)) + axes[0].set_yticks(onp.arange(-12.5, 53, 25)) + axes[0].set_xticklabels(onp.arange(0, 80, 25) * 25) + axes[0].set_yticklabels(onp.arange(0, 53, 25) * 25) + axes[0].set_xlabel('') + axes[0].set_ylabel('Distance in y-direction [m]') + cmap = copy.copy(plt.cm.get_cmap('viridis_r')) + norm = mpl.colors.Normalize(vmin=0, vmax=50) + axl1 = fig.add_axes([0.81, 0.68, 0.02, 0.2]) + cb1 = mpl.colorbar.ColorbarBase(axl1, cmap=cmap, norm=norm, + orientation='vertical', + ticks=[0, 10, 20]) + cb1.ax.set_yticklabels(['0', '10', '20']) + cb1.set_label(r'Precipitation [mm/day]') + + axes[1].imshow(onp.where(mask, ds_sim['theta'].isel(Time=t).values.T, onp.nan), origin="lower", cmap='viridis_r', vmin=0.15, vmax=0.35) + axes[1].set_xticks(onp.arange(-12.5, 80, 25)) + axes[1].set_yticks(onp.arange(-12.5, 53, 25)) + axes[1].set_xticklabels(onp.arange(0, 80, 25) * 25) + axes[1].set_yticklabels(onp.arange(0, 53, 25) * 25) + axes[1].set_xlabel('') + axes[1].set_ylabel('Distance in y-direction [m]') + cmap = copy.copy(plt.cm.get_cmap('viridis_r')) + norm = mpl.colors.Normalize(vmin=0.15, vmax=0.35) + axl1 = fig.add_axes([0.81, 0.4, 0.02, 0.2]) + cb1 = mpl.colorbar.ColorbarBase(axl1, cmap=cmap, norm=norm, + orientation='vertical', + ticks=[0.15, 0.25, 0.35]) + cb1.ax.set_yticklabels(['<0.15', '0.25', '>0.35']) + cb1.set_label(r'Soil water content [-]') + + axes[2].imshow(onp.where(mask, ds_sim['q_ss'].isel(Time=t).values.T, onp.nan), origin="lower", cmap='viridis_r', vmin=0, vmax=5) + axes[2].set_xticks(onp.arange(-12.5, 80, 25)) + axes[2].set_yticks(onp.arange(-12.5, 53, 25)) + axes[2].set_xticklabels(onp.arange(0, 80, 25) * 25) + axes[2].set_yticklabels(onp.arange(0, 53, 25) * 25) + axes[2].set_xlabel('Distance in x-direction [m]') + axes[2].set_ylabel('Distance in y-direction [m]') + cmap = copy.copy(plt.cm.get_cmap('viridis_r')) + norm = mpl.colors.Normalize(vmin=0, vmax=10) + axl1 = fig.add_axes([0.81, 0.12, 0.02, 0.2]) + cb1 = mpl.colorbar.ColorbarBase(axl1, cmap=cmap, norm=norm, + orientation='vertical', + ticks=[0, 2.5, 5]) + cb1.ax.set_yticklabels(['0', '2.5', '>5']) + cb1.set_label(r'Percolation [mm/day]') + fig.suptitle(f't = {t}d', fontsize=9) + fig.subplots_adjust(wspace=0, hspace=0.2, left=0.2, right=0.8, bottom=0.1, top=0.9) + file = base_path_figs / f"t{t}.png" + fig.savefig(file, dpi=250) + plt.close('all') + img = imageio.v2.imread(file) + frames.append(img) + + file = base_path_figs / "prec_theta_perc.gif" + imageio.mimsave(file, + frames, + fps = 7) + + + return + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/examples/catchment_scale/eberbaechle/svat_distributed/plot_average.py b/examples/catchment_scale/eberbaechle/svat_distributed/plot_average.py new file mode 100644 index 000000000..7a4d91124 --- /dev/null +++ b/examples/catchment_scale/eberbaechle/svat_distributed/plot_average.py @@ -0,0 +1,81 @@ +import os +from pathlib import Path +import xarray as xr +import numpy as onp +import matplotlib.pyplot as plt +import yaml +import click +import matplotlib as mpl +import seaborn as sns +mpl.use("agg") +import matplotlib.pyplot as plt # noqa: E402 +mpl.rcParams['font.size'] = 8 +mpl.rcParams['axes.titlesize'] = 8 +mpl.rcParams['axes.labelsize'] = 9 +mpl.rcParams['xtick.labelsize'] = 8 +mpl.rcParams['ytick.labelsize'] = 8 +mpl.rcParams['legend.fontsize'] = 8 +mpl.rcParams['legend.title_fontsize'] = 9 +sns.set_style("ticks") +sns.plotting_context("paper", font_scale=1, rc={'font.size': 8.0, + 'axes.labelsize': 9.0, + 'axes.titlesize': 8.0, + 'xtick.labelsize': 8.0, + 'ytick.labelsize': 8.0, + 'legend.fontsize': 8.0, + 'legend.title_fontsize': 9.0}) + +_lab_unit2 = { + "theta": r"$\theta$ [-]", +} + + +@click.command("main") +def main(): + base_path = Path(__file__).parent + # directory of output + base_path_output = base_path / "output" + if not os.path.exists(base_path_output): + os.mkdir(base_path_output) + # directory of figures + base_path_figs = base_path / "figures" + if not os.path.exists(base_path_figs): + os.mkdir(base_path_figs) + + # load configuration file + file_config = base_path / "config.yml" + with open(file_config, "r") as file: + config = yaml.safe_load(file) + + # load hydrologic simulations + states_hm_file = base_path_output / f"{config['identifier']}.nc" + ds_sim = xr.open_dataset(states_hm_file, engine="h5netcdf") + mask_file = base_path / "mask.nc" + ds_mask = xr.open_dataset(mask_file, engine="h5netcdf") + mask = ds_mask["MASK"].values + + # assign date + days_sim = (ds_sim['Time'].values / onp.timedelta64(24 * 60 * 60, "s")) + ds_sim = ds_sim.assign_coords(date=("Time", days_sim)) + + vars_sim = ["theta"] + for var_sim in vars_sim: + vals1 = onp.where(mask == 1, onp.mean(ds_sim[var_sim].values, axis=-1).T, onp.nan) + vals = onp.where(vals1 < 0, onp.nan, vals1) + fig, ax = plt.subplots(figsize=(6,4)) + grid_extent = (0, 80*25, 0, 53*25) + im = ax.imshow(vals, extent=grid_extent, cmap='viridis', zorder=2) + plt.colorbar(im, ax=ax, shrink=0.7, label=_lab_unit2[var_sim]) + plt.xlabel('Distance in x-direction [m]') + plt.ylabel('Distance in y-direction [m]') + plt.grid(zorder=-1) + plt.tight_layout() + file = base_path_figs / f"{var_sim}_avg.png" + fig.savefig(file, dpi=300) + plt.close(fig) + + return + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/examples/catchment_scale/eberbaechle/svat_distributed/plot_sum.py b/examples/catchment_scale/eberbaechle/svat_distributed/plot_sum.py index 8c5095790..5f9bf6eb3 100644 --- a/examples/catchment_scale/eberbaechle/svat_distributed/plot_sum.py +++ b/examples/catchment_scale/eberbaechle/svat_distributed/plot_sum.py @@ -1,16 +1,33 @@ import os from pathlib import Path import xarray as xr -import pandas as pd import numpy as onp import matplotlib.pyplot as plt import yaml import click -import roger.tools.labels as labs -import roger.tools.evaluation as eval_utils +import matplotlib as mpl +import seaborn as sns +mpl.use("agg") +import matplotlib.pyplot as plt # noqa: E402 +mpl.rcParams['font.size'] = 8 +mpl.rcParams['axes.titlesize'] = 8 +mpl.rcParams['axes.labelsize'] = 9 +mpl.rcParams['xtick.labelsize'] = 8 +mpl.rcParams['ytick.labelsize'] = 8 +mpl.rcParams['legend.fontsize'] = 8 +mpl.rcParams['legend.title_fontsize'] = 9 +sns.set_style("ticks") +sns.plotting_context("paper", font_scale=1, rc={'font.size': 8.0, + 'axes.labelsize': 9.0, + 'axes.titlesize': 8.0, + 'xtick.labelsize': 8.0, + 'ytick.labelsize': 8.0, + 'legend.fontsize': 8.0, + 'legend.title_fontsize': 9.0}) _lab_unit2 = { "q_ss": "PERC [mm]", + "aet": "AET [mm]", "transp": "TRANSP [mm]", "evap_soil": "$EVAP_{soil}$ [mm]", } @@ -44,7 +61,7 @@ def main(): days_sim = (ds_sim['Time'].values / onp.timedelta64(24 * 60 * 60, "s")) ds_sim = ds_sim.assign_coords(date=("Time", days_sim)) - vars_sim = ["transp", "evap_soil", "q_ss"] + vars_sim = ["transp", "evap_soil", "aet", "q_ss"] for var_sim in vars_sim: vals1 = onp.where(mask == 1, onp.sum(ds_sim[var_sim].values, axis=-1).T, onp.nan) / 3 vals = onp.where(vals1 < 0, onp.nan, vals1) diff --git a/examples/catchment_scale/eberbaechle/svat_distributed/svat.py b/examples/catchment_scale/eberbaechle/svat_distributed/svat.py index 8d46b10da..db2ad1406 100644 --- a/examples/catchment_scale/eberbaechle/svat_distributed/svat.py +++ b/examples/catchment_scale/eberbaechle/svat_distributed/svat.py @@ -159,6 +159,24 @@ def set_parameters_setup(self, state): vs.ks = update(vs.ks, at[2:-2, 2:-2], self._read_var_from_nc("ks", self._base_path, "parameters.nc")) # hydraulic conductivity of bedrock/saturated zone (mm/h) vs.kf = update(vs.kf, at[2:-2, 2:-2], 2500) + # weight factor of precipitation (-) + vs.prec_weight = update( + vs.prec_weight, + at[2:-2, 2:-2], + self._read_var_from_nc("prec_weight", self._base_path, "parameters.nc"), + ) + # offset of air temperature (degC) + vs.ta_offset = update( + vs.ta_offset, + at[2:-2, 2:-2], + self._read_var_from_nc("ta_offset", self._base_path, "parameters.nc"), + ) + # weight factor of potential evapotranspiration (-) + vs.pet_weight = update( + vs.pet_weight, + at[2:-2, 2:-2], + self._read_var_from_nc("pet_weight", self._base_path, "parameters.nc"), + ) @roger_routine def set_parameters(self, state): @@ -226,6 +244,9 @@ def set_forcing_setup(self, state): "year", "month", "doy", + "prec_weight", + "pet_weight", + "ta_offset", ], ) def set_forcing(self, state): @@ -248,21 +269,21 @@ def set_forcing(self, state): at[2:-2, 2:-2, :], self._read_var_from_nc("PREC", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, at[2:-2, 2:-2, :], self._read_var_from_nc("TA", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, at[2:-2, 2:-2, :], self._read_var_from_nc("PET", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/examples/catchment_scale/ruetlibach/svat_distributed/svat.py b/examples/catchment_scale/ruetlibach/svat_distributed/svat.py index af0978333..876447434 100644 --- a/examples/catchment_scale/ruetlibach/svat_distributed/svat.py +++ b/examples/catchment_scale/ruetlibach/svat_distributed/svat.py @@ -217,6 +217,9 @@ def set_forcing_setup(self, state): "year", "month", "doy", + "prec_weight", + "pet_weight", + "ta_offset", ], ) def set_forcing(self, state): @@ -239,21 +242,21 @@ def set_forcing(self, state): at[2:-2, 2:-2, :], self._read_var_from_nc("PREC", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, at[2:-2, 2:-2, :], self._read_var_from_nc("TA", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, at[2:-2, 2:-2, :], self._read_var_from_nc("PET", self._input_dir, "forcing.nc")[ :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + ] * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/examples/hillslope_scale/oneD_distributed_routing_tutorial/oneD.py b/examples/hillslope_scale/oneD_distributed_routing_tutorial/oneD.py index da53e27a7..e270fa32f 100644 --- a/examples/hillslope_scale/oneD_distributed_routing_tutorial/oneD.py +++ b/examples/hillslope_scale/oneD_distributed_routing_tutorial/oneD.py @@ -309,8 +309,8 @@ def set_parameters_setup(self, state): 1, ) # weight factor of air temperature (-) - vs.ta_weight = update( - vs.ta_weight, + vs.ta_offset = update( + vs.ta_offset, at[2:-2, 2:-2], 1, ) @@ -414,10 +414,10 @@ def set_parameters_setup(self, state): # ), # ) # # weight factor of air temperature (-) - # vs.ta_weight = update( - # vs.ta_weight, + # vs.ta_offset = update( + # vs.ta_offset, # at[2:-2, 2:-2], - # self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape( + # self._read_var_from_csv("ta_offset", self._base_path, "parameters.csv").reshape( # settings.nx, settings.ny # ), # ) @@ -500,11 +500,6 @@ def set_forcing_setup(self, state): vs.PREC = update(vs.PREC, at[:], self._read_var_from_nc("PREC", self._input_dir, "forcing.nc")[0, 0, :]) vs.TA = update(vs.TA, at[:], self._read_var_from_nc("TA", self._input_dir, "forcing.nc")[0, 0, :]) vs.PET = update(vs.PET, at[:], self._read_var_from_nc("PET", self._input_dir, "forcing.nc")[0, 0, :]) - # vs.PREC = update(vs.PREC, at[:], 0) - # vs.PREC = update(vs.PREC, at[18:24], 1) - # vs.PREC = update(vs.PREC, at[24:30], 20) - # vs.PREC = update(vs.PREC, at[30:512], 0.01) - # vs.PET = update(vs.PET, at[:], 0.05) @roger_routine def set_forcing(self, state): @@ -526,19 +521,19 @@ def set_forcing(self, state): vs.prec_day, at[:, :, :], vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.prec_weight[:, :, npx.newaxis], + * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.ta_weight[:, :, npx.newaxis], + + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.pet_weight[:, :, npx.newaxis], + * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/examples/hillslope_scale/oneD_distributed_tutorial/oneD.py b/examples/hillslope_scale/oneD_distributed_tutorial/oneD.py index 8d4835c99..747906a38 100644 --- a/examples/hillslope_scale/oneD_distributed_tutorial/oneD.py +++ b/examples/hillslope_scale/oneD_distributed_tutorial/oneD.py @@ -215,10 +215,10 @@ def set_parameters_setup(self, state): ), ) # weight factor of air temperature (-) - vs.ta_weight = update( - vs.ta_weight, + vs.ta_offset = update( + vs.ta_offset, at[2:-2, 2:-2], - self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape( + self._read_var_from_csv("ta_offset", self._base_path, "parameters.csv").reshape( settings.nx, settings.ny ), ) @@ -316,19 +316,19 @@ def set_forcing(self, state): vs.prec_day, at[:, :, :], vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.prec_weight[:, :, npx.newaxis], + * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.ta_weight[:, :, npx.newaxis], + + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.pet_weight[:, :, npx.newaxis], + * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/examples/hillslope_scale/oneD_event_distributed_routing_tutorial/oneD_event.py b/examples/hillslope_scale/oneD_event_distributed_routing_tutorial/oneD_event.py index 260dd0be4..a8398176c 100644 --- a/examples/hillslope_scale/oneD_event_distributed_routing_tutorial/oneD_event.py +++ b/examples/hillslope_scale/oneD_event_distributed_routing_tutorial/oneD_event.py @@ -286,8 +286,8 @@ def set_parameters_setup(self, state): 1, ) # weight factor of air temperature (-) - vs.ta_weight = update( - vs.ta_weight, + vs.ta_offset = update( + vs.ta_offset, at[2:-2, 2:-2], 1, ) @@ -385,10 +385,10 @@ def set_parameters_setup(self, state): # ), # ) # # weight factor of air temperature (-) - # vs.ta_weight = update( - # vs.ta_weight, + # vs.ta_offset = update( + # vs.ta_offset, # at[2:-2, 2:-2], - # self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape( + # self._read_var_from_csv("ta_offset", self._base_path, "parameters.csv").reshape( # settings.nx, settings.ny # ), # ) @@ -478,7 +478,7 @@ def set_forcing(self, state): vs.prec = update( vs.prec, at[2:-2, 2:-2, vs.tau], vs.PREC[npx.newaxis, npx.newaxis, vs.itt] * vs.prec_weight[2:-2, 2:-2] ) - vs.ta = update(vs.ta, at[2:-2, 2:-2], vs.TA[npx.newaxis, npx.newaxis, vs.itt]) + vs.ta = update(vs.ta, at[2:-2, 2:-2], vs.TA[npx.newaxis, npx.newaxis, vs.itt] * vs.ta_offset[2:-2, 2:-2]) vs.event_id = update(vs.event_id, at[vs.tau], 1) @roger_routine diff --git a/examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/oneD.py b/examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/oneD.py index 1b1ed21df..98b550614 100644 --- a/examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/oneD.py +++ b/examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/oneD.py @@ -221,11 +221,11 @@ def set_parameters_setup(self, state): settings.nx, settings.ny ), ) - # weight factor of air temperature (-) - vs.ta_weight = update( - vs.ta_weight, + # offset of air temperature (degC) + vs.ta_offset = update( + vs.ta_offset, at[2:-2, 2:-2], - self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape( + self._read_var_from_csv("ta_offset", self._base_path, "parameters.csv").reshape( settings.nx, settings.ny ), ) @@ -340,19 +340,19 @@ def set_forcing(self, state): vs.prec_day, at[:, :, :], vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.prec_weight[:, :, npx.newaxis], + * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.ta_weight[:, :, npx.newaxis], + + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] - * vs.pet_weight[:, :, npx.newaxis], + * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.z_gw = update(vs.z_gw, at[2:-2, 2:-2, vs.tau], vs.Z_GW[npx.newaxis, npx.newaxis, vs.itt_forc]) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/examples/hillslope_scale/svat_distributed_tutorial/svat_distributed/svat.py b/examples/hillslope_scale/svat_distributed_tutorial/svat_distributed/svat.py index af0978333..ca6d7d8a5 100644 --- a/examples/hillslope_scale/svat_distributed_tutorial/svat_distributed/svat.py +++ b/examples/hillslope_scale/svat_distributed_tutorial/svat_distributed/svat.py @@ -217,6 +217,9 @@ def set_forcing_setup(self, state): "year", "month", "doy", + "prec_weight", + "pet_weight", + "ta_offset", ], ) def set_forcing(self, state): @@ -236,24 +239,21 @@ def set_forcing(self, state): ) vs.prec_day = update( vs.prec_day, - at[2:-2, 2:-2, :], - self._read_var_from_nc("PREC", self._input_dir, "forcing.nc")[ - :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + at[:, :, :], + vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] + * vs.prec_weight[2:-2, 2:-2, npx.newaxis], ) vs.ta_day = update( vs.ta_day, - at[2:-2, 2:-2, :], - self._read_var_from_nc("TA", self._input_dir, "forcing.nc")[ - :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + at[:, :, :], + vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] + + vs.ta_offset[2:-2, 2:-2, npx.newaxis], ) vs.pet_day = update( vs.pet_day, - at[2:-2, 2:-2, :], - self._read_var_from_nc("PET", self._input_dir, "forcing.nc")[ - :, :, vs.itt_forc : vs.itt_forc + 6 * 24 - ], + at[:, :, :], + vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] + * vs.pet_weight[2:-2, 2:-2, npx.newaxis], ) vs.itt_forc = vs.itt_forc + 6 * 24 diff --git a/roger/variables.py b/roger/variables.py index 081047431..ba9c20f61 100644 --- a/roger/variables.py +++ b/roger/variables.py @@ -820,12 +820,12 @@ def remove_ghosts(array, dims): "air temperature", time_dependent=True, ), - "ta_weight": Variable( - "weight factor of air temperature", + "ta_offset": Variable( + "offset of air temperature", CATCH_GRID, - "-", - "weight factor of air temperature", - initial=1.0, + "degC", + "offset of air temperature", + initial=0, active=lambda settings: not settings.enable_offline_transport, ), "ta_day": Variable(