 In this tutorial, you will learn how to do stream and catchment delineation using PCRaster Python. Here we are in the tutorials environment in the Anaconda prompt, and I'm looking at the files that come with this tutorial. And we see a bounding box shape file and four tiles that are downloaded from Earth Explorer, SRTM tiles. And that's all we need to create the streams and the catchment in this area. The first step is to mosaic the SRTM tiles, the four tiles together, and make sure that you have pointed here to the right folder where these files are stored. And we can use some libraries import globe and from OSGO import gdoll. And let's make a list of input files that are the TIFF files. So we want to only list the TIFF files and we do it like this. And you see here already in this explorer that input files is a list with the size of four with these files in it. You can also print the list and there we see the file names of the different TIFFs in the Python list. Now we're going to use that list to create a mosaic in a virtual file, a VRT file. A virtual file is very efficient in this case because it doesn't create the whole dataset again, but just describes how the tiles are connected. We use gdoll build VRT for that, with the default settings. It only needs the output file name, mosaic.vrt, and the list of input files. We need to add flush cache to write the layer to disk, so let's do it. mosaic equals gdoll.build VRT, and then the output file name mosaic.vrt as a string, and then the list with the input files. And then we do mosaic.flush cache. And then we will find the mosaic.vrt on our disk. Let's check that. Here it is, and we can open that in QGIS. And there it is. Next step is to reproject and subset the mosaic, so we have a smaller area to continue with. So we'll create a function to reproject and clip the mosaic using gdoll-warp. Here we see the function. We import gdoll from Osgeo. We have the function reproject and clip. Inputs are the input raster, the output raster, projection, the shapefile, and the resolution. So here we define the gdoll-warp options. The cutline dsname is the shapefile. That's what we use to crop to the cutline and to clip. The output format is a geotiff. This is the projection, and here the x and y dimensions of the pixels. So the pixel size. And here the out image, where we perform the warp with the input raster, the output raster, and the options. Then in the main script, the mosaic is the vrt file that we just created. The polygon is the bounding box. The epsg code also a string. The spatial resolution 30 meters. And the dmsubset. Let's save this. I'm going to save this in the data folder. And this will be clip reproject.py. Let's run it. And there it is. So let's see if the result is on our disk. And there it is. dmsubset.tiff. And we can open that also in QGIS. We do refresh. Here it is. If I remove the mosaic, we see that the subset is there. And we can check the projection epsg 32632. Let's change it also here in the project. There it is. So the next step is to convert it to the PC raster format. And we'll use the same procedure as we've done before. I'm going to make a new script. This one, convert to PC raster. And there we have the tiff file as an input. dm.map as an output. It's a flow32 and a scalar. And I'm going to save this. I'm going to call this tiff to PC raster. And let's run it. It has run. Let's see if we have it also on our disk. We can import PC raster. We can read the map. And then with Aguila, we can visualize it. And this is our dm. And it's exactly the same as the subset. So no difference. The next step is to calculate the flow direction. In the stream and catchment delineation procedure, we need to first remove the pits in the dm. Pits are pixels surrounded by only higher pixels. The catchment can only have one pit, the outlet. So the other pits have to be removed by a procedure called fill sinks. In PC raster, the lddcreate operation will both fill the dm and derive the flow direction. The lddcreate operation needs the dm as input and has four arguments to control thresholds for the filling algorithm. These can be found in the PC raster documentation. Here we want to remove all pits. So we set the thresholds to a very high number. So let's calculate the flow direction using lddcreate. It needs as an input the dm. And then we use these very high threshold values to make sure that all pits are filled. This will take a bit because it's a big file, so have some patience. When it's done, we can visualize it. If we zoom in, we can see the flow directions and we can query them. And pits are these black dots. They're mostly on the sides. Although in this procedure we can proceed with the flow direction map, you might want to create a dm that is filled. And then you can use the lddcreate dm operation for that. It uses the same arguments as lddcreate. So let's get that back with arrow up and change it a bit. dm filled, lddcreate dm. For the rest it's the same. And that will also take a bit, but that will make a hydrologically corrected dm. And we can report the files to this. The next step is to delineate the streams. And we are going to look at two methods. The first method is the strahler-order method and the second one is the flow-accumulation method. Let's start with the strahler-order. Strahler-orders can be calculated using stream-order from PC Raster. And the input is the flow direction. This will order the streams using the methods developed by Strahler. Let's have a look at the result. Here we see the result. It makes more sense if we change the color. We do it by double clicking and I want it low orders, yellow and then to more blue for higher orders. So the more blue, the higher the chance that it's a real river. And that's exactly the point. We need to calibrate these rivers in order to find out what is the threshold value for strahler from which we consider that it's a river. So we need to determine the value from which we consider these streams a river. The data type of the strahler-order map is ordinal and it starts at one. So we need to determine after which strahler-order we can consider the flow big enough to call it a river. We do that through calibration. Let's calculate maps with strahler-order 5 to the maximum and compare them with open-street map. First we need to determine the maximum strahler-order in this map using map-maximum. The map operations are the global operations in PC Raster. So let's define it. The maximum strahler-order equals map-maximum of strahler-orders. And then I can visualize the result. And it is 10. So what these global functions in PC Raster do is they write the result of the global function to each pixel. So each pixel here is value 10, but 10 is the highest strahler-order found in the strahler-order map. So the problem is that the result is a PC Raster map with for each pixel the maximum value of the map. If we want to iterate over the strahler-orders, we need to get the value. And we can use the cell-value operation. The cell-value operation needs the raster for which it has to give the cell-value as an input, as well as a row and column index number. Here all cells are the same, so we can simply use 0 for the index, referring to the first row and first column of the input raster. So maximum strahler-order tuple because there will be the result equals cell-value of maximum strahler-order and then with index values 0. And then we can print the result. And we see it's a tuple within the first element. It has 10 in the maximum value and the second one true, so that is a non-missing value. So the maximum strahler-order value is then maximum strahler-order tuple, element 0. Let's print that. And now we see that's value 10. So we can use that now in the iteration. That is 4 order in range from 1 to maximum strahler-order value, plus 1 because the range doesn't include the last value, only the first value of the range. Then do for each order stream equals if then strahler-orders a larger or equal to order, then give a boolean of 1 and write it to disk. Let's write it as stream and we concatenate a stream of the order and we have to add .map. So here in the loop we say for each order in this range from 1 to the maximum order value plus 1 because it's not included, do stream equals if the strahler-orders are larger equal to the order number, then give it boolean 1, otherwise no data, and then we write it to the disk for each order. So this will give us 10 maps to order 1 to 10. And we can then evaluate all these orders in QGIS. So I'm going to open the OpenStreetMap and I'm going to open strahler-orders and you can for example check order number 8, so I've added it here. And here I can compare if order 8 matches well the amount of rivers that I see on this map. So that's the calibration. We see that it fits quite well and that it seems to correspond more or less with the rivers that we see on the map. So you can tune that by playing around with that threshold value and comparing it with OpenStreetMap or with a satellite image. And that's the calibration, so for now we will just stick to strahler-order 8 or larger to be considered as a river. Then another method is to use the flow accumulation method to derive the rivers. And there we can use the AccuFlux function for flow accumulation. So flow accumulation equals AccuFlux and then the input is the flow direction, the LDD and we use one unit of water, all the pixels that has to be accumulated. And then we can visualize the result and it is on a logarithmic scale because there are many pixels with a low value and a few with a very high value, the rivers. So we can change it here to shifted logarithmic. That's needed because there are zeros. And let's also change this again to yellow to blue. And here also again the blue ones are the rivers. And then we need to determine a threshold and you do that in the same way. So you compare it with OpenStreetMap in QGIS and then you can determine the threshold. For example, we can say river flow if then the flow accumulation is larger than 9000. Let's say that that's what we determined with our calibration then give us Boolean1 else no data. And then you can look at the river flow map where true is the river and false does not exist and there's no data. So you can use that also in the calibration. Here we'll use the Strahler method. So stream8 is already at the disk so that's our river. So the next step is to delineate the catchment belonging to a selected outlet. First we need to identify an outlet on the delineated stream of the previous step. You can use the result from the Strahler order method or the flow accumulation method for finding the pixels that are part of the river. You can use QGIS or Aguila to find the coordinate. In Aguila you can do that using the crosshair tool. So I'm going to read the river using stream8 and I'm going to visualize it. So we said this was our river and then I can define an outlet on the river. Let me use this here. Make sure it's on the river. That it's really in a pixel. And then we can read the coordinates from here and pieceRuster comes with a tool called 2MAP which reads a text file in the format xCoordinate, space, yCoordinate, space, id. And that's what we're going to use here. So what we're going to do is to use these coordinates and write them into a file. So I use copyCon, location.text and then first the xCoordinate to 88934 and the yCoordinate 5, 6, 7, 5, 1, 5, 0 and then id1, ctrl z and then I can use call2MAP to find the data type nominal location.text, location.map and then the clone is stream8.map and here we can see that one record is read and that it's apparently worked. Let's check it. idlocation.map There should be one value which we can't easily find but that somewhere there we can check it and we see that the minimum maximum value is 1 so that's okay. So let's read the outlet location.map and then we can use this to derive the catchment using the catchment operation. So this catchment is the roar catchment so I'm going to call it roar catchment equals then catchment and it needs the flow direction and it needs the outlet and let's check the result. This is our catchment. There's a hole in it because there's a deep mine there that the water flows around. It has its own sub catchment. Sometimes you don't want to have the catchment of a specified outlet where you measure the discharge but you want to have all the catchments that you can find in the DEM. In that case we need to find all the pits and these pits are then considered as outlets so we can create outlets using the pit operation which will find all the pits in the flow direction map and there are a lot of them there they are and most of them are at the side of the map and we'll use these pits to create the catchments so let's call it catchments equals and then the catchment operator uses the flow direction and now we'll use outlets and that will derive all the catchments and here they are so just based on all the pits we find here all the catchments in the study area Most are at the boundary because these are boundary effects where other catchments extend beyond the study area Now that's all nice stepwise but you also would like to automate this so also here we can compile everything in one piece of code and that is done here so written in different functions a function to Mosaik a function to then reproject and clip a function to convert to PCRuster then a function to calculate the flow direction function to do the stream delineation and then to define the location with the outlet using code to map here we use it with a system call and now this is the main script where we'll just save this in the data folder so I'm gonna remove this the extension of the tiles diff the name of the Mosaik name of the bounding box the projection, the resolution then the name of the subset the name of the PCRuster file for the DEM the output of flow direction the threshold so this needs to be calibrated so therefore you can't do this automatically here you can put the outlet coordinates that you want to use define the clone and then it can calculate all these things visualize it and report it so let's save this and run it let's collect delineation and let's run it with our coordinates 288934 and 567515 so it's safe to delete all the files before you run the script so you only have the input files then there is no... so let's run it and there's the result