Include KsatHorFrac estimation method from Awad Ali's MSc thesis

[RubenImhoff]

Kind of request

Adding new functionality

Enhancement Description

In his MSc thesis, Awad Ali tested Random Forest (RF) and Boosted Regression Trees (BRT) methods to estimate the KsatHorFrac parameter of wflow_sbm. The methods show good results and may be a first step to an alternative for the now uniform KsatHorFrac parameter (with value 100). Awad has already prepared two global maps (one based on RF and one based on the BRT algorithm), which can be found on the wflow P-drive (hydromt_staging/kstathorfrac_thesis_awadali/). Because he trained the RF and BRT models on the CAMELS-UK dataset in the UK, the method may not yet be working well in regions with different geology and soil types. But this is a point of further research at the moment.

I think the results are ready for usage in hydromt as an alternative option to the standard uniform KsatHorFrac derivation. In the additional context below, I've provide a function that we use to include one of the KsatHorFrac maps to update an existing wflow_sbm model. Let's see if we can include this in hydromt_wflow in such a way that we can do this up front through a function. Can you assist me in creating a hydromt function that can make use of other KsatHorFrac maps? I am not familiar enough with all functions and procedures in the current version of hydromt to directly open a PR myself, so let's first see where the function belongs and what it should look like. If we make it general enough, you can switch between the 'standard' derivation, the RF and BRT methods (or possibly something else in the future). I can imagine that we update the RF and BRT maps at some point based on future research, but that shouldn't change the hydromt functionality, I assume.

Use case

The KsatHorFrac parameter is currently derived as a uniform parameter with value 100, which has to be calibrated later on for improved model performance. Some sort of transfer function for this parameter would be ideal. With Awad's work, we have maps of gridded KsatHorFrac values that can be used in wflow_sbm. These values are based on a ML-based transfer functions between soilgrids data and the KsatHorFrac value. With the current global maps, we expect good results in regions with geology and soil types similar to the UK. For other regions, this has still to be tested.

Additional Context

Below the function we use to update the wflow model with the new KSatHorFrac maps. The maps are stored at the soilgrids spatial resolution and are store as log10(KSatHorFrac).

import os
import numpy as np

import hydromt
from hydromt_wflow import WflowModel

# A function adjust_ksathorfrac is defined that takes in a model_folder path. 
# The function adjusts the KsatHorFrac parameter grid based on either a 
# random forest or boosted regression trees based global parameter map.
def adjust_ksathorfrac_and_toml_file(
    model_folder, 
    ksathorfrac_path=None, 
    pred_method=None,
    ):
    """
    Adjust the KsatHorFrac parameter grid with empirical relationship.

    Args:
        model_folder (str): The path to the wflow model.
        ksathorfrac_path (str): optional. the path to the .tif file with the global KsatHorFrac 
                                maps based on the random forest or boosted regression trees
                                methods. Standard set to None.
        pred_method (str): optional. 'RF' for Random Forest, and 'BRT' for Boosted Regression 
                                Trees. If None, the KsatHorFrac value will not be changed. 
                                Standard set to None.

    Returns:
        None.    
    """
    # Open the wflow model
    print(f"Reading wflow model at {model_folder}")
    mod = WflowModel(model_folder, config_fn="wflow_sbm.toml", mode="r+")

    # Change KsatHorFrac if requested
    if pred_method is not None:
        if pred_method == "RF" or pred_method == "BRT":
            var_name = pred_method + '_250_KsatHorFrac_global'

            # read (clipped by extent) KsatHorFrac with data_catalog.get_rasterdataset
            da_KsatHorFrac = mod.data_catalog.get_rasterdataset(
                data_like = ksathorfrac_path,
                geom = mod.staticgeoms['basins']
                )

            # tranform to a logarithmic scale
            da_KsatHorFrac.values = np.log10(da_KsatHorFrac.values)

            # resample KsatHorFrac with grid_from_rasterdataset
            resampled_KsatHorFrac = hydromt.workflows.grid.grid_from_rasterdataset(
                grid_like = mod.staticmaps['KsatHorFrac'],
                ds = da_KsatHorFrac,
                reproject_method = ['average'])

            # tranform back the log KsatHorFrac
            resampled_KsatHorFrac[var_name] = 10**resampled_KsatHorFrac[var_name]

            # fill missing values
            filled_KsatHorFrac = resampled_KsatHorFrac
            filled_KsatHorFrac[var_name] = filled_KsatHorFrac[var_name].interpolate_na(dim = 'lon', method='linear')

            # fill NA's around the boundary
            filled_KsatHorFrac[var_name] = filled_KsatHorFrac[var_name].interpolate_na(dim = 'lon', method='nearest', fill_value="extrapolate")

            # add and write to staticmap
            mod.set_staticmaps(filled_KsatHorFrac[var_name], 
                            name="KsatHorFrac_"+pred_method)

    # Edit model config
    setting_toml = {}
    if pred_method == "RF":
        setting_toml["input.lateral.subsurface.ksathorfrac"] = "KsatHorFrac_RF"
    elif pred_method == "BRT":
        setting_toml["input.lateral.subsurface.ksathorfrac"] = "KsatHorFrac_BRT"

    # Loop through each setting defined in setting_toml and update 
    # it in the model configuration
    for option in setting_toml:
        mod.set_config(option, setting_toml[option])

    #Write staticmaps and new TOML config
    if pred_method is not None:
        mod.write_staticmaps()
    mod.write_config(config_name="wflow_sbm_calibrated.toml")

    return None

[JoostBuitink]

Thanks for addressing this @RubenImhoff!

I think it would be best if we create a new setup function (setup_ksathorfrac or something?). I think we can largely include the code in your example to build this function. We also need to add this data to the data_catalog. As it is quite specific wflow_sbm related data, we need to check with the HydroMT-core team about how they feel adding this data to the deltares_data catalog.

[hboisgon]

Hi @RubenImhoff ! Nice work :) And indeed interesting to add to HydroMT. Would this function already do the trick: https://deltares.github.io/hydromt_wflow/latest/_generated/hydromt_wflow.WflowModel.setup_grid_from_raster.html Or you think that to help out users we can wrap it around a setup_ksathorfrac emthod that calls this function in the background?

Indeed next is to add the maps to a data catalog. This is quite wflow specific so we can decide to add to deltares data or wflow also has its own paramter_data catalog that is automatically read when building a model: https://github.com/Deltares/hydromt_wflow/blob/main/hydromt_wflow/data/parameters_data.yml

[hboisgon]

Ah I see in your function that you go to log before resampling and back to normal again. What is the reason for that? If needed then we definitely need to create a new function

[RubenImhoff]

Same as for the KsatVer parameter; the idea is to upscale from the log, because Ksat has a log-normal distribution.

[RubenImhoff]

It would indeed also fit well in the wflow parameter_data catalog. :)