4  Exercise 1: Delineation of the Zwalm catchment

4.1 Assignment

The delineation of a catchment based on a digital elevation model is often one of the first steps in a hydrological modelling study. The procedure is relatively straightforward, but the digital elevation model needs to be properly processed prior to the analysis. The pre-processing of the digital elevation model at least includes the removal of data gaps (i.e. spatial interpolation) and the removal of sinks (or depressions). Based on the processed DTM, the catchment can then be easily delineated using the simple principle that water always flows in the direction of the largest gradient. Different algorithms exist to do so (see course notes); here, only the gradient of every cell to its eight direct neighbours, the deterministic 8 (D8) algorithm, will be used.

Tip 4.1

The outlet of the Zwalm catchment is located in the vicinity of \(x\)=101673 m and \(y\)=17537 m

For this geospatial analysis, GRASS will be used, a powerful computational engine for raster, vector, and geospatial processing. It supports terrain and ecosystem modelling, hydrology, data management, and imagery processing. There are multiple ways of using GRASS. In this exercise, you can choose between 2 ways of interacting with GRASS:

  • QGIS With QGIS as Graphical User Inteface (GUI), all GRASS functionality can be accessed via the Procesing Toolbox
  • Python The Python GRASS interface allows programmatic access to GRASS functionality. For an in depth introduction, we refer to the documentation. In this exercise, we will use the grass.tools interface. The initial setup of the GRASS environment is provided in scripts/1_catchment_delineation.py, to which you can add your own code. All tools are accessible via the tools variable defined there. For example if you want ot use the r.import function, you would call tools.r_import(input="PATH_TO_FILE",output="OUTPUT_NAME", overwrite=True). Note that:
    • It is recommended to set the overwrite=True argument for all functions.
    • It is a general convention that the . defined in the GRASS function name (as found in the tool manual)is replaced by _ in the Python function name.
    • The output of a previous function (e.g. OUTPUT_NAME in the example above) can be used as input for a subsequent function. For example, if you want to use the output of r.import as input (called elevation) for a slope calculation with r.slope.aspect, you would call tools.r_slope_aspect(elevation="OUTPUT_NAME", slope="SLOPE", overwrite=True).
    • In scripts/helper_functions.py, two helper functions are defined: get_data_array_from_raster gives you an xarray.DataArray given the name (e.g. "OUTPUT_NAME") of a GRASS raster layer and plot_grass_raster makes a plot of a GRASS raster layer (given its name as again). Use these functions to check intermediate results.

Irrespective of the choice of interface, an overview of all GRASS functions can be found in tool manual. Specific tips for each interface will be denoted by their logo in the steps below.

  1. Start by loading the DTM (see Section 3.2). Check if it is clipped to a reasonable region of interest around the Zwalm river (bounding box of approximately \(x_{\mathrm{min}}=\) 98000 m, \(x_{\mathrm{max}}\) = 116000 m, \(y_{\mathrm{min}}\) = 160000 m, \(y_{\mathrm{max}}\) = 180000 m).
    • Python Use the r.import function to load in raster data.
  2. Spatially gap-fill the DTM (i.e. make sure there are no missing values).
    • Python Use the r.fillnulls function to fill in missing values.
    • QGIS Use the Fill nodata tool (from GDAL).
  3. Filter the digital elevation model for sinks (or depressions) and fill them. Use the r.fill.dir function.
  4. Calculate the flow direction and accumulation into each grid cell of the domain. Use the r.watershed function. Make sure that D8 is used as flow direction algorithm. Further more, ensure that the minimum size of an exterior watershed basins is set to 500. Keep both the accumulation (i.e. how many cells drain into each cell) and the flow direction map as outputs.
  5. Check if the accumulation map aligns with the actual river network of the Zwalm catchment from the Vlaamse Hydrografische Atlas (see Section 3.1). In the attributes of the shapefile, you can filter in "NAMEN" for "Zwalmbeek - Dorenbosbeek" to only withhold the main Zwalm river.
    • Python Use get_data_array_from_raster to get the accumulation map as an xarray.DataArray. Use geopandas.read_file to read in the river network shapefile.
  6. Find the point with the largest accumulation that is located within the Zwalm river network. Rasterise the river network to the same grid as the accumulation map to do so. This point is the outlet of the catchment.
    • Python Use the make_geocube function from the geocube package. Use the like argument.
    • QGIS Use the Rasterize (vector to raster) tool (from GDAL).
  7. Given the outlet determined in the previous step, delineate the catchment using r.water.outlet. Vectorize the resulting raster map of the catchment.
    • Python Convert the output raster to an xarray.DataArray using a helper function. Subsequently, the vectorize function from the geocube package can be applied.
    • QGIS Use the Polygonize (raster to vector) tool (from GDAL).
  8. Repeat step 7 for the coordinates provided in Tip 4.1. Compare the resulting catchment boundary with the one obtained in step 7 and with the provided catchment boundary from the Oppervlaktewaterlichamen en hun afstroomgebieden, 2022-2027 dataset (see Section 3.1). What are the differences between these catchment boundaries? What could be the reasons for these differences?
  9. Calculate the catchment area.

In your results and discussion below, make sure to include an overview figure of the Zwalm catchment including the DTM, river network and catchment boundary.

4.2 Results and discussion