Adding snow model

CHANGELOG.md

@@ -35,4 +35,6 @@ Adding `.finalize()` method - clears up the directory. Especially useful for DA.
 ### 1.6.0
   - now compatible with ewatercycle V2.1 `LumpedMakkinkForcing` which generates evaporation from era5/CMIP.
 #### 1.6.1
-  - bug fix occuring when loading makkink data
+  - bug fix occuring when loading makkink data
+### 1.7.0
+  - new version of HBV bmi which adds snow 

pyproject.toml

@@ -10,7 +10,7 @@ name = "ewatercycle-HBV"
 description = "Implementation of HBV for eWaterCycle"
 readme = "README.md"
 license = "Apache-2.0"
-version = "1.6.1"
+version = "1.7.0"
 authors = [
   { name = "David Haasnoot", email = "[email protected]" },
 ]

src/ewatercycle_HBV/__init__.py

@@ -1 +1 @@
-__version__ = "1.6.1"
+__version__ = "1.7.0"

src/ewatercycle_HBV/forcing.py

@@ -13,10 +13,13 @@
 
 
 RENAME_CAMELS = {'total_precipitation_sum':'pr',
-                      'potential_evaporation_sum':'pev',
-                      'streamflow':'Q'}
+                 'potential_evaporation_sum':'pev',
+                 'streamflow':'Q',
+                 'temperature_2m_min':'tasmin',
+                 'temperature_2m_max':'tasmax',
+                 }
 
