Tip

For an interactive online version click here:

Build a model from Python

In this example a simple SFINCS compound flood model will be made, using the underlying Python functions of HydroMT-SFINCS to build a model.

The model is situated in Northern Italy, where a small selection of topography and bathymetry data has already been made available for you to try the examples.

import os
import numpy as np
import pandas as pd
import geopandas as gpd

from hydromt_sfincs import SfincsModel
from hydromt_sfincs import utils

This example shows how to build a SFINCS model containing elevation data and spatially varying roughness (together processed into subgrid tables), spatially varying infiltration and a local floodwall. In addition, multiple forcing conditions are set-up, and this is all done using Python scripting.

In case you want to adjust this example to build a SFINCS model anywhere else in the world, you will have to add your own datasets to HydroMT’s data catalog. For more info on that, check-out:

• Prepare data catalog

• Example: Prepare data catalog

• Example: Datasources

Steps followed in this notebook to build your SFINCS model:

1. Open SfincsModel class, set data library and output folder

2. Specify characteristics of the wanted grid

3. Load in wanted elevation datasets

4. Make mask of active and inactive cells

5. Update mask with water level and outflow boundary cells

6. Add spatially varying roughness data

7. Make subgrid derived tables

8. Add spatially varying infiltration data

9. Add water level time-series as forcing

10. Add an upstream discharge time-series as forcing

11. Add spatially varying rainfall data

12. Add weirfile

13. Add observation points

14. Show model

15. Save all files

Let’s get started!

1. Initialize SfincsModel class, set data library and output folder:

# Initialize SfincsModel Python class with the artifact data catalog which contains publically available data for North Italy
# we overwrite (mode='w+') the existing model in the root directory if it exists
sf = SfincsModel(data_libs=["artifact_data"], root="tmp_sfincs_compound", mode="w+")

2. Specify characteristics of the wanted grid and generate grid:

For more info about how to define a grid, click here.

Hint: if you only have a bounding box or geometry, you can also use ``SfincsModel.setup_grid_from_region``.

# Specify an input dictionary with the grid settings x0,y0,dx,dy,nmax,mmax,rotation and epsg code.
# create SFINCS model with regular grid and characteristics of the input dictionary:
sf.setup_grid(
    x0=318650,
    y0=5040000,
    dx=50.0,
    dy=50.0,
    nmax=107,
    mmax=250,
    rotation=27,
    epsg=32633,
)

# the input file is automatically updated. Uncomment to displayed below:
print(sf.config)

{'mmax': 250, 'nmax': 107, 'dx': 50.0, 'dy': 50.0, 'x0': 318650, 'y0': 5040000, 'rotation': 27, 'latitude': 0.0, 'tref': datetime.datetime(2010, 2, 1, 0, 0), 'tstart': datetime.datetime(2010, 2, 1, 0, 0), 'tstop': datetime.datetime(2010, 2, 2, 0, 0), 'tspinup': 60.0, 'dtout': 3600.0, 'dthisout': 600.0, 'dtrstout': 0.0, 'dtmaxout': 99999.0, 'trstout': -999.0, 'dtwnd': 1800.0, 'alpha': 0.5, 'theta': 1.0, 'huthresh': 0.01, 'manning': 0.04, 'manning_land': 0.04, 'manning_sea': 0.02, 'rgh_lev_land': 0.0, 'zsini': 0.0, 'qinf': 0.0, 'rhoa': 1.25, 'rhow': 1024.0, 'dtmax': 60.0, 'advection': 2, 'baro': 0, 'pavbnd': 0, 'gapres': 101200.0, 'stopdepth': 100.0, 'crsgeo': 0, 'btfilter': 60.0, 'viscosity': 1, 'inputformat': 'bin', 'outputformat': 'net', 'cdnrb': 3, 'cdwnd': [0.0, 28.0, 50.0], 'cdval': [0.001, 0.0025, 0.0015], 'epsg': 32633}

# show the model grid outline
# sf.region.boundary.plot(figsize=(6,6))
_ = sf.plot_basemap(plot_region=True, bmap="sat", zoomlevel=12)

3. Load in wanted elevation datasets:

# In this example we want to combine 2 elevation datasets, merit_hydro as elevation and gebco as bathymetry, in that order.

# NOTE: from the 1st dataset (merit_hydro) only elevation above ("zmin":0.001) meters is used;
# the 2nd elevation dataset (gebco) is used where the 1st dataset returned nodata values
datasets_dep = [{"elevtn": "merit_hydro", "zmin": 0.001}, {"elevtn": "gebco"}]

# Add depth information to modelgrid based on these chosen datasets
dep = sf.setup_dep(datasets_dep=datasets_dep)

