modified examples

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]

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()

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()

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)

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
 

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