Link LAI and landuse classes (#273)

* add new methods to derive lai from landuse mapping

* update docs

* update changelog

* updates after review

docs/api.rst

@@ -36,6 +36,7 @@ Setup components
    WflowModel.setup_glaciers
    WflowModel.setup_lulcmaps
    WflowModel.setup_laimaps
+   WflowModel.setup_laimaps_from_lulc_mapping
    WflowModel.setup_ksathorfrac
    WflowModel.setup_rootzoneclim
    WflowModel.setup_soilmaps
@@ -155,6 +156,7 @@ Setup components
    WflowSedimentModel.setup_reservoirs
    WflowSedimentModel.setup_lulcmaps
    WflowSedimentModel.setup_laimaps
+   WflowSedimentModel.setup_laimaps_from_lulc_mapping
    WflowSedimentModel.setup_canopymaps
    WflowSedimentModel.setup_soilmaps
    WflowSedimentModel.setup_riverwidth
@@ -247,6 +249,9 @@ Wflow workflows
    workflows.river_bathymetry
    workflows.pet
    workflows.landuse
+   workflows.lai
+   workflows.create_lulc_lai_mapping_table
+   workflows.lai_from_lulc_mapping
    workflows.ksathorfrac
    workflows.soilgrids
    workflows.soilgrids_sediment

docs/changelog.rst

@@ -22,13 +22,15 @@ Added
 - new utils method **get_grid_from_config** to get the right wflow staticmaps variable based on the TOML configuration (e.g. detects name in netcdf, value, scale and offset). Only applied now to prepare cold states (e.g. not yet in read_grid). PR #252
 - Added support for the "GLCNMO" land-use dataset, with a default parameter mapping table (similar to the existing tables). PR #272
 - Added the `alpha_h1` parameter (based on land use maps). This parameter represents whether root water uptake reduction at soil water pressure head h1 occurs or not. By default, it is set  to 0.0 for all "non-natural" vegetation (crops) and to 1.0 for all "natural vegetation" PR #272
+- New function **setup_laimaps_from_lulc_mapping** to set leaf area index (LAI) climatology maps per month based on landuse mapping. PR #273
 
 Changed
 -------
 - Basins geometry (**basins**) is now based on the actual wflow model resolution basins, instead of based on the supplied DEM `PR #266 <https://github.com/Deltares/hydromt_wflow/pull/266>`
 - **setup_soilmaps**: the user can now supply variable sbm soil layer thicknesses. The Brooks Corey coefficient is then computed as weighted average over the sbm layers. PR #242
 - **setup_outlets**: the IDs of wflow_subcatch are used to define the outlets IDs rather than [1:n]. PR #247
 - wflow forcing data type should always be float32. PR #268
+- **setup_laimaps**: if a landuse map if provided, setup_laimaps can also prepare landuse mapping tables for LAI. PR #273
 
 Fixed
 -----

docs/user_guide/sediment_model_setup.rst

@@ -44,6 +44,8 @@ a specific method see its documentation.
       - This component derives several wflow maps are derived based on landuse- landcover (LULC) data.
     * - :py:func:`~WflowSedimentModel.setup_laimaps`
       - This component sets leaf area index (LAI) climatology maps per month.
+    * - :py:func:`~WflowSedimentModel.setup_laimaps_from_lulc_mapping`
+      - This component sets leaf area index (LAI) climatology maps per month based on landuse mapping.
     * - :py:func:`~WflowSedimentModel.setup_canopymaps`
       - Setup sediments based canopy height maps.
     * - :py:func:`~WflowSedimentModel.setup_soilmaps`
@@ -52,17 +54,17 @@ a specific method see its documentation.
       - This component sets the river width parameter based on a power-lay relationship with a predictor.
     * - :py:func:`~WflowSedimentModel.setup_riverbedsed`
       - Setup sediments based river bed characteristics maps.
-    * - :py:func:`~WflowModel.setup_outlets`
+    * - :py:func:`~WflowSedimentModel.setup_outlets`
       - This method sets the default gauge map based on basin outlets.
-    * - :py:func:`~WflowModel.setup_gauges`
+    * - :py:func:`~WflowSedimentModel.setup_gauges`
       - This method sets the default gauge map based on a gauges_fn data.
-    * - :py:func:`~WflowModel.setup_areamap`
+    * - :py:func:`~WflowSedimentModel.setup_areamap`
       - Setup area map from vector data to save wflow outputs for specific area.
-    * - :py:func:`~WflowModel.setup_config_output_timeseries`
+    * - :py:func:`~WflowSedimentModel.setup_config_output_timeseries`
       - This method add a new variable/column to the netcf/csv output section of the toml based on a selected gauge/area map.
     * - :py:func:`~WflowSedimentModel.setup_constant_pars`
       - Setup constant parameter maps.
-    * - :py:func:`~WflowModel.setup_grid_from_raster`
+    * - :py:func:`~WflowSedimentModel.setup_grid_from_raster`
       -  Setup staticmaps from raster to add parameters from direct data.
 
 

docs/user_guide/wflow_model_setup.rst

@@ -50,6 +50,8 @@ a specific method see its documentation.
       - This component derives several wflow maps are derived based on landuse- landcover (LULC) data.
     * - :py:func:`~WflowModel.setup_laimaps`
       - This component sets leaf area index (LAI) climatology maps per month.
+    * - :py:func:`~WflowModel.setup_laimaps_from_lulc_mapping`
+      - This component sets leaf area index (LAI) climatology maps per month based on landuse mapping.
     * - :py:func:`~WflowModel.setup_soilmaps`
       - This component derives several (layered) soil parameters based on a database with physical soil properties using available point-scale (pedo)transfer functions (PTFs) from literature with upscaling rules to ensure flux matching across scales.
     * - :py:func:`~WflowModel.setup_ksathorfrac`

hydromt_wflow/wflow.py

@@ -953,13 +953,27 @@ def setup_lulcmaps(
             if wflow_param is not None:
                 self.set_config(wflow_param, name)
 
-    def setup_laimaps(self, lai_fn: Union[str, xr.DataArray]):
+    def setup_laimaps(
+        self,
+        lai_fn: Union[str, xr.DataArray],
+        lulc_fn: Optional[Union[str, xr.DataArray]] = None,
+        lulc_sampling_method: str = "any",
+        lulc_zero_classes: List[int] = [],
+        buffer: int = 2,
+    ):
         """
         Set leaf area index (LAI) climatology maps per month [1,2,3,...,12].
 
         The values are resampled to the model resolution using the average value.
         Currently only directly cyclic LAI data is supported.
 
+        If `lulc_fn` is provided, mapping tables from landuse classes to LAI values
+        will be derived from the LULC data. These tables can then be re-used later if
+        you would like to add new LAI maps derived from this mapping table and new
+        landuse scenarios. We advise to use a larger `buffer` to ensure that LAI values
+        can be assigned for all landuse classes and based on a lage enough sample of the
+        LULC data.
+
         Adds model layers:
 
         * **LAI** map: Leaf Area Index climatology [-]
@@ -974,10 +988,62 @@ def setup_laimaps(self, lai_fn: Union[str, xr.DataArray]):
             * Required variables: 'LAI' [-]
 
             * Required dimensions: 'time' = [1,2,3,...,12] (months)
+        lulc_fn : str, xarray.DataArray, optional
+            Name of RasterDataset source for landuse-landcover data.
+            If provided, the LAI values are mapped to landuse classes and will be saved
+            to a csv file.
+        lulc_sampling_method : str, optional
+            Resampling method for the LULC data to the LAI resolution. Two methods are
+            supported:
+
+            * 'any' (default): if any cell of the desired landuse class is present in
+              the resampling window (even just one), it will be used to derive LAI
+              values. This method is less exact but will provide LAI values for all
+              landuse classes for the high resolution landuse map.
+            * 'mode': the most frequent value in the resampling window is
+              used. This method is less precise as for cells with a lot of different
+              landuse classes, the most frequent value might still be only a small
+              fraction of the cell. More landuse classes should however be covered and
+              it can always be used with the landuse map of the wflow model instead of
+              the original high resolution one.
+            * 'q3': only cells with the most frequent value (mode) and that cover 75%
+              (q3) of the resampling window will be used. This method is more exact but
+              for small basins, you may have less or no samples to derive LAI values
+              for some classes.
+        lulc_zero_classes : list, optional
+            List of landuse classes that should have zero for leaf area index values
+            for example waterbodies, open ocean etc. For very high resolution landuse
+            maps, urban surfaces and bare areas can be included here as well.
+            By default empty.
+        buffer : int, optional
+            Buffer around the region to read the data, by default 2.
         """
         # retrieve data for region
         self.logger.info("Preparing LAI maps.")
-        da = self.data_catalog.get_rasterdataset(lai_fn, geom=self.region, buffer=2)
+        da = self.data_catalog.get_rasterdataset(
+            lai_fn, geom=self.region, buffer=buffer
+        )
+        if lulc_fn is not None:
+            self.logger.info("Preparing LULC-LAI mapping table.")
+            da_lulc = self.data_catalog.get_rasterdataset(
+                lulc_fn, geom=self.region, buffer=buffer
+            )
+            # derive mapping
+            df_lai_mapping = workflows.create_lulc_lai_mapping_table(
+                da_lulc=da_lulc,
+                da_lai=da.copy(),
+                sampling_method=lulc_sampling_method,
+                lulc_zero_classes=lulc_zero_classes,
+                logger=self.logger,
+            )
+            # Save to csv
+            if isinstance(lulc_fn, str):
+                df_fn = f"lai_per_lulc_{lulc_fn}.csv"
+            else:
+                df_fn = "lai_per_lulc.csv"
+            df_lai_mapping.to_csv(join(self.root, df_fn))
+
+        # Resample LAI data to wflow model resolution
         da_lai = workflows.lai(
             da=da,
             ds_like=self.grid,
@@ -987,6 +1053,53 @@ def setup_laimaps(self, lai_fn: Union[str, xr.DataArray]):
         rmdict = {da_lai.dims[0]: "time"}
         self.set_grid(da_lai.rename(rmdict), name="LAI")
 
+    def setup_laimaps_from_lulc_mapping(
+        self,
+        lulc_fn: Optional[Union[str, xr.DataArray]],
+        lai_mapping_fn: Union[str, pd.DataFrame],
+    ):
+        """
+        Derive cyclic LAI maps from a LULC data source and a LULC-LAI mapping table.
+
+        Adds model layers:
+
+        * **LAI** map: Leaf Area Index climatology [-]
+            Resampled from source data using average. Assuming that missing values
+            correspond to bare soil, these are set to zero before resampling.
+
+        Parameters
+        ----------
+        lulc_fn : str, xarray.DataArray
+            Name of RasterDataset source for landuse-landcover data.
+        lai_mapping_fn : str, pd.DataFrame
+            Path to a mapping csv file from landuse in source name to
+            LAI values. The csv file should contain rows with landuse classes
+            and LAI values for each month. The columns should be named as the
+            months (1,2,3,...,12).
+            This table can be created using the :py:meth:`setup_laimaps` method.
+        """
+        self.logger.info(
+            "Preparing LAI maps from LULC data using LULC-LAI mapping table."
+        )
+
+        # read landuse map to DataArray
+        da = self.data_catalog.get_rasterdataset(
+            lulc_fn, geom=self.region, buffer=2, variables=["landuse"]
+        )
+        df_lai_mapping = self.data_catalog.get_dataframe(
+            lai_mapping_fn,
+            driver_kwargs={"index_col": 0},  # only used if fn_map is a file path
+        )
+        # process landuse with LULC-LAI mapping table
+        da_lai = workflows.lai_from_lulc_mapping(
+            da=da,
+            ds_like=self.grid,
+            df=df_lai_mapping,
+            logger=self.logger,
+        )
+        # Add to grid
+        self.set_grid(da_lai, name="LAI")
+
     def setup_config_output_timeseries(
         self,
         mapname: str,