Tutorial: Stream and Catchment Delineation

1. Introduction

One of the most important uses of GIS in hydrology is the delineation of catchments and streams. This lesson presents a generic workflow for stream and catchment delineation for areas where only open data is available. The Shuttle Radar Topography Mission (SRTM) 1 Arc Second DEM will be used. The workflow will be applied to the Rur catchment in Germany.

2. Theory

2.1. Digital Elevation Models in GIS

Watch this video on Digital Elevation Models in GIS.

2.2. Stream and Catchment Delineation

Watch this video on the theory of stream and catchment delineation.

3. Download DEM Tiles

For the Rur study area, we will download the tiles from the SRTM 1 Arc-Second global data set. Since the end of 2014 a 1-arc second global digital elevation model (approximately 30 meters at the equator) has been released as open data. Most parts of the world have been covered by this dataset ranging from 54 degrees south to 60 degrees north latitude. A description of the SRTM data products can be found here.

The following steps explain how to download the SRTM DEM tiles for the study area using the SRTM-Downloader plugin.

1. Start QGIS Desktop. We’ll start with a new project. In the Browser panel, add the folder with the data for this lesson to the Favorites as you learned in the previous lesson.

2. Drag the boundingbox.shp layer from the Browser panel to the map canvas.

This polygon defines the boundary of our initial analysis. Normally you have to create such a polygon yourself or use the map canvas coordinates. We need this to clip the DEM to the study area. Note that the projection of the project is now set to the projection of the boundingbox layer.

3. Install the SRTM-Downloader plugin from the Plugins Manager.

4. Click the SRTM-Downloader button in the toolbar.

6. Choose for Output-Path the folder where you're storing the files for this lesson. Use the  button to locate the folder. Don't keep the default temporary folder! You'll loose the data.

4. Reproject DEM

The DEM is in its original Lat/Lon Geographic Coordinate System (GCS) with datum WGS 84 (EPSG: 4326). It is not recommended to use a GCS for DEM analysis, because the z units (e.g. meters) are different than the x and y units (degrees). We need to choose a projection for our project. If the project covers one country we can choose a national projection. In our case, however, the project covers multiple countries (Germany, the Netherlands, and Belgium). Therefore we will use a global projection: UTM Zone 32 North, with WGS-84 as datum.

We can find the EPSG codes at http://www.spatialreference.org

1. Use the website to search for UTM 32N wgs 84. You can leave QGIS running and open a browser.

The figure below shows the result.

We need to take a look at the EPSG codes. We will use EPSG: 32632 throughout this project - clicking on it will provide more details. Note that this is the same projection the boundingbox layer has and was used for setting the projection of our project when we loaded the boundingbox layer as our first layer.

Now we are going to reproject the DEM from unprojected GCS (Lat/Lon WGS 84 - EPSG: 4326) to UTM Zone 32 North / WGS 84 (EPSG: 32632).

2. From the main menu choose Raster | Projections | Warp (Reproject).

3. In the Warp (reproject) window, choose the Project CRS: EPSG:32632 - WGS 84/UTM zone 32N as Target CRS.

4. Now complete the dialog:

The dialog should now look like the figure below.

Note the gdalwarp command that will be executed in the background. Under Advanced Parameters, we could have specified the extent of the output file based on the boundingbox layer. In that way the dem_mosaic layer is reprojected and clipped in one step.

5. Click Run to run the algorithm. After running, click Close to close the window.

The re-projected DEM now appears in the map canvas.

This is a good time to save your project.

5. Create a Subset of the DEM

In order to reduce the calculation time of the algorithms, we will subset (or clip) the raster layer to the boundingbox polygon.

1. Now, from the main menu select Raster | Extraction | Clip Raster by Mask Layer.

2. In the Clip Raster by Mask Layer dialog, choose dem_reprojected for the Input layer. For Mask Layer, choose boundingbox. Check the box to Keep resolution of input raster and keep the defaults for the other options. Choose dem_subset for the Clipped (mask) output and click Run. Click Close when done.