# Make a plot of the merged topobathy, here colour limits are set between an elevation of -5 to 5 meters
_ = sf.plot_basemap(variable="dep", bmap="sat", zoomlevel=12)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

4. Make mask of active and inactive cells:

Here we generate the mask of active (msk=1) and inactive cells (msk=0), determining what cells on your grid should be used. For more info about the msk-file, click here.

# Choosing how to choose you active cells can be based on multiple criteria, here we only specify a minimum elevation of -5 meters
sf.setup_mask_active(zmin=-5, reset_mask=True)

# Make a plot of the mask file
_ = sf.plot_basemap(variable="msk", plot_bounds=True, bmap="sat", zoomlevel=12)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

NOTE: - The given output of HydroMT says “3 gaps outside valid elevation range < 10 km2”. HydroMT does some smart filtering that if small groups of inactive cells are found, surrounded by active cells, these are still included as active, in this case 3 gaps. You can control the size of these gaps to filter by adding fill_area = 10 in setup_mask_active(). - A similar argument exists to neglect a group of active cells surrounded by inactive cells: drop_area - reset_mask=True means that every time you start with an ampy mask while defining your msk cells, if reset_bounds=False (default) you continue on the prevous time you ran that function

5. Update mask with water level and outflow boundary cells - including use of polygons:

# Loading a shapefile clicked by user:
file_name = "data//compound_example_outflow_boundary_polygon.geojson"
gdf_include = sf.data_catalog.get_geodataframe(file_name)

# Example of the same, but how to load an existing ascii .pol file with x&y-coordinates, e.g. coming out of Delft Dashboard, here assumed to be in the CRS of the SFINCS model:
# file_name = "XX.pol"
# gdf_include = utils.polygon2gdf(feats=utils.read_geoms(fn=file_name), crs=sf.crs)

In SFINCS you can specify cells where you want to force offshore water levels (msk=2), or outflow boundaries (msk=3)

# Here we add water level cells along the coastal boundary, for cells up to an elevation of -5 meters
sf.setup_mask_bounds(btype="waterlevel", zmax=-5, reset_bounds=True)

# Here we add outflow cells, only where clicked in shapefile along part of the lateral boundaries
sf.setup_mask_bounds(btype="outflow", include_mask=gdf_include, reset_bounds=True)

# Make a plot of the mask file
fig, ax = sf.plot_basemap(variable="msk", plot_bounds=True, bmap="sat", zoomlevel=12)
gdf_include.to_crs(sf.crs).boundary.plot(ax=ax, color="k", lw=1, ls="--")

<GeoAxesSubplot:title={'center':'SFINCS msk map'}, xlabel='x coordinate UTM zone 33N [m]', ylabel='y coordinate UTM zone 33N [m]'>

NOTE: - As you can see now, also msk=2 values (red line) have been added along the coastal boundary - As you can see now, also msk=3 values (purple line) have been added along the lateral inland boundaries within the gdf_include shapefile - reset_bounds=True means that you start without initial boundary cells (of the specified type), if restet_bounds=False (default) you build on the existing boundary cells (if available)

6. Add spatially varying roughness data:

# --> this is used in making subgrid tables

# read river shapefile and add manning value to the attributes
gdf = sf.data_catalog.get_rasterdataset("rivers_lin2019_v1", geom=sf.region).to_crs(
    sf.crs
)
gdf["geometry"] = gdf.buffer(50)
gdf["manning"] = 0.03

# rasterize the manning value of gdf to the  model grid
da_manning = sf.grid.raster.rasterize(gdf, "manning", nodata=np.nan)

# uncomment to plot either the raster or the vector data
# da_manning.plot(vmin=0, x='xc', y='yc', cmap='viridis')
# gdf.plot()

# use the river manning raster in combination with vito land to derive the manning roughness file
# NOTE that we can combine in-memory data with data from the data catalog
datasets_rgh = [{"manning": da_manning}, {"lulc": "vito"}]

sf.setup_manning_roughness(
    datasets_rgh=datasets_rgh,
    manning_land=0.04,
    manning_sea=0.02,
    rgh_lev_land=0,  # the minimum elevation of the land
)
_ = sf.plot_basemap(variable="manning", plot_bounds=False, bmap="sat", zoomlevel=12)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

7. Make subgrid derived tables:

You can specify multiple settings about how the subgrid derived tables should be made.

Every single grid cell of the flux grid of the size inp.dx by inp.dy is defined into subgrid pixels (default nr_subgrid_pixels = 20). For every subgrid pixel the topobathy data is loaded, ideally this consists of high-resolution DEM datasets that you specify as user.

In this example with dx=dy=50m, having nr_subgrid_pixels = 20 means we are loading data onto a 2.5 m subpixel grid However, the input data of Gebco and Merit_hydro is way coarser, therefore let’s set the ratio to 5 for now.

