 Hello, my name is Hans van der Quest, lecturer at IHE Delft Institute for Water Education. In this video I'm going to explain how streams and catchments can be delineated from a digital elevation model using GIS. After this lecture you'll be able to define what is a catchment, to describe the generic GIS procedure to delineate streams and catchments from a digital elevation model, to explain how slopes are calculated in GIS from digital elevation models, to explain corrections that have to be made to use a DEM for hydrological applications, and to explain the Strahler order method. There is some confusion in the terminology of the word catchment in English. In British English the word catchment refers to the surface area for which all the water drains to one outlet. This is also called a drainage basin. The term watershed refers to the drainage divide, which is indicated with a red dotted line in the picture. In US English this is a bit different. Watershed refers to the area for which all the water drains to one outlet, and the red dotted line is called the drainage divide. A drainage basin, which is used in both British and US English, can be defined as an area of land where all the water that originates from different kinds of precipitation converges to one single outlet, and this can be another river, a bigger river, or it can be a lake, a reservoir, or a wetland, or the sea, or the ocean. Because a lot of GIS software is in US English you will often see the word watershed in the tools for what we would call catchments in British English. Here we see a catchment, this is the Vestra catchment in Belgium, and the size of a catchment depends much on the ruggedness of the area, the elevations. Here we exaggerated a bit the elevations, and so in British English we would call the area inside the red polygon the catchment, and we call the red line itself the watershed. So how do we derive streams and catchments from a digital elevation model using GIS? On this slide I'm going to show you the general workflow that's applicable to many GIS desktop software. In the next slide I'm going into detail of every step in the procedure. First we need to acquire the DEM tiles, then we need to mosaic the DEM tiles if our study area is spread over different tiles. We need to re-project the DEM, we need to subset, then if necessary interpolate the voids, we need to fill the sinks and remove the spikes, we need to burn in the stream network if we have a stream network available, we calculate the flow direction map, we derive the streams, define the outlet of the catchment and then derive the catchment, and in the end we can convert the dataset to the format that our hydrological model needs. The first step is to download the DEM, normally it comes in tiles and if we don't have access to local data because of the cost or the unavailability we can use global open access data as explained in another video. So most of the time we will use the SRTM 1 arc second global product, which is around 30 meters on the equator. There are also some other products, the SRTM Voidfield product, it's around 30 meters for the USA and 90 meters globally, and there's the experimental Aster global DEM or GDEM product which is also around 30 meters spatial resolution. You have to keep in mind that resolution is not the same as accuracy, these terms should not be confused. Not the higher the resolution, the higher the accuracy, it is a different concept so you can have a very high spatial resolution but then even the error per pixel can be also very high. And it's also not the case that you always need the highest resolution for each purpose. So for catchment modeling, generally the SRTM 1 arc second global product is fine, but if you want to plan a DEM in a river, you probably need very high spatial resolution. In both cases you need a good accuracy. You can download the SRTM DEM from the USGS Earth Explorer website or you can use a plugin in QGIS, the SRTM download plugin. Then usually your study area is not located nicely in the center of a tile but is spread over two or more tiles. So the next step is to mosaic the DEM tiles, which means that you simply merge them. The next step is to re-project the DEM. This is necessary because the global data sets are usually in the geographic coordinate system, latitude, longitude, and the units are in degrees. For DEM analysis we need to use a metric projection, so we need to re-project to another coordinate system. If we cannot use the local projection of a country because the area is too large, it's a trans-boundary catchment, then we can use another more global projection such as the UTM. Now I'm going to show you how slope is calculated and why it does matter that the units are in meters instead of degrees like we have in latitude-longitude coordinates. In this graph we see on the x-axis the distance and on the y-axis the elevation. The graph shows how the elevation changes with the distance. If we want to know the slope at this point, we know from mathematics that we have to draw the tangent, and then the slope can be calculated as delta z over delta x, and if we take the arc tangent then we have it in degrees. You can imagine with this equation if the x coordinates are in degrees and the z coordinates in meters, the software will still calculate the slope by applying this equation, but the result does not make sense. The z units need to be the same as the x and the y units. Now this is for our graph, but how do we calculate the slope in a raster? Therefore we need to use a so-called moving window because we need to calculate the slope to the surrounding pixels. So we consider the pixel in the middle of a 3x3 moving window and we assign the steepest slope to the pixel that we consider. Then we move the window to the next pixel and we calculate again the slopes to all the directions and we assign the steepest slope every time to the center pixel that we consider. So this is a focal operation that we discussed in another video. Now there's one problem with this method. You can think of what would happen to the first row and column and the last row and column. There's not enough information to calculate the slope in a kernel or in a window and that will result in no data at the borders. So if you want to calculate slopes of a study area, be sure that you make your study area a little bit larger than strictly your area of interest because otherwise you will lose the borders. The next step is to subset the DEM. Subsetting is a technical term for clipping. Remember that in the previous steps we mosaic the different tiles, but our study area is often not covering the whole mosaic, but only part of it. The consequence is that if we keep all the pixels of the mosaic in the calculations that the calculation times will be too long and sometimes even your computer will run out of memory. So it's better to subset the area to your study area, but it should also not be too small because then it will cut off the boundaries of your catchment in the delineation. So a good way is to understand a little bit the elevation differences in your area of interest to determine where the outlet is and to determine where the source areas are and where the divide could be by looking at the elevations. And then you make it a little bit larger to have a safety margin. The next step is to interpolate the voids. Voids are no data pixels in your DEM which can be a result of the acquisition procedure. DEM procedures don't handle well areas that are covered by snow or that are in the shade of other mountains and therefore we need to use the surrounding information to interpolate these areas. And this is then the result. We can never create the real data, but we can make it possible to continue with the data set after interpolation of the voids. The next step is to fill the sinks. We should not be confused with the voids. Voids are no data and sinks are depressions, artificial depressions in the landscape that are caused by the DEM creation process. So it's one or more cells which don't have a downstream cell around it. So the water is not able to escape to the outlet. And we can remove those using functionality in GIS software which is called fill sink algorithms. If the landscape contains real sinks such as lakes or other depressions they need to be added after the removal of the sinks and sinks are also called pits. Here we see a digital elevation model where water is trapped in the center cell which is much lower than the surrounding cells. So there is no possibility to drain towards the outlet. The fill sinks algorithm will either fill up those depressions or cut through. So here you see the side of elevation, a transect. And the outlet is downstream of the flow direction. And the algorithm can either increase the surface of the depression or decrease the surface after the depression so the water is routed to the outlet. Use its necessary to burn in the stream network to force the water to follow the real rivers. This can only be done of course when the river network layer exists. The next step is to calculate the flow direction. Most often we use the D8 algorithm. Similar to the slope algorithm it will look in a 3 by 3 window and will calculate the slope to all the directions around the cell that we consider. And instead of assigning the slope it will assign the direction of the steepest slope. An alternative is the D-infinite algorithm which uses not the eight surrounding cells but it will use continuous directions. Of course that is more calculation intensive and it will not always result in better results. Let's have a closer look how to calculate the flow directions using the D8 algorithm. On the left side we see the elevations. Each pixel has the elevation in meters and we are looking at the center cell and we need the information of the surrounding cells because it's a focal operation like with the slope calculation. The spatial resolution is 30 meters so the pixels are 30 by 30 meters. Now the question is which direction are we going to assign to the center pixel. We know that we can calculate the slope by evaluating the surrounding pixels of the center pixel by calculating the delta z over delta x, the drop divided by the distance. And for the D8 algorithm we assign then the direction of the steepest slope. So you need to look at this and select one of the direction values that are displayed on the right side of this slide. Let's do the calculation. If you would consider the cell to the south, so direction number 4, the calculation will be 67 minus 52 divided by the distance 30 equals 0.5. If you would consider the diagonal, the drop will be 67 minus 48 but the distance is not anymore 30 because we have to apply Pythagoras. The diagonal distance is longer than horizontal or vertical distance therefore we use 30 square root 2 and this results in 0.45. So the correct answer of the flow direction was number 4, it's the direction to the south. So we do that for every cell. And then we can link the lines that flow in the same direction and in this way we can construct the stream network. The map on the right side is called the stream link map. Now the stream links are for every pixel but not every pixel contributes to a river. Not all water that falls on a pixel becomes a river so we have two methods to derive the stream. We can look at the flow accumulation or we can look at the strahler order. Let's have a closer look at the flow accumulation. In this map we see the amount of precipitation in each cell. We can superimpose the stream link that has been derived in the previous step. The stream link determines based on the flow direction how the amount of precipitation is moving through the catchment. When the precipitation follows the stream link towards the outlet the flow is accumulating and we can define a threshold above which we consider the accumulated flow as being part of a river. Let's assume that we consider 10 units as the threshold in this case. So that means that everything that accumulates larger than 10 is considered a river and this is the way that the stream network can be derived. To determine this threshold value is a bit difficult and it needs a bit of trial and error. There is no rule of thumb and it depends on other properties of your catchment. So here on the left we see the example if we have 500 units as a threshold and on the right side we see 1000. So we see if we make the threshold lower it results in more tributaries and subcatchments than if we make it higher. The second method to derive streams is to use the strahler order. This method orders the reaches and we can set the threshold for which we consider it as a stream. It starts with the smallest ones and they are ordered with number one. When two of the same order join it increases the order so two of one join into order two. When two of order two join it becomes order three. However, when a lower order joins a higher order the order is not changed so here it remains three. Then two of order three become four and when we have order one joining again it stays order four. Based on this ordering we can determine a threshold value for which we consider the reaches as being part of a river. For example all the orders larger than three. To determine this threshold we need to do a calibration. We can use a satellite image or an existing map and see which threshold value fits best with the knowledge we have of the area. Of course this will not always be a perfect fit because we filled the DEM which makes the DEM a model and there can be a lot of human influences such as mining or urbanization or channels. The next step is to define the outflow point of our catchment. We need to define this on the delineated river so we cannot use a background map where a river has been defined in a different way because our DEM and the delineated streams are now part of a model. And if we don't define the outlet on the model it will not result in the catchment. Now what can we define as an outlet that is a location in a river where we have discharge measurements or the outlet of a tributary. So here we define the red cell as the outlet and the algorithm can then define the catchment based on the flow direction, the stream link. So the blue area is the catchment area and orange the drainage divide and all the water that falls within that area drains towards the outlet to the red cell. So you've just learned the general procedure for catchment and stream delineation which consists of downloading the DEM tiles to Mosa Icta tiles. If your study area is spread over multiple tiles then we project the DEM to subset the DEM, interpolate the voids if necessary, fill the things and remove the spikes, burn the stream network if you find it necessary and if you have a stream network available, calculate the flow direction map, derive the streams, define the outflow point and then derive the catchment. And in the end this normally feeds into your model so it needs to be converted to your model or your tool for which you want to use this data. There are some limitations to this method. The methods are based on elevation differences. So in areas that are flat this does not work. Also in areas where there's no gravity, natural gravity flow but are more human controlled this does not work. So for the Netherlands this doesn't work because we don't have much elevation difference and also the flow is very much human controlled in a large part of the country where we pump up the water to keep us dry. So keep that in mind when you use this method. Eventually you will use this data for hydrological modeling together with other GIS data. So the catchment delineation procedure provides very important inputs to catchment models. The delineation of the study area, the streams and the flow direction can be used.