Often you'll not have a bounding box shapefile. In that case you can choose Raster | Extraction | Clip Raster by Extent. Then you can drag a box in the Map Canvas and use that for clipping. While using that, make sure your on-the-fly reprojection is similar to the layer that you want to clip, because the map canvas coordinates are used by the algorithm.

3. Now you can remove dem_reprojected from the layers list as we have done before for other layers that are no longer needed.

Watch this video to check the steps in this section:

5.1. Styling the DEM

By default, QGIS styles raster layers with the Singleband gray renderer. As you have learned in Lesson 3, the software does not know the raster data type (e.g. boolean, discrete or continuous) and how the user wants to present the information. To help interpret the results, it is good practice to intuitively style your layers.

The DEM is a continuous raster. Continuous rasters represent gradients and can therefore contain real numbers ((also called decimal numbers or floating point). Continuous rasters are styled in QGIS using ramps from the Singleband pseudocolor renderer in the Layer Styling panel.

1. Select the dem_subset layer and click  to open the Layer Styling panel (or press F7).

2. In the Layer Styling panel, choose Singleband pseudocolor from the drop-down list.

3. Click the arrow at Color ramp and choose Create New Color Ramp from the dropdown menu.

4. In the popup Color ramp type dialog choose Catalog: cpt-city from the drop-down list.

The cpt-city catalog has a lot of useful preset color ramps.

5. Choose Topography | Elevation. Note that cd-a and sd-a are also nice choices.

6. Click OK to close the dialog.

This gives us more intuitive colors in the DEM where we can clearly distinguish higher and lower areas (if you don’t see the colors applied, click Classify).

Now we will further improve the vizualisation.

8. Right-click on the dem_subset layer in the Layers panel and select Duplicate Layer.

This creates a copy of the dem_subset layer called dem_subset copy.

9. Uncheck the dem_subset layer, check the dem_subset copy layer and rename it to hillshade.

10. In the Layer Styling panel, which should still be open, make sure that the hillshade layer is now selected. In the drop-down list change Singleband pseudocolor to Hillshade.

Now the hillshade layer is vizualised with a shading.

• Which direction is the illumination coming from?
  
• Is this possible in reality?

Hillshade gives the best results with an artificial illumination in the northwest, which in reality can not exist in the Northern Hemisphere. If you move the dial in the Layer Styling panel to the southwest, you will see an inverted relief. Also note that there is a Resampling section. The default resampling method for both Zoomed in and Zoomed out is Nearest Neighbor. This method is fine for categorical data however, elevation is considered continuous data. You should therefore choose a Zoomed in resampling method of Bilinear and a Zoomed out resampling method of Cubic.

Next, we're going to blend the DEM with the hillshade layer.

9. Switch on the dem_subset layer by checking the box.

10. In the Layer Styling panel make sure the dem_subset layer is selected. In the Layer Rendering block of the panel, change the Blending mode to Multiply.

As you can see, blending gives a much nicer effect than transparency. With transparency the colors will fade. Now we can clearly see the elevation differences: the gradient from south to north and the valleys where we expect the streams.

Blending modes allow for more elegant rendering between GIS layers. They can be much more powerful than simply adjusting layer Opacity. Blending modes allow for effects where the full intensity of an underlying layer is still visible through the layer above. There are 13 blending modes available. See here for more info.

Watch this video to check the steps of this section:

6. Fill Sinks and Calculate Flow Direction

Raw, unprocessed DEMs have artifacts such as depressions. Artifacts are a result of the DEM acquisition process and must be removed before a DEM can be used for hydrological analysis, like catchment and stream delineation or hydrological modelling. There are several algorithms for filling sinks. Here we'll use the lddcreate tool from the PCRaster Tools plugin. This tool fills the sinks and creates a flow direction map (also called local drain direction map) at the same time. The resulting flow direction raster can be used in the next steps.

6.1. Install the PCRaster Tools plugin

Before you can install the PCRaster Tools plugin you need to install PCRaster.

On Windows you can install PCRaster with the OSGeo4W installer. Alternatively you can install PCRaster and QGIS in conda. Below are the steps for installation on Windows, which you should only do if you have missed the steps when following the installation instructions of QGIS for this course. Otherwise continue from Step 13.

1. Save your QGIS project and close QGIS
   
2. Run the OSGeo4W setup which comes with your QGIS version
   
3. Choose Advanced Install, click Next
4. Choose Install from Internet, click Next
5. Select the root install directory or keep the defaults, click Next
6. Select local package directory or keep the defaults, click Next
7. Select your internet connection, click Next
8. Choose one of the download sites, click Next
9. In the Select Packages window search for pcraster
10. Click the arrows icon to change from skip to a PCRaster version to install:
    
    
11. Click Next to complete the installation
12. Click Finished to close the OSGeo4W Setup wizard
    

6.2. Convert GeoTIFF to PCRaster Format

In order to use the PCRaster tools we need to convert raster to the PCRaster map format, which is a GDAL supported format. This is needed, because PCRaster is strict with data types.

PCRaster distinguishes the following raster data types:

• Boolean: cells are 0 (False) or 1 (True).
• Nominal: cells have positive integer values representing classes (discreet raster) without a specific order, for example land-use classes.
  
• Ordinal: cells have positive integer values representing classes (discreet raster) with a specific order, for example stream order.
• Scalar: cells have real values (decimal, positive, negative). This is used for continuous rasters.
• Directional: for rasters with a compass direction (e.g. aspect).
• LDD: local drain direction, a specific data type for flow direction rasters.

6.3. Calculate the Flow Direction

The next step is to fill the sinks (artificial depressions in the DEM) and calculate the flow direction raster. PCRaster does both in one step with the lddcreate tool. Check the documentation for the details.

1. In the Processing Toolbox go to PCRaster | Hydrological and material transport operations | lddcreate

2. In the Lddcreate dialog choose dem as DEM layer. Keep all other defaults so it will fill the sinks completely. Save the results to flowdirection.map.

3. Click Run. This can take a few minutes before the result appears in the map canvas. Click Close to close the dialog.

The lddcreate tool fills the sinks in the DEM and results in a local drain direction (flow direction) map. If you need the filled DEM, you can use the lddcreatedem tool.

In the next step we'll style the flow direction map.

6.4. Styling the Flow Direction Layer using a Circular Color Ramp

1. Open the Layer Styling panel for the flowdirection layer.

2. Change the renderer to Paletted/Unique values, because the flow directions are encoded in discrete values from 1 to 9. Click Classify to assign random colors to the cell values.

3. Right-click on flowdirection in the Layers panel and choose Create Raster Attribute Table.

4. In the New Attribute Table pop-up choose the default Managed by the data provider, which is the GDAL auxiliary XML format, and click OK.

Now you'll see the new raster attribute table.

5. Toggle on the editing of the attribute table by clicking  .

6. Edit the values in the Class field to reflect the compass directions as text.

7. Change the colors in the Color field by clicking on the colors and using the RGB values from this figure:

8. Toggle off the editing and save the changes. Click Classify to apply the changes and click Close to close the attribute table.

When you blend the flowdirection layer with the hillshade layer the result looks like this:

6.5. Styling the Flow Direction Layer using Arrows

We can further improve the styling of the flowdirection layer by adding arrows. This can be done using the mesh styling functionality of QGIS. To use that functionality, we need to convert the PCRaster LDD to a mesh format. We can do that with the Crayfish plugin.

1. Install the Crayfish plugin from the Plugins Manager.

2. In the Processing Toolbox, go to Crayfish | Conversions | PCRaster LDD to GRIB.

3. In the PCRaster LDD to GRIB dialog, choose flowdirection as Input raster and flowdirmesh.grb as Output file (GRIB).

4. Click Run and Close to close the dialog.

5. In the Browser panel, expand the flowdirmesh.grb group (you might need to refresh the Browser panel with the  button) and drag the flowdirmesh layer with the mesh icon  to the map canvas.

This might take some time. If the file is too large for your computer’s memory, you can get errors. In that case, you can clip the flowdirection layer to a smaller area and repeat the steps to convert the file to the mesh format.

When the map canvas shows a completely yellow layer, the flowdirmesh layer has been loaded and we can start styling it.

6. Select the flowdirmesh layer in the Layers panel and open the Layer Styling panel.

7. In the Layer Styling panel, go to the Datasets tab and click on  to disable contours and click on the arrow to enable vectors.

Now so many arrows are drawn in the map canvas that it turns black. Let’s tune the settings to improve this.

8. Go to the Vectors tab and change the Arrow Length settings to Fixed and the Length to 2.00. Change the Color to dark blue.

9. Zoom in to see the flow direction with arrows. 

10. Check the box to Display on User Grid to show the arrows fixed to a grid, e.g. with an X and Y Spacing of 10 px.

The settings are depending on your zoom level. Play with the settings to get a nice result. You can also try the other Symbology settings for visualization as Streamlines and Traces.

Watch this video to check the steps until this point:

6.6. Visualize Flow Direction in 3D

We can also visualize the flow direction using the QGIS 3D Map View.

1. In the main menu, choose View | New 3D Map View.

2. In the 3D Map view, click the Options button and choose Configure....

3. In the 3D Configuration dialog, stay in the Terrain tab and change the Type to DEM (Raster Layer) and choose dem as Elevation. Click OK to apply the settings and return to the 3D Map view.

The 3D view will now start rendering. Try to navigate the scene with the different mouse buttons and the compass. You can change the Vertical scale, Tile resolution and Skirt height in the 3D Configuration to improve the visualization.

7. Delineate Streams

Now we have the flow direction we can proceed with delineating the streams.

First we'll calculate the Strahler Orders, then we'll calibrate the Strahler threshold to determine the streams.

7.1. Calculate Strahler Orders

Before we can derive the streams from the DEM, we need to determine what we consider streams. For this purpose we use the Strahler order. The higher the order, the bigger the stream.

1. Close the 3D Map view and show only dem_subset with blended hillshade by unchecking the other layers in the Layers panel. Also zoom out to the full extent.

2. Search for streamorder in the Processing Toolbox and select PCRaster | Hydrological and transport operations | streamorder.

3. In the Streamorder dialog, select flowdirection for the Local Drain Direction layer, use strahler.map as the Stream Order output filename, and click Run. Click Close when the algorithm is done.

To make more sense of the strahler layer we are going to style it.

• Is the Strahler order layer a Boolean, discrete, or continuous raster?

The Strahler order layer is a discrete raster, but the order of the classes is important. Therefore, for PCRaster, it has the ordinal data type. For discrete rasters in QGIS we use the Paletted/Unique values styling. The higher the Strahler order, the bigger the stream. So we will use a color ramp from white to blue.

4. Open the Layer Styling panel if you closed it before and make sure that the strahler layer is selected.

5. Choose Paletted/Unique values from the drop-down menu.

6. Click Classify.

This assigns a unique, random color to each unique value in the raster.

7. Right-click on Random colors and choose Blues.

The strahler layer now shows in an intuitive way that the higher orders will be larger streams than the lower orders.

Watch this video to check the steps in this section:

7.2. Calibrate Strahler Threshold to Determine Streams

7. Click on the color for value 1 and change it to dark blue. Select the two other classes with the Ctrl button pushed. Then click the Remove selected row(s) button  to remove these classes.

7.3. Calculate Channel Network

The next step is to calculate the channel network.

First we're going to create a layer with the Strahler orders for the rivers by selecting the cells with order equal to or higher than the threshold determined in the previous step. Here we'll use order 8 as threshold.

With the PCRaster Tools plugin, we can only use raster layers in map algebra. If we want to calculate a new boolean layer with strahler >= 8 as True (our channels) and < 8 as False, we need to first create an ordinal raster that consists of cells with only value 8. This can then be used in map algebra.

1. In the Processing Toolbox go to PCRaster | Data management | spatial

2. In the Spatial dialogue type 8 for Input nonspatial. Choose Ordinal for Output data type and strahler as the Mask layer. Call the Output raster layer ordinal8.map.

3. Click Run and Close the dialogue after completion.

Now we're going to create a boolean layer with 1  (True) for all strahler cells that are >= 8:

4. In the Processing Toolbox go to PCRaster | Conditional and boolean operators | comparison operators.

5. In the Comparison Operators dialogue choose strahler as Input raster, >= as Comparison operator and ordinal8 as the second Input raster.

6. Click Run and Close the dialogue after completion.

Now we have a boolean layer with the channels that we can use now to assign the river Strahler orders.

7. In the Processing Toolbox go to PCRaster | Conditional and boolean operators | ifthen.

8. In the If Then dialog, choose channels as the Input Boolean Condition Raster and strahler as the Input True Raster. Call the Output raster channelsstrahler.map.

This means: if the channels layer has cells that are True (1), then give those cells the value of the strahler layer. All other cells get "nodata".

9. Click Run and Close the dialog after processing.

Now we have a raster with the Strahler orders for the channels. This can be visualized better as a vector, which will be explained in the next section.

7.4. Convert Raster to Vector Lines

For better visualisation of the channel network, we need to convert the channelsstrahler layer to a line vector.

We'll use the GRASS tools from the Processing Toolbox for this.

First we need to thin the raster lines, so they are only 1 cell wide.

1. From the Processing Toolbox choose GRASS | Raster (r.*) | r.thin

2. In the r.thin dialogue choose channelsstrahler as Input raster layer to thin. Keep the defaults and save the Thinned raster as channelsthin.tif. Make sure to choose the GeoTIFF format instead of the PCRaster .map format.

3. Click Run to perform the thinning and click Close when the processing is completed.

4. Compare the result with the channelsstrahler layer to understand what r.thin does.

Now we can convert the channelsthin layer to a line vector.

5. In the Processing Toolbox go to GRASS | Raster (r.*) | r.to.vect

6. In the r.to.vect dialog

• choose channelsthin as Input raster layer
• choose line from the drop down menu as Feature type
  
• check the box to Use raster values as categories instead of a unique sequence
• under Advanced Parameters change the v.out.ogr output type to line
• name the Vectorized output channels.shp
  

7. Keep the rest as default. Click Run to perform the conversion. Click Close when the processing is done.

The result has some geometrical errors that can be fixed using v.clean from the Processing Toolbox or the Topology Checker plugin. This is not covered in this tutorial.

In the next section we'll style the channels layer.

7.5. Styling the channels

Now we can style the channels vector layer to get a better understanding of the results. We will apply a styling in such a way that the lines get thicker with higher Strahler orders.

2. Click the Field calculator button

3. In the Field Calculator dialog, create a new field with the name ORDER. It should have an Output field type of Whole number (integer) and Output field length 1. Under Expression write the code as given below:

The CASE expression is used to evaluate a series of conditions. In our case, that means that it first checks if cat is equal to 8. If that’s the case, it assigns value 1, etc. The CASE expression ends with the END statement.

4. Click OK to apply the calculation.

5. Check the result in the attribute table. Toggle off editing by clicking and and confirm that you want to save the edits. You can now close the attribute table.

6. Open the Layer Styling panel and set the target layer to channels.

7. Keep the Single Symbol renderer, but change the Symbol layer type to Interpolated Line.

8. Change the Stroke Width to Varying Width. Set the Start value and End value to the ORDER field. Click the refresh button to get the correct Min. value and Max. value. Change the Min. width to 0.30 Points and the Max. width to 1.00 Points.

9. Under Color, keep the Single Color but change it to blue with an RGB value of 15 |66 | 220.

Watch this video to check the steps:

8. Define Outflow Point

A catchment is an extent or an area of land where surface water from rain, melting snow, or ice converges to a single point at a lower elevation, usually the exit of the basin, where the waters join another water body, such as a river, lake, reservoir, estuary, wetland, sea, or ocean. In order to delineate a catchment we need to have:

• the coordinates of our outlet in the same coordinate system as the map we are using
• the channel network that matches the flow directions as calculated from a hydrologically correct DEM

The outflow point of the Rur catchment is in Roermond, where the Rur enters the Meuse river (Maas in Dutch). The channel network that has been derived in the previous step is in the channels layer. We will however use the channelsstrahler raster layer, because we need to define the outlet exactly on a river pixel.

1. Make sure you have the channelsstrahler layer on top of the OSM Standard layer from the QuickMapServices plugin. You can style the channelsstrahler layer with a blue ramp, so that the main channel appears in dark blue.

2. Look for the location where the Rur river flows into the Meuse. Note that on the OpenStreetMap the Dutch names are used, because it lies in the Netherlands. Rur is spelled as Roer and Meuse is spelled as Maas.

Note that the delineated channels are not corresponding well with the channels on OpenStreetMap. This can be for the following reasons: (1) Incorrect automatic delineation of streams, which can be caused by errors in the DEM or areas that are too flat, (2) Distortion due to (on-the-fly) reprojection and resampling, and (3) Human influence on the natural course of the channels. The catchment delineation, however, only works when the outlet is defined on a delineated channel, because that corresponds with our flow direction layer. Results can be improved by using different (parameters of) fill sinks algorithms or burning an existing vector layer of the stream network in the DEM.

3. Choose a pixel on the delineated channel that is close to the real outlet of the Rur in the Meuse (step 2).

4. Right-click on the pixel and copy the Map CRS coordinate.

5. Paste the coordinates in a text editor (e.g. Notepad) and add a comma and a 1:

We've added the ID 1 to the coordinate. You can add more outlets to the file with unique ID's if you want to delineate (sub)catchments for multiple outlets.

6. Save the file as outlet.txt.

7. In the Processing Toolbox go to PCRaster | Data management | Column file to PCRaster Map.

8. In the Column File to PCRaster Map dialog, browse to the text file with the outlet, choose flowdirection as the Raster mask layer, choose Nominal (small) as the Output data type and save the file as outlet.map.

9. Click Run and Close after processing.

Now the outlet point has been added to the map canvas.

10. style the point with a clear colour using the Paletted/unique values renderer.

Watch this video to check the steps:

In the next section we'll delineate the catchment of this outlet.

9. Delineate the Catchment

3. Click Run and Close the dialog after processing.

The result should look like the screenshot in the figure below if you zoom to the extent of the layer. The result is a nominalraster with a cell value 1 for the catchment and 0 for outside the catchment. Nominal, because our outlet raster was defined as nominal data type. If the outlet.txt file had more coordinates of outlets, each catchment would have had cells with the value of the corresponding outlet cell. If you have subcatchments with outlets upstream of a higher order catchment, you can use the PCRaster subcatchment tool to avoid the overlap.

In order to overlay the catchment boundary with other data, it is better to convert it from raster to vector (polygon).

4. To convert the raster layer to vector, go to the main menu and choose Raster | Conversion | Polygonize (Raster to vector).

5. Make sure you choose the right input and call the output Rur_catchment.shp. Click Run. Click Close to get back to the main screen.

6. Look at the result. Also check the attribute table (right-click on layer name and choose Open attribute table).

In the catchment calculation, cells belonging to the catchment get a value of 1, while the other cells get a value of 0. During the conversion to polygons, it can happen that geometry errors are introduced. If you find more than one feature with a value of 1 this indicates a geometry error (incorrect topology), because the boundary of the polygon makes a loop (see figure below). This can give errors when we use the polygon for geoprocessing. You can fix those errors with the fix geometry tool in the Processing Toolbox.

7. In the attribute table, toggle to editing mode using the button, then select the record that you want to remove by clicking on the row. The selection will be highlighted in yellow on the map. Click the button to delete the selected feature.

You might have noticed that the polygon has a hole. This is caused by a huge open pit lignite mine, which results in its own subcatchment in our procedure. However, we would like to see the mine as part of the delineated catchment and will correct this.

8. Enable the Advanced Digitizing Toolbar that we have also used in Lesson 1 (right-click on a toolbar and check the box before the toolbar).

9. Click the Delete Ring button and click on the mine.

The hole has now disappeared.

10. Toggle off editing by clicking again, and save the changes.

11. ow remove all unnecessary layers from the layers list so that we have only channels, Rur_catchment, dem_subset, hillshade and OSM Standard (in that order).

Watch this video to check the steps in this section:

10. Clipping layers to the catchment boundary

1. Let’s first clip the channels vector layer to see only the streams that are inside the catchment. Go to the main menu and select Vector | Geoprocessing Tools | Clip...

2 .Fill in the dialog as in the figure below to use the catchment layer as a "cookie cutter" to clip the Input layer channels to the boundary of the Overlay layer Rur_catchment. Call the Clipped layer Rur_channels.shp. Click Run to run the tool. Click Close to return the main screen.

If you get an error related to the geometry, you can fix it. Go to the Processing Toolbox and choose Vector geometry | Fix geometries. The geometric errors are caused by the vectorisation process.

We can easily copy the styles from channels to Rur_channels.

4. Right-click on the channels layer and choose Styles | Copy Style | All Style Categories.

5. Now right-click on Rur_channels and choose  Styles | Paste Style | All Style Categories.

6. Remove channels from the layers list and make sure that Rur_channels layer is at the top of the list.

11. Styling the Resulting Catchment Area

To show the results of your analysis you can use a technique named Inverted Polygon Shapeburst Fills to focus attention on the Rur Catchment.

1. Open the Layer Styling panel by clicking the button. Set the target layer to Rur_Catchment.

2. Change from a Single symbol renderer to an Inverted polygons renderer. This renders the data as the inverse of its geometry, i.e. you'll style everything outside of the polygon. This creates a mask around the Rur valley.

3. Next, select the Simple Fill component. Then change the Symbol layer type to Shapeburst fill. In the Gradient colors section use the default Two color method. Change the first color to an RGB value of 135 | 135 | 135. Change the second color to white with an opacity of 65%.

4. In the Shading style section, click the Set distance option and set the distance to 4 mm and increase the Blur strength to 10.

5. Finally at the top of the Layer Styling Panel select the Shapeburst fill component.

7. Click the Add symbol layer button. Change the new Simple fill renderer to a Symbol layer type of Outline: simple line. Give it a Color of black and a Stroke width of 0.46 mm.

This gives us the nicely styled basin as shown in in the figure below.

Watch this video to check the steps:

12. Storing the Data in a GeoPackage

To keep the data together and enable easy distribution, it is good to save the layers as a GeoPackage.

1. In the Processing Toolbox look for the Package Layers tool.

2. In the Package Layers dialogue click the and select all layers. Note that these are only the vector layers.

3. Save it as Rur_data.gpkg and click Run and Close when it's done. By default it will also save the styles.

4. We can add the raster layers from the Browser panel. Simply drag the raster layers (in our case dem_subset) to the Rur_data.gpkg. You might need to refresh the Browser panel to see the GeoPackage.

The layers in our project still refer to the the GeoTIFFs and the shapefiles. Let's change that in such a way that they will be linked to the layers in our GeoPackage.

5. Select a layer in the Layers panel and right-click. From the context menu, select Change Data Source.... Select the corresponding layer in the GeoPackage.

6. Repeat this for all layers.

Now we can save the project to the Rur_data.gpkg with the correct paths.

7. In the main menu choose Project | Save to | GeoPackage...

8. In the pop-up, browse to the Rur_data.gpkg at Connection and name the Project Rur. Click OK to save the project in the GeoPackage.

The end result should now look like this:

Watch this video to check the steps of this section:

13. Conclusions

In this lesson you have learned to:

• find and download SRTM DEM tiles from the EarthExplorer website
• mosaic raster layers into a virtual raster
• reproject rasters
• create subsets of rasters
• fill sinks in a DEM
• calculate and style flow direction raster layers
• calculate Strahler orders
• delineate streams
• delineate catchments
• style layers for visualising catchments