Water demand

docs/api.rst

@@ -36,6 +36,8 @@ Setup components
    WflowModel.setup_glaciers
    WflowModel.setup_lulcmaps
    WflowModel.setup_laimaps
+   WflowModel.setup_allocation
+   WflowModel.setup_non_irigation
    WflowModel.setup_rootzoneclim
    WflowModel.setup_soilmaps
    WflowModel.setup_outlets
@@ -238,6 +240,8 @@ Wflow workflows
 .. autosummary::
    :toctree: _generated
 
+   workflows.allocate
+   workflows.non_irigation
    workflows.hydrography
    workflows.topography
    workflows.river

docs/changelog.rst

@@ -12,6 +12,10 @@ Unreleased
 Added
 -----
 - **setup_rivers**: Add river depth based on rivdph columns in river_geom_fn rather than only computed from qbankfull column.
+- Added setup for water demand for wflow
+- Added setup method **setup_allocation** to **WflowModel**
+- Added setup method **setup_non_irigation** to **WflowModel**
+- Added workflow functions **allocation** & **non_irigation** for water demand setup
 
 Changed
 -------

hydromt_wflow/wflow.py

@@ -5,6 +5,7 @@
 import glob
 import logging
 import os
+from itertools import product
 from os.path import basename, dirname, isdir, isfile, join
 from pathlib import Path
 from typing import List, Optional, Union
@@ -65,6 +66,9 @@ class WflowModel(GridModel):
         "glacareas": "wflow_glacierareas",
         "glacfracs": "wflow_glacierfrac",
         "glacstore": "wflow_glacierstore",
