hydroMT provides a simple interface to model schematization and forcing data from which we can make beautiful plots:
staticmaps
attribute as a xarray.Dataset
staticgeoms
attribute as a geopandas.GeoDataFrame
. Note that in case of Wflow these are not used by the model engine, but only for analysis and visualization purposes.forcing
attribute as a dict
of xarray.DataArray
We rely on the cartopy package to plot maps. This packages provides a simple interface to plot geographic data and add background satellite imagery.
import pandas as pd
import xarray as xr
import numpy as np
from os.path import join, dirname
import os
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import cm, colors
import hydromt
root = 'wflow_piave_subbasin'
mod = hydromt.WflowModel(root, mode='r')
Here we plot the model basemaps (topography map with rivers, lakes, reservoirs, glaciers and gauges geometries).
# plot maps dependencies
import matplotlib.patches as mpatches
import cartopy.crs as ccrs
import descartes # required to plot polygons
import cartopy.io.img_tiles as cimgt
# read and mask the model elevation
da =mod.staticmaps['wflow_dem'].raster.mask_nodata()
da.attrs.update(long_name='elevation', units='m')
# read/derive river geometries
gdf_riv = mod.rivers
# read/derive model basin boundary
gdf_bas = mod.basins
plt.style.use('seaborn-whitegrid') # set nice style
# we assume the model maps are in the geographic CRS EPSG:4326
proj = ccrs.PlateCarree()
# adjust zoomlevel and figure size to your basis size & aspect
zoom_level = 10
figsize=(10, 8)
shaded= False # shaded elevation (looks nicer with more pixels (e.g.: larger basins))!
# initialize image with geoaxes
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(projection=proj)
extent = np.array(da.raster.box.buffer(0.02).total_bounds)[[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\
# create nice colormap
vmin, vmax = da.quantile([0.0, 0.98]).compute()
c_dem = plt.cm.terrain(np.linspace(0.25, 1, 256))
cmap = colors.LinearSegmentedColormap.from_list("dem", c_dem)
norm = colors.Normalize(vmin=vmin, vmax=vmax)
kwargs = dict(cmap=cmap, norm=norm)
# plot 'normal' elevation
da.plot(transform=proj, ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=.8), **kwargs)
# plot elevation with shades
if shaded:
ls = colors.LightSource(azdeg=315, altdeg=45)
dx, dy = da.raster.res
_rgb = ls.shade(
da.fillna(0).values,
norm=kwargs["norm"],
cmap=kwargs["cmap"],
blend_mode="soft",
dx=dx,
dy=dy,
vert_exag=200,
)
rgb = xr.DataArray(
dims=("y", "x", "rgb"), data=_rgb, coords=da.raster.coords
)
rgb = xr.where(np.isnan(da), np.nan, rgb)
rgb.plot.imshow(transform=proj, ax=ax, zorder=2)
# plot rivers with increasing width with stream order
kwargs = dict()
for strord in np.unique(gdf_riv['strord']):
if strord == np.unique(gdf_riv['strord']).max():
kwargs.update(label='river')
gdf_riv[gdf_riv['strord']==strord].plot(ax=ax, linewidth=strord/5, color='blue', zorder=3, **kwargs)
# plot the basin boundary
gdf_bas.boundary.plot(ax=ax, color='k', linewidth=0.3)
# plot various vector layers if present
if 'gauges' in mod.staticgeoms:
mod.staticgeoms['gauges'].plot(ax=ax, marker='d', markersize=25, facecolor='k', zorder=5, label='gauges')
patches = [] # manual patches for legend, see https://github.com/geopandas/geopandas/issues/660
if 'lakes' in mod.staticgeoms:
kwargs = dict(facecolor='lightblue', edgecolor='black', linewidth=1, label='lakes')
mod.staticgeoms['lakes'].plot(ax=ax, zorder=4, **kwargs)
patches.append(mpatches.Patch(**kwargs))
if 'reservoirs' in mod.staticgeoms:
kwargs = dict(facecolor='blue', edgecolor='black', linewidth=1, label='reservoirs')
mod.staticgeoms['reservoirs'].plot(ax=ax, zorder=4, **kwargs)
patches.append(mpatches.Patch(**kwargs))
if 'glaciers' in mod.staticgeoms:
kwargs = dict(facecolor='grey', edgecolor='grey', linewidth=1, label='glaciers')
mod.staticgeoms['glaciers'].plot(ax=ax, zorder=4, **kwargs)
patches.append(mpatches.Patch(**kwargs))
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"wflow base map")
legend = ax.legend(
handles=[*ax.get_legend_handles_labels()[0], *patches],
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")
Here we plot the model basin average forcing.
# read wflow forcing; mask region outside the basin and compute the basin average
# NOTE: only very limited forcing data is available from the artifacts
ds_forcing = xr.merge(mod.forcing.values()).where(mod.staticmaps['wflow_subcatch']>0)
ds_forcing = ds_forcing.mean(dim=[ds_forcing.raster.x_dim, ds_forcing.raster.y_dim])
# plot axes labels
_ATTRS = {
"precip": {"standard_name": "precipitation", "unit": "mm.day-1", "color": 'darkblue'},
"pet": {"standard_name": "potential evapotranspiration", "unit": "mm.day-1", "color": 'purple'},
"temp": {"standard_name": "temperature", "unit": "degree C", "color": 'orange'},
}
plt.style.use('seaborn') # set nice style
n = len(ds_forcing.data_vars)
kwargs0 = dict(sharex=True, figsize=(6, n * 3))
fig, axes = plt.subplots(n, 1, **kwargs0)
axes = [axes] if n == 1 else axes
for i, name in enumerate(ds_forcing.data_vars):
df = ds_forcing[name].squeeze().to_series()
attrs = _ATTRS[name]
longname = attrs.get("standard_name", "")
unit = attrs.get("unit", "")
if name == 'precip':
axes[i].bar(df.index, df.values, facecolor=attrs['color'])
else:
df.plot.line(ax=axes[i], x="time", color=attrs['color'])
axes[i].set_title(longname)
axes[i].set_ylabel(f'{longname}\n[{unit}]')
# save figure
# fn_out = join(mod.root, "figs", "forcing.png")
# plt.savefig(fn_out, dpi=225, bbox_inches="tight")