Plot Delft3D FM model using HydroMT

HydroMT provides a simple interface to model schematization from which we can make beautiful plots:

• Mesh layers are saved to the model mesh.data component as a xugrid.UgridDataset

• Vector layers are saved to the model geoms.data component as a geopandas.GeoDataFrame. Note that in case of Delft3D FM these are not used by the model kernel, but only for analysis and visualization purposes.

• Gridded data like bedlevels or infiltration capacity are saved to the model inifield.data component as a xarray.DataArray. Here the maps are regular grid in the same CRS as the Delft3D FM model but not necessarily the same resolution or grid, as Delft3D FM kernel can do the interpolation.

We use the cartopy package to plot maps. This packages provides a simple interface to plot geographic data and add background satellite imagery.

Load dependencies

import numpy as np
from hydromt_delft3dfm import DFlowFMModel

/home/runner/work/hydromt_delft3dfm/hydromt_delft3dfm/.pixi/envs/default/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

# plot maps dependencies
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt

Read the model

root = "dflowfm_piave"
mod = DFlowFMModel(root, mode="r")

/home/runner/work/hydromt_delft3dfm/hydromt_delft3dfm/.pixi/envs/default/lib/python3.12/site-packages/xugrid/ugrid/conventions.py:252: UserWarning: No standard_name of ('projection_x_coordinate', 'longitude', 'projection_y_coordinate', 'latitude') in ['mesh1d_node_branch', 'mesh1d_node_offset'].
Using mesh1d_node_branch and mesh1d_node_offset as projected x and y coordinates.
  vardict[name] = _infer_xy_coords(ds, candidates)
/home/runner/work/hydromt_delft3dfm/hydromt_delft3dfm/.pixi/envs/default/lib/python3.12/site-packages/xugrid/ugrid/conventions.py:252: UserWarning: No standard_name of ('projection_x_coordinate', 'longitude', 'projection_y_coordinate', 'latitude') in ['mesh1d_edge_branch', 'mesh1d_edge_offset'].
Using mesh1d_edge_branch and mesh1d_edge_offset as projected x and y coordinates.
  vardict[name] = _infer_xy_coords(ds, candidates)
/home/runner/work/hydromt_delft3dfm/hydromt_delft3dfm/.pixi/envs/default/lib/python3.12/site-packages/xugrid/ugrid/ugrid1d.py:170: UserWarning: No CRS or recognizable standard_name found for topology 'mesh1d'. Assuming projected coordinates.
  crs, is_projected = cls._extract_crs(ds, topology)
Replacing geom: crosssections
Replacing geom: manholes
Replacing geom: rivers
Replacing geom: pipes
Replacing geom: branches

# read the mesh to get the grids as geodataframes for plotting
mesh1d = mod.mesh.mesh_gdf["mesh1d"]
mesh2d = mod.mesh.mesh_gdf["mesh2d"]

# Get the different types of branches in mesh1d
rivers = mod.rivers
pipes = mod.pipes
# Additional geometry and structures
manholes = mod.geoms.data["manholes"]
crosssections = mod.geoms.data["crosssections"]

# Read the elevation from maps
elv = mod.inifield.data["elevtn"]

Plot model schematization base maps

Here we plot the model base mesh information as well as the topography map.

# we assume the model maps are in the projected CRS EPSG:3857
proj = ccrs.epsg(3857)
# adjust zoomlevel and figure size to your basis size & aspect
zoom_level = 12
figsize = (10, 8)

# initialize image with geoaxes
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(projection=proj)
bbox = elv.raster.box.to_crs(3857).buffer(3e3).to_crs(elv.raster.crs).total_bounds
extent = np.array(bbox)[[0, 2, 1, 3]]
ax.set_extent(extent, crs=proj)

# add sat background image
ax.add_image(cimgt.QuadtreeTiles(), zoom_level, alpha=0.5)

## plot elevation\
elv.plot(transform=proj, ax=ax, zorder=0.5, cmap="BrBG", add_colorbar=True)

# plot rivers
rivers.plot(ax=ax, linewidth= 1, color="blue", zorder=3, label="rivers")
# plot pipes
pipes.plot(ax=ax, color="k", linewidth=1, zorder=3, label="pipes")

## plot mesh
mesh1d.plot(ax=ax, color="r", zorder=2, label="mesh1d")
mesh2d.plot(ax=ax, facecolor="none", edgecolor="gray", linewidth=0.5, zorder=2, label="mesh2d")

ax.xaxis.set_visible(True)
ax.yaxis.set_visible(True)
ax.set_ylabel(f"latitude [degree north]")
ax.set_xlabel(f"longitude [degree east]")
_ = ax.set_title(f"dflowfm base map")
legend = ax.legend(
    title="Legend",
    loc="lower right",
    frameon=True,
    framealpha=0.7,
    edgecolor="k",
    facecolor="white",
)

# save figure
# NOTE create figs folder in model root if it does not exist
# fn_out = join(mod.root, "figs", "basemap.png")
# plt.savefig(fn_out, dpi=225, bbox_inches="tight")

/tmp/ipykernel_3228/947199666.py:34: UserWarning: Legend does not support handles for PatchCollection instances.
See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler
  legend = ax.legend(

And zoom in to see better some of the structures:

# we assume the model maps are in the projected CRS EPSG:3857
proj = ccrs.epsg(3857)
# adjust zoomlevel and figure size to your basis size & aspect
zoom_level = 14
figsize = (10, 8)
#extent = np.array([1385000, 1390000, 5861000, 5863000])
extent = np.array([1385000, 1390000, 5857000, 5860000])

# initialize image with geoaxes
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(projection=proj)
ax.set_extent(extent, crs=proj)

# add sat background image
ax.add_image(cimgt.QuadtreeTiles(), zoom_level, alpha=0.5)

# plot rivers
rivers.plot(ax=ax, linewidth= 1, color="blue", zorder=2, label="rivers")
# plot pipes
pipes.plot(ax=ax, color="k", linewidth=1, zorder=2, label="pipes")

# plot structures
manholes.plot(ax=ax, facecolor="y", markersize=2, zorder=4, label="manholes")
crosssections.plot(ax=ax, facecolor="grey", marker = '|',  markersize=15, zorder=3, label="cross-sections")

ax.xaxis.set_visible(True)
ax.yaxis.set_visible(True)
ax.set_ylabel(f"latitude [degree north]")
ax.set_xlabel(f"longitude [degree east]")
_ = ax.set_title(f"dflowfm base map")
legend = ax.legend(
    title="Legend",
    loc="lower right",
    frameon=True,
    framealpha=0.7,
    edgecolor="k",
    facecolor="white",
)