+        "paddy_area": "paddy_irrigation_areas",
+        "nonpaddy_area": "nonpaddy_irrigation_areas",
+        "crop_factor": "crop_factor",
     }
     _FOLDERS = [
         "staticgeoms",
@@ -2494,6 +2498,292 @@ def setup_rootzoneclim(
         # update config
         self.set_config("input.vertical.rootingdepth", update_toml_rootingdepth)
 
+    def setup_allocation(
+        self,
+        min_area: float | int = 0,
+        admin_bounds_fn: str = "gadm",
+        admin_level: int = None,
+    ):
+        """_summary_.
+
+        _extended_summary_
+
+        Parameters
+        ----------
+        min_area : float | int
+            _description_
+        admin_bounds_fn : str, optional
+            _description_, by default "gadm"
+        admin_level : int, optional
+            _description_, by default 0
+        """
+        self.logger.info("Preparing water demand allocation map.")
+        # Will be fixes but for know this is done like this
+        # TODO fix in the future
+        admin_bounds = None
+        if admin_bounds_fn is not None:
+            if admin_level is not None:
+                admin_bounds_fn += f"_level{admin_level}"
+            admin_bounds = self.data_catalog.get_geodataframe(
+                admin_bounds_fn,
+                geom=self.region,
+            )
+            # Add this identifier for usage in the workflow
+            admin_bounds["admin_id"] = range(len(admin_bounds))
+
+        # Create the allocation grid
+        alloc = workflows.demand.allocate(
+            ds_like=self.grid,
+            min_area=min_area,
+            admin_bounds=admin_bounds,
+            basins=self.geoms["basins"],
+            rivers=self.geoms["rivers"],
+        )
+        self.set_grid(alloc)
+
+        # Update the settings toml
+        self.set_config("input.vertical.waterallocation.areas", "Allocation_id")
+
+    def setup_non_irigation(
+        self,
+        non_irigation_fn: str = "pcr_globwb",
+        non_irigation_vars: list = ["dom", "ind", "lsk"],
+        non_irigation_year: int = None,
+        non_irigation_method: str = "nearest",
+        population_fn: str = "worldpop_2020_constrained",
+        population_method: str = "sum",
+    ):
+        """_summary_.
+
+        _extended_summary_
+
+        Parameters
+        ----------
+        non_irigation_fn : str, optional
+            _description_, by default "pcr_globwb"
+        non_irigation_vars : list, optional
+            _description_, by default ["dom", "ind", "lsk"]
+        non_irigation_year : int, optional
+            _description_, by default None
+        non_irigation_method : str, optional
+            _description_, by default "nearest"
+        population_fn : str, optional
+            _description_, by default "worldpop_2020_constrained"
+        population_method : str, optional
+            _description_, by default "sum"
+
+        Raises
+        ------
+        ValueError
+            _description_
+        """
+        self.logger.info("Preparing non irigation demand maps.")
+        if not all([item in ["dom", "ind", "lsk"] for item in non_irigation_vars]):
+            raise ValueError("")
+
+        # Set flag for cyclic data
+        _cyclic = False
+
+        # Selecting data
+        if non_irigation_year is None:
+            non_irigation_year = 2005
+        non_irigation_raw = self.data_catalog.get_rasterdataset(
+            non_irigation_fn,
+            geom=self.region,
+            buffer=2,
+            variables=[
+                f"{var}_{mode}"
+                for var, mode in product(non_irigation_vars, ["gross", "net"])
+            ],
+            version=non_irigation_year,
+        )
+        if "time" in non_irigation_raw.coords:
+            _cyclic = True
+            non_irigation_raw["time"] = non_irigation_raw.time.astype("int32")
+
+        if _cyclic and self.get_config("input.cyclic") is None:
+            self.set_config("input.cyclic", [])
+
+        # Get population data
+        pop_raw = self.data_catalog.get_rasterdataset(
+            population_fn,
+            geom=self.region,
+            buffer=2,
+        ).raster.mask_nodata()
+
+        # Create static water demand rasters
+        non_irigation, popu = workflows.demand.non_irigation(
+            non_irigation_raw,
+            ds_like=self.grid,
+            ds_method=non_irigation_method,
+            popu=pop_raw,
+            popu_method=population_method,
+        )
+        self.set_grid(non_irigation)
+        self.set_grid(popu)
+
+        # Update the settings toml with non irigation stuff
+        for var in non_irigation.data_vars:
+            sname, suffix = var.split("_")
+            lname = workflows.demand.map_vars[sname]
+            self.set_config(
+                f"input.vertical.{lname}.demand_{suffix}",
+                var,
+            )
+
+            # Also for the fact that these parameters are cyclic
+            if _cyclic:
+                self.config["input"]["cyclic"].append(
+                    f"input.vertical.{lname}.demand_{suffix}",
+                )
+
+    def setup_irrigation(
+        self,
+        irrigated_area_fn: str,
+        landuse_fn: str,
+        paddy_class: int = 12,
+        area_threshold: float = 0.6,
+        crop_irrigated_fn: str = "mirca_irrigated_data",
+        crop_rainfed_fn: str = "mirca_rainfed_data",
+        crop_info_fn: str = "mirca_crop_info"
+        # crop_info_fn: str, TODO, required when adding the support for rootingdepth and
+        # cropfactor
+    ):
+        """
+        Add required information to simulate irrigation water demand.
+
+        The function requires data that contains information about the location of the
+        irrigated areas (`irrigated_area_fn`). This, combined with a landuse data that
+        contains a class for paddy (rice) land use (`landuse_fn`), determines which
+        locations are considered to be paddy irrigation (based on the `paddy_class`),
+        and which locations are considered to be non-paddy irrigation.
+
+        Next, these maps are reprojected to the model resolution, where a threshold
+        (`area_threshold`) determines when pixels are considered to be classified as
+        irrigation cells (both paddy and non-paddy). It adds the resulting maps to the
+        input data.
+
+
+        Adds model layers:
+
+        * **paddy_irrigation_areas**: Irrigated (paddy) mask [-]
+        * **nonpaddy_irrigation_areas**: Irrigated (non-paddy) mask [-]
+
+        Parameters
+        ----------
+        irrigated_area_fn: str
+            Name of the (gridded) dataset that contains the location of irrigated areas
+            (as a mask), `irrigated_area` for example
+        landuse_fn: str
+            Name of the landuse dataset that contains a classification for paddy/rice,
+            use `glcnmo` for example
+        paddy_class: int
+            Class in the landuse data that is considered as paddy or rice, by default 12
+            (matching the glcmno landuse data)
+        area_threshold: float
+            Fractional area of a (wflow) pixel before it gets classified as an irrigated
+            pixel, by default 0.6
+        crop_irrigated_fn: str
+            Name of dataset that contains information about irrigated crops, by default
+            "mirca_irrigated_data",
+        crop_rainfed_fn: str
+            Name of dataset that contains information about rainfed crops, by default
+            "mirca_rainfed_data",
+        crop_info_fn: str
+            Name of dataframe that contains rootingdepth and cropfactor values for the
+            rice/paddy crop class"mirca_crop_info"
+
+
+        See Also
+        --------
+        workflows.demand.find_paddy
+        workflows.demand.classify_pixels
+        """
+        self.logger.info("Preparing irrigation maps.")
+        # Extract irrigated area dataset
+        irrigated_area = self.data_catalog.get_rasterdataset(
+            irrigated_area_fn, bbox=self.grid.raster.bounds, buffer=3
+        )
+
+        # Extract landcover dataset (that includes a paddy/rice class)
+        landuse_da = self.data_catalog.get_rasterdataset(
+            landuse_fn, bbox=self.grid.raster.bounds, buffer=3
+        )
+
+        # Get paddy and nonpaddy masks
+        paddy, nonpaddy = workflows.demand.find_paddy(
+            landuse_da=landuse_da,
+            irrigated_area=irrigated_area,
+            paddy_class=paddy_class,
+        )
+
+        # Get paddy and non paddy pixels at model resolution
+        wflow_paddy = workflows.demand.classify_pixels(
+            da_crop=paddy,
+            da_model=self.grid[self._MAPS["basins"]].raster.mask_nodata(),
+            threshold=area_threshold,
+        )
+        wflow_nonpaddy = workflows.demand.classify_pixels(
+            da_crop=nonpaddy,
+            da_model=self.grid[self._MAPS["basins"]].raster.mask_nodata(),
+            threshold=area_threshold,
+        )
+
+        # Add maps to grid
+        self.set_grid(wflow_paddy, name=self._MAPS["paddy_area"])
+        self.set_grid(wflow_nonpaddy, name=self._MAPS["nonpaddy_area"])
+
+        # Update config
+        self.set_config(
+            "input.vertical.paddy.irrigation_areas", self._MAPS["paddy_area"]
+        )
+        self.set_config(
+            "input.vertical.nonpaddy.irrigation_areas", self._MAPS["nonpaddy_area"]
+        )
+
+        # TODO: Include this support for adjusted crop_factor and rooting depth maps, or
+        # move to seperate function?
+        self.logger.info("Preparing crop factor maps.")
+        crop_rainfed_ds = self.data_catalog.get_rasterdataset(
+            crop_rainfed_fn,
+            bbox=self.grid.raster.bounds,
+            buffer=3,
+        )
+        crop_irrigated_ds = self.data_catalog.get_rasterdataset(
+            crop_irrigated_fn,
+            bbox=self.grid.raster.bounds,
+            buffer=3,
+        )
+
+        df = self.data_catalog.get_dataframe(crop_info_fn)
+        # TODO: Make more flexible
+        rice_value = df.loc["Rice", "kc_mid"]
+
+        cropfactor = workflows.demand.add_crop_maps(
+            ds_rain=crop_rainfed_ds,
+            ds_irri=crop_irrigated_ds,
+            paddy_value=rice_value,
+            mod=self,
+            default_value=1.0,
+            map_type="crop_factor",
+        )
+
+        # Add maps to grid and update config
+        self.set_grid(cropfactor, name=self._MAPS["crop_factor"])
+        self.set_config("input.vertical.kc", self._MAPS["crop_factor"])
+
+        # # TODO: Make more flexible?
+        # rice_value = df.loc["Rice", "rootingdepth_irrigated"]
+
+        # rootingdepth = workflows.demand.add_crop_maps(
+        #     ds_rain=mirca_rain_ds,
+        #     ds_irri=mirca_irri_ds,
+        #     paddy_value=rice_value,
+        #     mod=self,
+        #     default_value=self.grid["RootingDepth"],
+        #     map_type="rootingdepth",
+        # )
+
     # I/O
     def read(self):
         """Read the complete model schematization and configuration from file."""

hydromt_wflow/workflows/__init__.py

@@ -1,6 +1,7 @@
 """HydroMT wflow workflows."""
 
 from .basemaps import *
+from .demand import *
 from .gauges import *
 from .glaciers import *
 from .landuse import *