Burning river in DEM using map algebra in QGIS

This video demonstrates how to use map algebra for burning a river into a DEM prior to catchment delineation. In this way we force the drainage in a DEM to follow the river. In the presented method the DEM is proportionally lowered in order to preserve the flow direction in the DEM. The method has been described by Wesseling et al. (1997) and developed for PCRaster. Here QGIS is used.
First the river is rasterized. Then a standardized DEM is calculated with values between 0 and 1. Next at river cells the standardized DEM is reduced with 1, resulting in a range from -1 to 0 for those cells. Finally we multiply the standardized DEM with the original DEM to get the values back in the original range, but with the river burned into the grid.
Credits for OpenStreetMap data: © OpenStreetMap contributors
SRTM data available from USGS Earth Explorer: earthexplorer.usgs.gov
Reference to the method:
Wesseling, C. G., W. P. A. Van Deursen, and M. De Wit. "Large scale catchment delineation: a case study for the river Rhine basin." Proc. JEC-GI 97 (1997): 16-18.
This video is part of IHE Delft OpenCourseWare. More free resources can be found at www.gisopencourseware.org
Subscribe to the KZbin Channel for updates!

[@AbdullahSalimalig]

Niice one

[@Zoe-kw6we]

Hi Hans! Thanks for the great video. Unfortunately, the 'reclassify values' option doesn't show up in my toolbox. I do have SAGA, but it's not there. I tried achieving the same result with QGIS's 'reclassify by table', but get 0 & 241 as options when I choose byte as data output & two categories of 1 when I choose Float32 as data output. I am new to QGIS so mainly working through trial & error, but I think not having an understanding of what I'm doing doesn't help here. Do you know why the 'reclassify values' option does not show up? & could I work around it with the QGIS 'reclassify by table' option if I choose different settings?

[@HansvanderKwast]

The "reclassify values" is indeed no longer included. The "reclassify by table" tool can indeed be used instead. The data type should be Byte and you can use the same settings. Check with the Paletted/Unique values renderer if the result then has only value 1. I'll soon upload a new version of the video where we use the new PCRaster Tools plugin, which might be easier.

[@Zoe-kw6we]

@@HansvanderKwast Thanks a lot for your answer, I've worked it out! Looking forward to more videos, your tutorials are the best! :)

[@benmissimer]

@@HansvanderKwast looking forward to the PCRaster tutorial. Did that get generated yet by chance?

[@JuanManuelCotrino]

Hello, I'm using the carve function but the output raster generates no data values in the last section of the river. I've played with the parameters but always end up getting the same no data pixels. That region is a valley, and because of that the river has no define "axis" and it's surrounded by lots of swamps. I don't know if I'm getting this no value pixels because of the flatness of that zone. Do you know the reason for this? thank you very much in advance, I hope you can help me.

[@HansvanderKwast]

Not sure what the problem is. You can also try this method: kzbin.info/www/bejne/oHm3kIaqhceSpa8

[@josephrojaszamora9601]

Hello Hans, great method but I have a problem when I try to rasterize the rivers as the Error shows: Attempt to create 0x0 dataset is illegal, sizes must be larger than zero. I tried changing the pixel length but still the same issue, what should I do to overcome this problem?

[@Zoe-kw6we]

Have you solved it? I'm having the same issue!

[@josephrojaszamora9601]

@@Zoe-kw6we Hey, no I could not, I ended up using arcgis to rasterize the rivers and export it as a shape and following the resto of the process on qgis

[@Zoe-kw6we]

@@josephrojaszamora9601 Thanks for your reply! I hope I find a way to work it out as well

[@yoaslauren8085]

Thank you for the video. Can we burn the stream into DEM using PCRaster in Python? Do you have references about it if we can do it?

[@HansvanderKwast]

Check some examples in the Resource Sharing repository: github.com/jvdkwast/qgisrepository/tree/main/collections/qgis_pcrasteruserscripts/processing

[@max7164]

Thanks for the video. Why is there a need to go through the std dtm step? Why not just subtract a set value for river depth from the raw dtm? Thanks.

[@HansvanderKwast]

According to the original publication of this algorithm this is done to "maintain the lowest neighbour relation between cells except for those cells where one neighbour is in the river network and the other one is not."

[@max7164]

@@HansvanderKwast thanks for your response. There seems to be a lot of good content on the channel. I will check it out!

[@dibyajyotinath3071]

How to do this in ArcGIS ?

[@guofengwang2614]

Can I burn raster images to DEM, such as MNDWI images？

[@HansvanderKwast]

You can use the same approach as in the video to burn any raster in then DEM.

[@guofengwang2614]

I burned an NDVI image to the DEM, which resulted in negative values ​​for the corrected DEM.

[@HansvanderKwast]

Interesting. Why would you burn NDVI in a DEM? Also negative elevation is normal where I live below sea level.

[@guofengwang2614]

MNDWI image.

[@guofengwang2614]

@@HansvanderKwast The corrected DEM range is not the same as the original. There should be a problem in the step of generating demBurn. The Mndwi image is subtracted from stdDEM. Although the purpose of burning is achieved, this step leads to a different range of corrected DEM. Is there any good solution?

[@keithbjohnson1]

Dear Hans, I am deeply grateful for all the teaching you have done on QGIS. It has been a huge part of my learning!. The last step in which you multiply the burned standardized DEM with the original DEM would seem to introduce a non-linearity in the DEM values, and not give back the original DEM values outside the river burned pixels.. if stdBurnedDEM = (DEM - base)/(top - base) then multiplying DEM*stdBurnedDEM gives (DEM*DEM - DEM*base) / (top - base). Instead, as the final scaling step, why not apply the inverse of the scaling that was used to create the stdDEM? i.e., Final Full Range burned DEM = stdBurnedDEM*(top - base) + base.

[@HansvanderKwast]

Dear Keith. Good point. I'll check it and get back to you.

[@SuperDanisse]

@@HansvanderKwast did you check it?

[@smrutiranjanpatra3234]

I think this step should be rechecked. I find erroneous negative elevation values while simply multiplying the standardized burner dem with the original dem.

[@smrutiranjanpatra3234]

Any advice on this step?

[@HansvanderKwast]

Agree. I think that step is incorrect. In the original paper they've used the relative DEM directly to derive the flow direction.

[@salpiseyheng6089]

Dear sir, can we also burn road network in dem by this method? Thanks

[@HansvanderKwast]

With this method you can burn anything in your DEM.

[@salpiseyheng6089]

@@HansvanderKwast Dear sir, one more question, i saw you use rive line to burn in dem. why you not use river in polygon to burn because it can also identify the width of river? Thanks

[@HansvanderKwast]

Sure. That's better in case you have a polygon.

[@salpiseyheng6089]

@@HansvanderKwast sir, one more, we fill sink just for delineate watershed? how about we use it for hydraulic simulation like flood modelling in hec ras, we need to fill sinks for dem or no need? Thanks

[@HansvanderKwast]

For hydrological applications with water flow you need to fill the sinks, otherwise the water will not reach the outlet. Both artifacts and real sinks will be removed, so take into account that if you need the real depressions to stay in you elevation model.