 Hello, in this video we are going to delineate a catchment using PC Raster in Python. This video follows a tutorial that you can find on GISopenCourseWare.org and is part of the course Programming for Geospatial Hydrological Applications. You can find the link in the description of this video. The tutorial is provided in a Jupyter notebook, but because we want to visualize the results in QGIS we are going to run the instructions from the QGIS Python console. This will work with a PC Raster in Python in a Conda environment, but if you have installed a PC Raster for the PC Raster Tools plugin you can also run QGIS without a Conda environment and use a PC Raster that comes with QGIS. If you have installed it through the OSDL4W installer for example. So the PC Raster Tools plugin provides a graphical user interface to the PC Raster Python operations and you can easily install it in QGIS and have a look. There are other courses and videos available and this will add these tools to the processing toolbox and these tools are named exactly the same as the operations in Python and if you click on a tool you can always click on the link to the documentation. You can find there how to write the syntax and all the parameters are documented in their examples on how to write this in Python code. But for this tutorial we're going to use the Python console in QGIS. We can launch the Python console by clicking the Python icon and we can open an editor by clicking this icon. We can rearrange the screen a little bit and here you see a Python prompt as you also have in the Jupyter Lab or if you just type Python at the terminal. We can import packages or libraries from the prompt and the first thing we need to do is to mosaic the tiles from SRTN and therefore we need to use GLOB and we need to use GDAL and GDAL is from the OS Geo package and we also need to import the OS package because our current working directory is where QGIS started and we can get that information with os.getcwd and we see that we are in OS Geo for W bin while we need to be in our data folder so I'm going to use os.change directory and there we give in a string with forward slashes the location which is pzraster tutorials slash pzraster catchment delineation slash data. Let's check and indeed we are now in the correct folder and we're going to create a list of input files with the SRTM TIFF files. So using GLOB.GLOB we can search for file extensions or patterns and here we look for TIFF files. These are the four SRTM tiles that you can see and this will generate a list so if I print the input files variable the result is a list with each element the string with the file name of the Geo TIFF. Now we're going to create a virtual raster we call the variable mosaic and we use GDAL build vrt and it needs the destination name we call it mosaic dot vrt and that will be our output it's a string and as an input we use our list of files the input files variable. We type enter it runs but to save it to disk we need to use mosaic dot flush cache. You see I made a typo here because it's case sensitive so I need to type it exactly in a correct way where the C has to be a capital. Now when I look in the browser panel I can see the vrt file and I can drag it to the map canvas so here's our mosaic. This is very similar to what QGIS does through the graphical user interface with creating a virtual raster. GDAL is used in the back end. The next step is to reproject and clip the mosaic to a smaller area fitting our bounding box polygon and we use this little script here which uses a function reproject and clip which needs an input raster an output raster projection shape file and resolution and it uses GDAL warp which has GDAL warp options needs the cut line ds name which is our shape file crop to cut line is true the format that we want is a geotiff we add a projection and spatial resolution for the output and you see then under below the function that the mosaic file name will be mosaic dot vrt and the polygon is bounding box dot shp that's provided epsg is then the projection epsg code utm zone 32 north on wgs84 the spatial resolution is 30 meters and the dm subset that will be created will have the name dmsubset.tiff and then we run the function with these variables as an input let's run the script and check if the output is correct it should appear in the browser panel sometimes you need to click the refresh button and there we see dmsubset.tiff and if we remove the mosaic we can clearly see the result of clipping a reprojection it's tilted because the projection of our project on the fly reprojection is in the geographic coordinate system epsg 4226 and i'm changing it now to utm 32 north and now we see that it's nicely clipped to the bounding box that was provided because we're going to use pc raster the next step is to convert the geotiff to the pc raster format and we use this script here which uses gdall translate i can copy the script from the tutorial from the jupyter notebook after running this script you will see a copy of the dmsubset which should be exactly the same but then in the pc raster format so far we have been using only gdall so now we have the dm in pc raster format and we can use that for further calculations first we need to read the dm to memory from the file but we see that we get an error and that's because we didn't import yet the pc raster package so the readmap operation is not known but after importing pc raster this works the next step is to calculate the flow direction after filling the sinks in pc raster this is done with one operation which is called ldd create and it has four parameters to control the filling of the sinks you can find that in the documentation here we will use very high values to make sure that all the sinks are filled so the input is the dm and then we use the high numbers here for all the parameters to make sure that all sinks are filled running this will take a while and you need to wait until you get back the python prompt here now it's back and you can visualize the results or continue the calculations we can also run the agila operation here from qdis and you will see what happens if we do agila flow direction it will create a pop-up with the map and the nice thing about pc raster is that the ldd data type is visualized as these lines so we can follow the flow direction and we can query the values here it says that pixel flows to the east and one to the southwest west north etc in qdis there are also tools to convert the ldd format to a grip file and to visualize it using the mesh styling functionality there's another video showing that sometimes you're also interested in the filled dm which is skipped in this step if you want the filled dm you need to use exactly the same function but then change ldd create to ldd create dm so i create a new variable dm filled and i change this to ldd create dm again this takes a while and then the prompt returns you can also easily calculate the difference between the dm and the filled dm and you see that in this way the python console is a nice raster calculator if you use the pc raster operators and you can visualize it quickly in aguila or report it to disk and open the file in qdis in the map canvas so here we see the areas that are filled and how much they are filled and these are obviously the depressions caused by the mines i'm going to save the flow direction to disk and i call it ldd.map ldd local drain direction that's the pc raster terminology for flow direction and then i can drag it to the map canvas and there is our ldd map you can see that by using the python console we lose our projections when we use pc raster operations when you use the pc raster tools plugin with the interface then we make sure that the projections are kept with the layers that you calculate the next step is to determine the rivers and we'll first use the strahler order method so to create the strahler orders we use the stream order operation input is the flow direction map let's write it to disk for visualization and we drag it to the map canvas this is an ordinal raster so we need to visualize it using the pelleted unique values renderer and we can choose a ramp from white to blue so the more blue the higher the strahler order and the bigger the stream and then we need to calibrate this with a reference layer to find out what the streams are when we find out how many strahler orders there are we can also create a loop to generate for each strahler order a raster so let's calculate the maximum strahler order by using the map maximum operation that will find the highest value in the strahler orders layer and when i visualize it using agila you'll see that this operation results in the same value for all pixels so the maximum value found in a raster is assigned to all the pixels in the raster that's a global function now with this information we can create a script that loops over the strahler orders so the maximum strahler order is this map maximum operation that we just used but we need to have a value because we cannot create a range with a raster so there is the cell value operation from pc raster which has as an input a raster layer and then the x and the y row column value and it will report the value at that position so in our case the entire raster has the same values so 0 comma 0 will do the job but the result is a tuple and we need the first element of the tuple so we also add a line of code maximum strahler order value equals maximum strahler order tuple and then element 0 and then we can create a loop for order in range from one because strahler orders always start with order one to maximum strahler order value but the range in python does not include the upper boundary therefore we need to do plus one so for each of these orders in this range do stream equals if then strahler orders larger equal than the current order then write boolean one for those pixels and no data for the other pixels and write this to disk and we create a string for the file names it's called stream and then we add the order number and dot map the first one for pixels with strahler order one or larger the second one for strahler order two and larger etc etc and we can use that then for the calibration so we can run the script and we will see in the browser panel that this results in stream with strahler order numbers and we can drag it to the map canvas we need to hide a few layers and we can add open street map in order to see how well this matches with the rivers on the map so under xyz tiles you can find open street map or use the one from the quick map services plugin and here you can compare the results with the rivers on the map and determine which strahler order threshold we need to use to determine or to calculate our rivers you can see that they better match upstream which is more natural than downstream another way of deriving the rivers is using the flow accumulation instead of the strahler orders so in that case you can use the aquiflux operation from pc ruster which needs as an input the flow direction layer and material layer and in our case we use a value of one which means that we accumulate the amount of pixels we can write this to disk and visualize it in the map canvas the values are quite extreme from very high to very low values so i changed to single them pseudo color and choose a way to display this better also here i use a ramp from white to blue and i'm going to change the bin max values and what works very well with these extreme values is a cumulative count cut where you see now the bigger streams in darker blue we also need to calibrate this result with the reference layer like open street map and then if you have determined the accumulation threshold above which you consider something a river pixel belonging to a river then you can use the if then operation and in this case we do it with flow accumulation larger than 9000 those pixels will get boolean one and we consider those as a river let's write it to disk and visualize it and we can also compare it with the straller 8 layer here we see the flow accumulation result and you see that it's a bit different from the straller order result but the larger streams are similar the next step is to define the outlet on the stream and we can find the outlet by tracing the stream downstream and see where in this case the rule gets into the most and we copy the coordinates from qjs we can open a terminal from qjs which ends up in the folder where we store our data and there we can create a text file you can also use notepad but here we use copy con creation dot text i paste the coordinates and i add an id value one so the pixel that we will create will have value one and use control z to close this editor to save the file and i can use the call to map tool i do dash n because i want a nominal output location dot text is the input table location dot map is the output roster in pc roster format and the clone which determines the extent and pixel size can be any of our maps because they all are the same dimension and here we see that it has been converted so we can check it in uh with map attributes and there we see that the minimum and maximum value is one and we can also visualize this in qjs to see if it ends up at the location that we determined we need to style it palleted unique values because it's a nominal data and we see there the red pixel which is our outlet so we will use that for the lineation of catchment so create a new variable outlet by reading the map from the disk and i can use that in the catchment operation to delineate the catchment i call the variable rural catchment because that's the name of this catchment use the catchment operation which needs the flow direction as an argument and the outlet pixel and i write the result to disk so we can visualize it in qjs when i zoom to the layer i can see here the entire delineated catchment sometimes you don't have an outlet coordinate but you want to delineate all the catchments in the dem then you can use the pit operation to find all the depressions that it will consider as an outlet and then you use the output of that to use the catchment operation to delineate the catchments belonging to those pits let's write it to disk and visualizing qjs you can use palleted unique values because they are nominal and here we see all the catchments and we can also clearly see the rural catchment in there there are also ways to delineate subcatchments and that's covered in other videos and tutorials there are also tools available additional tools available in the resource sharing repository of qji now we would also like to automate this whole procedure instead of doing all the steps with separate commands the command line so here we have all the different steps written in functions and we need to remove the change directory because we assume that we are already running this from the directory where the files are here we define the different inputs we need to adjust the coordinates to the ones that we have used you can simply use type and then in locational text to get those values I can copy them save the script make sure that you have removed all the outputs from the previous exercise so you start clean make sure then that you set the directory and then you can run the script and when it's done it produces the result with akila but you can also find the files on your hard disk and load it in the map canvas with catchment and channels so you see that by automating these things you can have workflows calculated really fast you can just play with some input variables to change the coordinates or the data layers and this is only a small step towards creating a tool in qji with an interface which is basically a python script with some additional elements belonging to the template of the framework of the processing toolbox in qji s