Notebook

Example: Working with models in Python¶

The main feature of HydroMT is to facilitate the process of building and analyzing spatial geoscientific models with a focus on water system models. It does so by automating the workflow to go from raw data to a complete model instance which is ready to run and to analyse model results once the simulation has finished.

This notebook will explore how to work with HydroMT models in Python.

In [ ]:

# import hydromt and setup logging
import hydromt
from hydromt.log import setuplog

# other imports
import matplotlib.pyplot as plt
import geopandas as gpd

logger = setuplog("working with models", log_level=10)

Available models in HydroMT¶

To know which models are available within your active environment, you can use global MODELS variable in hydromt

In [ ]:

# generic model classes
print(f"Generic model classes: {hydromt.MODELS.generic}")
# model classes from external plugin
print(f"Model classes from plugins: {hydromt.MODELS.plugins}")

Here you may only see the generic models grid_model, vector_model and network_model. There is one more generic model within HydroMT mesh_model which is only available if the additional python mesh dependency xugrid is available in the activated environment.

Model components¶

HydroMT defines any model through the model-agnostic Model API based on several general model components and computational unit components, see Model API. Below is a scheme representing the Model API and general model classes available in HydroMT (without any plugin):

Let's discover how models are constructed within HydroMT and take the example of grid_model. We will first instantiate a GridModel object called mod. The api property helps us discover the available components and their type. You can do the same with any other HydroMT Model class or plugin (give it a try!).

In [ ]:

mod = hydromt.GridModel()
mod.api

Here you see all the general components from the Model class like config, geoms, forcing etc. as well as the GridModel specific computational unit grid. You can see that most components are dictionaries of either xarray DataArray or Datasets or of geopandas GeoDataFrame. For now we are starting from an empty model so all these components will be empty but here is how you can access them:

In [ ]:

print(type(mod.grid))
mod.grid

In [ ]:

mod.geoms

Model setup_* methods¶

To fill in our model components with data, HydroMT uses setup_ methods. These methods go from reading input data using the DataAdapter, transforming the data using workflows (e.g. reprojection, deriving model parameters, etc...) and adding the new model data to the right model component. An overview of available setup methods can be found in the API reference for the GridModel, VectorModel, and MeshModel

Note that these methods for the generic model classes are still quite limited. To get an idea of potential setup_ methods, checkout the model plugins

Let's have a look at some examples of the setup_ functions to start populating our model like setup_grid. This method parses the HydroMT region option to define the geographic region of interest and grid of the GridModel to build and once done adds region into the geoms component and grid mask into the grid component. You can check the required arguments in the docs.

Let's now setup a region for our model using for example a subbasin for any point in the Piave basin. We first initialize a GridModel instance in writing mode at a model root folder. Data is sourced from the artifact_data catalog.

In [ ]:

root = "tmp_grid_model_py"
mod = hydromt.GridModel(
    root=root,
    mode="w",
    data_libs=["artifact_data=v0.0.8", "data/vito_reclass.yml"],
    logger=logger,
)

In [ ]:

xy = [12.2051, 45.8331]
region = {"subbasin": xy, "uparea": 50}
mod.setup_grid(
    region=region,
    res=1000,
    crs="utm",
    hydrography_fn="merit_hydro",
    basin_index_fn="merit_hydro_index",
)
print(mod.geoms)
print(mod.grid)

In [ ]:

# Plot
fig = plt.figure(figsize=(5, 6))
ax = plt.subplot()
mod.region.boundary.plot(ax=ax)
# grid mask
mod.grid["mask"].plot(ax=ax)
# grid vector cells using hydromt.raster.vector_grid method
mod.grid["mask"].raster.vector_grid().boundary.plot(ax=ax, color="black", linewidth=0.1)
# the outlet point we used to derive region
gdf_xy = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=[xy[0]], y=[xy[1]]), crs=4326)
gdf_xy.to_crs(mod.crs).plot(ax=ax, markersize=40, c="red", zorder=2)

Similarly, we can also populate the config component using the setup_config method. For HydroMT config represents the configuration of the model kernel, e.g. the file that would fix your model kernel run settings or list of outputs etc. For most models, this is usually a text file (for example .yaml, .ini, .toml, .inp formats) that can be ordered in sections. Within HydroMT we then use the dictionary object to represent each header/setting/value.

Let's populate our config with some simple settings:

In [ ]:

config = {
    "header": {"setting": "value"},
    "timers": {"start": "2010-02-05", "end": "2010-02-15"},
}