# Every single grid cell of the flux grid of the size inp.dx by inp.dy is defined into subgrid pixels (default is 20, nr_subgrid_pixels = 20).
# For every subgrid pixel the topobathy data is loaded, ideally this consists also of high-resolution DEM datasets that you specify as user.

sf.setup_subgrid(
    datasets_dep=datasets_dep,
    datasets_rgh=datasets_rgh,
    nr_subgrid_pixels=5,
    write_dep_tif=True,
    write_man_tif=False,
)

# NOTE: we turned on that the merged topobathy of the different (high-res) datasets is written to a geotiff

# NOTE: if you have a very large domain with 100,000s to millions of cells, and very high-resolution datasets, this step might take minutes to hours!!!
#       But good news; when finished succesfully, you can very quickly run very accurate SFINCS simulations!
#       The whole point of the subgrid functionality of SFINCS is that by derived subgrid tables based on high res elevation data,
#       you either have more accurate results or run on a coarser grid resolution (= much faster) or both

Now we can see what kind of subgrid-derived variables are created:

# uncomment to see the subgrid table variales
# sf.subgrid

# we can plot the 2D subgrid variables
_ = sf.plot_basemap(
    variable="subgrid.z_zmin", plot_bounds=False, bmap="sat", zoomlevel=12
)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

8. Add spatially varying infiltration data:

# independent from subgrid files
# curve number infiltration based on global CN dataset
sf.setup_cn_infiltration("gcn250", antecedent_moisture="avg")

# check all variables in the sf.grid dataset
sf.grid.data_vars.keys()

KeysView(Data variables:
    dep      (y, x) float32 -12.0 -12.0 -12.0 -12.0 ... 0.3033 0.3213 0.351
    msk      (y, x) uint8 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 1 1 1 1 1 1 1 1 3 3 3 3
    manning  (y, x) float64 0.02 0.02 0.02 0.02 0.02 ... 0.037 0.037 0.037 0.037
    scs      (y, x) float64 dask.array<chunksize=(107, 250), meta=np.ndarray>)

By now we have made all basic SFINCS spatial layers to make the mskfile, infiltrationfile and subgridfiles, now we’re going to add some forcing…

9. Add water level time-series as forcing:

# Change period of model simulation time, specified in yyyymmdd HHMMSS --> simulation time here is 24 hours
sf.setup_config(
    **{
        "tref": "20100201 000000",
        "tstart": "20100205 000000",
        "tstop": "20100207 000000",
    }
)

print(sf.config)

{'mmax': 250, 'nmax': 107, 'dx': 50.0, 'dy': 50.0, 'x0': 318650, 'y0': 5040000, 'rotation': 27, 'latitude': 0.0, 'tref': '20100201 000000', 'tstart': '20100205 000000', 'tstop': '20100207 000000', 'tspinup': 60.0, 'dtout': 3600.0, 'dthisout': 600.0, 'dtrstout': 0.0, 'dtmaxout': 99999.0, 'trstout': -999.0, 'dtwnd': 1800.0, 'alpha': 0.5, 'theta': 1.0, 'huthresh': 0.01, 'manning': 0.04, 'zsini': 0.0, 'rhoa': 1.25, 'rhow': 1024.0, 'dtmax': 60.0, 'advection': 2, 'baro': 0, 'pavbnd': 0, 'gapres': 101200.0, 'stopdepth': 100.0, 'crsgeo': 0, 'btfilter': 60.0, 'viscosity': 1, 'inputformat': 'bin', 'outputformat': 'net', 'cdnrb': 3, 'cdwnd': [0.0, 28.0, 50.0], 'cdval': [0.001, 0.0025, 0.0015], 'epsg': 32633, 'mskfile': 'sfincs.msk', 'indexfile': 'sfincs.ind', 'sbgfile': 'sfincs.sbg', 'scsfile': 'sfincs.scs'}

a. Specify water level input locations:

For more info about what the bndfile is, click here

# Here we specify at what x&y-locations we have measured/modelled input water level data in the bndfile of SFINCS:

# x&y-locations in same coordinate reference system as the grid:
x = [319526, 329195]
y = [5041108, 5046243]

# add to Geopandas dataframe as needed by HydroMT
pnts = gpd.points_from_xy(x, y)
index = [1, 2]  # NOTE that the index should start at one
bnd = gpd.GeoDataFrame(index=index, geometry=pnts, crs=sf.crs)

# show what has been created:
display(bnd)

	geometry
1	POINT (319526.000 5041108.000)
2	POINT (329195.000 5046243.000)

b. Make up some time-series:

For more info about what the bzsfile is, click here

# Here we specify at what times we are providing water level input, and afterwards what the values are per input location:

# In this case we will provide 3 values (periods=3) between the start (tstart=20100201 000000) and the end (tstop=20100201 120000) of the simulation:
time = pd.date_range(
    start=utils.parse_datetime(sf.config["tstart"]),
    end=utils.parse_datetime(sf.config["tstop"]),
    periods=3,
)

# and the actual water levels, in this case for input location 1 a water level rising from 0 to 2 meters and back to 0:
bzs = [[0, 0.25], [0.75, 1.0], [0, 0.25]]

bzspd = pd.DataFrame(index=time, columns=index, data=bzs)
display(bzspd)

	1	2
2010-02-05	0.00	0.25
2010-02-06	0.75	1.00
2010-02-07	0.00	0.25

# Actually add it to the SFINCS model class:
sf.setup_waterlevel_forcing(timeseries=bzspd, locations=bnd)

# NOTE: the waterlevel forcing data is now stored in the sf.forcing dictionary
sf.forcing.keys()

dict_keys(['bzs'])

10. Add an upstream discharge time-series as forcing

a. specify discharge input locations: srcfile

For more info about what the srcfile is, click here

b. specify discharge time-series: disfile

For more info about what the disfile is, click here

# We follow exactly the same steps as for the water level forcing, but now specify 1 location where we specify discharges in m3/s
x = [321732]
y = [5047136]

# add to Geopandas dataframe as needed by HydroMT
pnts = gpd.points_from_xy(x, y)
index = [1]  # NOTE that the index should start at one
src = gpd.GeoDataFrame(index=index, geometry=pnts, crs=sf.crs)

time = pd.date_range(
    start=utils.parse_datetime(sf.config["tstart"]),
    end=utils.parse_datetime(sf.config["tstop"]),
    periods=3,
)

# and the actual water levels, in this case for input location 1 a water level rising from 0 to 2 meters and back to 0:
dis = [[2.0], [5.0], [2.0]]

dispd = pd.DataFrame(index=time, columns=index, data=dis)

# now we call the function setup_discharge_forcing
sf.setup_discharge_forcing(timeseries=dispd, locations=src)

# NOTE: the discharge forcing data is now stored in the sf.forcing dictionary
sf.forcing.keys()

dict_keys(['bzs', 'dis'])

In case you want to add other types of forcing, see the example notebook example_forcing.ipynb for other types. Or read more about this in the SFINCS manual

11. Add spatially varying rainfall data:

sf.setup_precip_forcing_from_grid(precip="era5_hourly", aggregate=False)

# NOTE: the precipitation forcing data is now stored in the sf.forcing dictionary
sf.forcing.keys()

dict_keys(['bzs', 'dis', 'precip'])

# Plot combined forcing time-series:
_ = sf.plot_forcing()

12. Add weirfile:

In SFINCS, a weirfile is often used to explicity account for line-element features such as dikes, dunes or floodwalls. Read more about structures in the SFINCS manual

# In this example specify a 'line' style shapefile for the location of the weir to be added
# NOTE: optional: dz argument - If provided, for weir structures the z value is calculated from the model elevation (dep) plus dz.
sf.setup_structures(
    structures=r"data/compound_example_weirfile_input.geojson",
    stype="weir",
    dz=None,
)

# NOTE: the observation points are now stored in the sf.geoms dictionary
sf.geoms.keys()

dict_keys(['region', 'weir'])

13. Add observation points

For more info about what the obsfile is, click here

# Loading a point shapefile clicked by user:
# NOTE: merge=True makes HydroMT merge the new observation points with already existing observation points (if present)
sf.setup_observation_points(
    locations=r"data/compound_example_observation_points.geojson", merge=True
)

# NOTE: the observation points are now stored in the sf.geoms dictionary
sf.geoms.keys()

dict_keys(['region', 'weir', 'obs'])

14. Show model

# Use predefined plotting function 'plot_basemap' to show your full SFINCS model setup
_ = sf.plot_basemap(fn_out="basemap.png" ,bmap="sat", zoomlevel=12)

15. Save all files

sf.write()  # write all

# Show created files in folder:
dir_list = os.listdir(sf.root)
print(dir_list)

['sfincs.dis', 'sfincs.sbg', 'gis', 'sfincs.bnd', 'subgrid', 'sfincs.scs', 'sfincs.inp', 'sfincs.obs', 'sfincs.ind', 'precip.nc', 'sfincs.src', 'figs', 'sfincs.weir', 'sfincs.msk', 'hydromt.log', 'sfincs.bzs']

Your basemap and forcing figures are saved in the folder ‘figs’, GIS files (tiff & geojson) of your model setup in ‘gis’ and merged elevation and manning roughness on subgrid resolution in ‘subgrid’.

Now you have made a model, you can progress to the notebook: run_sfincs_model.ipynb