-REQUIRED_PARAMS = ["pr", "pev"]
+REQUIRED_PARAMS = ["pr", "pev", "tasmean"]
 class HBVForcing(DefaultForcing):
     """Container for HBV forcing data.
 
@@ -71,6 +74,7 @@ class HBVForcing(DefaultForcing):
     # or pr and pev are supplied seperately - can also be the same dataset
     pr: Optional[str] = ".nc"
     pev: Optional[str] = ".nc"
+    tasmean: Optional[str] = ".nc"
     alpha: Optional[float] = 1.26 # varies per catchment, mostly 1.26?
     test_data_bool: bool = False # allows to use self.from_test_txt()
 
@@ -111,6 +115,9 @@ def from_test_txt(self) -> xr.Dataset:
         df_in.index = df_in.apply(lambda x: pd.Timestamp(f'{int(x.year)}-{int(x.month)}-{int(x.day)}'), axis=1)
         df_in = df_in.drop(columns=["year", "month", "day"])
         df_in.index.name = "time"
+        # test data has no snow but let's add in synthetic temperatures to ensure there's no snow:
+        df_in['tasmean'] = 25
+
         # TODO use netcdf-cf conventions
         ds = xr.Dataset(data_vars=df_in,
                         attrs={
@@ -121,6 +128,8 @@ def from_test_txt(self) -> xr.Dataset:
         ds, ds_name = self.crop_ds(ds, "test")
         self.pev = ds_name
         self.pr = ds_name
+        self.tasmean = ds_name
+
         return ds
 
     def from_camels_txt(self) -> xr.Dataset:
@@ -183,7 +192,8 @@ def from_camels_txt(self) -> xr.Dataset:
         # add attributes
         attrs = {"title": "HBV forcing data",
                  "history": "Created by ewatercycle_HBV.forcing.HBVForcing.from_camels_txt()",
-                 "units": "daylight(s), precipitation(mm/day), mean radiation(W/m2), snow water equivalen(mm), temperature max(C), temperature min(C), vapour pressure(Pa)", }
+                 "units": "daylight(s), precipitation(mm/day), mean radiation(W/m2), snow water equivalen(mm), temperature max(C), temperature min(C), temperature mean(c),vapour pressure(Pa)",
+                 }
 
         # add the data lines with catchment characteristics to the description
         attrs.update(data)
@@ -200,9 +210,11 @@ def from_camels_txt(self) -> xr.Dataset:
                              ds.attrs['elevation(m)'],
                              ds.attrs['lat']
                              )
+        ds['tasmean'] = (ds["tasmin"] + ds["tasmax"]) / 2
         ds, ds_name= self.crop_ds(ds, "CAMELS")
         self.pev = ds_name
         self.pr = ds_name
+        self.tasmean = ds_name
         return ds
 
     def from_external_source(self):
@@ -211,24 +223,28 @@ def from_external_source(self):
             self.file_not_found_error()
 
         # often same file
-        if self.pr == self.pev:
+        if self.pr == self.pev == self.tasmean:
             ds = xr.open_dataset(self.directory / self.pr)
 
+            # make compatile with CARAVAN data style:
             if sum([key in ds.data_vars for key in RENAME_CAMELS.keys()]) == len(RENAME_CAMELS):
                 ds = ds.rename(RENAME_CAMELS)
                 ds = ds.rename_dims({'date': 'time'})
                 ds = ds.rename({'date': 'time'})
+                ds['tasmean'] = (ds["tasmin"] + ds["tasmax"]) / 2
 
             ds, ds_name = self.crop_ds(ds, "external")
             self.pev = ds_name
             self.pr = ds_name
+            self.tasmean = ds_name
             return ds
 
         else:
             # but can also seperate
             ds_pr = xr.open_dataset(self.directory / self.pr)
             ds_pev = xr.open_dataset(self.directory / self.pev)
-            combined_data_vars = list(ds_pr.data_vars) + list(ds_pev.data_vars)
+            ds_tasmean = xr.open_dataset(self.directory / self.tasmean)
+            combined_data_vars = list(ds_pr.data_vars) + list(ds_pev.data_vars) + list(ds_tasmean.data_vars)
             if sum([param in combined_data_vars for param in REQUIRED_PARAMS]) != len(REQUIRED_PARAMS):
                 raise UserWarning(f"Supplied NetCDF files must contain {REQUIRED_PARAMS} respectively")
 
@@ -238,7 +254,10 @@ def from_external_source(self):
             ds_pev, ds_name_pev = self.crop_ds(ds_pev, "external")
             self.pev = ds_name_pev
 
-            return ds_pr, ds_pev
+            ds_tasmean, ds_name_tasmean = self.crop_ds(ds_tasmean, "external")
+            self.tasmean = ds_name_tasmean
+
+            return ds_pr, ds_pev, ds_tasmean
 
     def crop_ds(self, ds: xr.Dataset, name: str):
         start = np.datetime64(self.start_time)

src/ewatercycle_HBV/model.py

@@ -43,6 +43,7 @@ class HBVMethods(eWaterCycleModel):
     _config: dict = {
         "precipitation_file": "",
         "potential_evaporation_file": "",
+        "mean_temperature_file": "",
         "parameters": "",
         "initial_storage": "",
                         }
@@ -68,6 +69,9 @@ def _make_cfg_file(self, **kwargs) -> Path:
             self._config["potential_evaporation_file"] = str(
                 self.forcing.directory / self.forcing.pev
             )
+            self._config["mean_temperature_file"] = str(
+                self.forcing.directory / self.forcing.tasmean)
+
 
         elif type(self.forcing).__name__ == 'GenericLumpedForcing':
                 raise UserWarning("Generic Lumped Forcing does not provide potential evaporation, which this model needs")
@@ -95,23 +99,27 @@ def _make_cfg_file(self, **kwargs) -> Path:
                 ds.to_netcdf(temporary_pr_file)
                 ds.close()
 
+            temporary_tasmean_file = self.forcing.directory / self.forcing.filenames['tas'].replace('tas', 'tasmean')
+            if not temporary_tasmean_file.is_file():
+                ds = xr.open_dataset(self.forcing.directory / self.forcing.filenames['tas'])
+                attributes = ds['tas'].attrs
+                ds['tasmean'] = ds['tas']
+                if ds['tasmean'].mean().values > 200: # adjust for kelvin units
+                    ds['tasmean'] -= 273.15
+                ds['tasmean'].attrs = attributes
+                ds.to_netcdf(temporary_tasmean_file)
+                ds.close()
+
             self._config["precipitation_file"] = str(
                 temporary_pr_file
             )
             self._config["potential_evaporation_file"] = str(
                 temporary_pev_file
             )
 
-        ## possibly add later for snow?
-        # self._config["temperature_file"] = str(
-        #     self.forcing.directory / self.forcing.tas
-        # )
-        # self._config["temperature_min_file"] = str(
-        #     self.forcing.directory / self.forcing.tasmin
-        # )
-        # self._config["temperature_max_file"] = str(
-        #     self.forcing.directory / self.forcing.tasmax
-        # )
+            self._config["mean_temperature_file"] = str(
+                temporary_tasmean_file
+            )
 
         for kwarg in kwargs:  # Write any kwargs to the config. - doesn't overwrite config?
             self._config[kwarg] = kwargs[kwarg]
@@ -153,6 +161,8 @@ def parameters(self) -> ItemsView[str, Any]:
 
             Ks (−): Similarly the slow flow is also modelled as 𝑄𝑆=𝐾𝑠∗𝑆𝑆.
 
+            FM (mm/deg/d): Melt Factor: mm of melt per deg per day
+
         """
         pars: dict[str, Any] = dict(zip(HBV_PARAMS, self._config["parameters"].split(',')))
         return pars.items()
@@ -166,6 +176,7 @@ def states(self) -> ItemsView[str, Any]:
             Su (mm): Unsaturated rootzone storage, water stored accessible to plants
             Sf (mm): Fastflow storage, moving Fast through the soil - preferential flow paths, upper level
             Ss (mm): Groundwater storage, moving Slowly through the soil - deeper grounds water.
+            Sp (mm): SnowPack Storage, amount of snow stored
 
         """
         pars: dict[str, Any] = dict(zip(HBV_STATES, self._config["initial_storage"].split(',')))
@@ -207,7 +218,7 @@ def finalize(self) -> None:
             self.unlink()
 
     def unlink(self):
-        for file in ["potential_evaporation_file", "precipitation_file"]:
+        for file in ["potential_evaporation_file", "precipitation_file","mean_temperature_file"]:
             path = self.forcing.directory / self._config[file]
             if path.is_file():  # often both with be the same, e.g. with camels data.
                 path.unlink()
@@ -217,5 +228,5 @@ def unlink(self):
 class HBV(ContainerizedModel, HBVMethods):
     """The HBV eWaterCycle model, with the Container Registry docker image."""
     bmi_image: ContainerImage = ContainerImage(
-        "ghcr.io/daafip/hbv-bmi-grpc4bmi:v1.3.2"
+        "ghcr.io/daafip/hbv-bmi-grpc4bmi:v1.4.0"
     )