mod.setup_config(**config)
mod.config

We can setup maps data using the setup_maps_from_raster and setup_maps_from_raster_reclass methods. Both methods add data to the maps component based on input raster data (RasterDataset type), but the second method additionally reclassifies the input data based on a reclassification table. The maps component gathers any raster input data without any requirements for a specific grid (CRS and resolution). It can contain, for example, direct model input data for models like Delft3D-FM that will interpolate input data on the fly to the model mesh, or auxiliary data that are not used by the model kernel but can be used by HydroMT to build the model (e.g. a gridded DEM), etc.

For models that require all their input data to be resampled to the exact computation grid (all raster at the same resolution and projection), then the input data would go into the grid component. The corresponding setup_grid_from_raster and setup_grid_from_raster_reclass functions for the grid components are also available.

But back to our example, let's add both a DEM map from the data source merit_hydro_1k and a manning roughness map based on reclassified landuse data from the vito dataset to our model grid object.

In [ ]:

mod.setup_grid_from_rasterdataset(
    raster_fn="merit_hydro_1k",
    variables="elevtn",
    fill_method=None,
    reproject_method="bilinear",
    rename={"elevtn": "DEM"},
)
mod.setup_grid_from_raster_reclass(
    raster_fn="vito",
    fill_method="nearest",
    reclass_table_fn="vito_reclass",  # Note: from local data catalog
    reclass_variables=["manning"],
    reproject_method="average",
)

In [ ]:

# check which maps are read
print(f"model grid: {list(mod.grid.data_vars)}")

mod.grid["manning"]

In [ ]:

# Plot
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))
mod.grid["DEM"].raster.mask_nodata().plot(ax=axes[0], cmap="terrain")
mod.region.boundary.plot(ax=axes[0], color="k")
gdf_xy = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=[xy[0]], y=[xy[1]]), crs=4326).to_crs(mod.crs)
gdf_xy.plot(ax=axes[0], markersize=40, c="red", zorder=2)
axes[0].set_title("Elevation")

mod.grid["manning"].plot(ax=axes[1], cmap="viridis")
mod.region.boundary.plot(ax=axes[1], color="k")
gdf_xy.plot(ax=axes[1], markersize=40, c="red", zorder=2)
axes[1].set_title("Manning roughness")

Model read & write methods¶

Once our model is filled up with data, we can then write it down using either the general write method or component specific write_ methods. Similarly, our model can be read back with the general read method or component specific ones.

Let's now write our model into a model root folder.

In [ ]:

mod.write(components=["grid", "geoms", "config"])

In [ ]:

# print MODEL_ROOT folder
import os


def print_dir(root):
    for path, _, files in os.walk(root):
        print(path)
        for name in files:
            if name.endswith(".xml"):
                continue
            print(f" - {name}")


print_dir(root)

And now let's read it back in a new GridModel instance:

In [ ]:

mod2 = hydromt.GridModel(root=root, mode="r", logger=logger)
mod2.read(components=["config", "geoms", "grid"])

In [ ]:

# check which grid are read
print(f"model grid: {list(mod2.grid.data_vars)}")

Building / updating a model with python¶

Using the same functionalities, it is also possible to build or update a model within python instead of using the command line, using the build and update methods. Let's see how we could rebuild our previous GridModel with the build method.

First let's start with writing a HydroMT build configuration (ini-file) with the GridModel (setup) methods we want to use.

In [ ]:

from hydromt.config import configread
from pprint import pprint

# Read the build configuration
config = configread("grid_model_build.yaml")
pprint(
    config,
)

And now let's build our model:

In [ ]:

# First we instantiate GridModel with the output folder and use the write mode (build from scratch)
root3 = "tmp_grid_model_py1"
mod3 = hydromt.GridModel(
    root=root3,
    mode="w",
    data_libs=["artifact_data=v0.0.8", "data/vito_reclass.yml"],
    logger=logger,
)

# Now let's build it with the config file
mod3.build(region=region, opt=config)

In [ ]:

print_dir(root3)

And check that the results are similar to our one-by-one setup earlier:

In [ ]:

mod3.config

In [ ]:

# Plot
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))
mod3.grid["elevtn"].raster.mask_nodata().plot(ax=axes[0], cmap="terrain")
mod3.region.boundary.plot(ax=axes[0], color="k")
axes[0].set_title("Elevation")

mod3.grid["manning"].plot(ax=axes[1], cmap="viridis")
mod3.region.boundary.plot(ax=axes[1], color="k")
axes[1].set_title("Manning roughness")