 Hello, this is Hans van der Klaas senior lecturer at ICH Doft Institute for Water Education. Many people have asked how to delineate multiple subcatchments automatically using PCRuster Python. In this video I'm going to show you how to do that. In a previous video and exercise I've showed how to delineate catchment. This video builds on that. So let's first import PCRuster. Then we need to read the DEM file into a variable DEM. So I use the readMap function from PCRuster and I read dm.map from disk. I'm going to add some print statements so the user can see where the script is and first we need to calculate the flow direction. And we can calculate the flow direction with the lddcreate function where the input is at DEM and then we have a few arguments that control the filling of the DEM which we all put on a very large value. Let's report that to disk because later we might want to just use that file. It will take some time for it to calculate for this large DEM that we use. It's the root catchment area with four SRTM tiles. So the next step is that we want to calculate the strahler orders and strahler orders can be calculated with the stream order function from PCRuster and it uses the flow direction. Not all strahler orders are rivers. So we need to use a threshold and only calculate the strahler orders larger or equal to the threshold. The user FIFDEN function here, if strahler orders larger equal than a threshold, which we're going to define, then strahler order and here at the inputs I'm going to add threshold equals, let's use 8 here, the higher number the larger your subcatchments will be. Now we need to find the outlets and the outlets are at those junctions where one strahler order gets into another one in the downstream direction. So we can use an FIFDEN function and then we use the downstream function which looks at the value downstream of the flow direction. So that's an input of a certain layer to use that as strahler. If that is unequal to strahler rivers, so if it's unequal to the current pixel, then we have a junction and then we give it boolean 1, which means true and then the outlets are those points but they need to have unique values and we use the unique ID function but later we use it in the catchment function so it looks at all values that are larger than zero, the non-zero pixels, so we need to fill the rest with zeros and then we use ordinal because we want to convert it from scalar to ordinal in this case and we can't use nominal because the map maximum function that we use later doesn't accept that. Now we can calculate the subcatchments map with the catchment function from PC raster input flow direction and the outlets, so for each unique value of an outlet it will calculate the catchment area. So let's write that also to disk with the report function, call it subcatchments.map and let's see it after running by using argila, save it, I can run it and there we see it but this is not exactly what we need because the more downstream junctions or outlets they have subcatchments that cover the ones that are more upstream because this is a nested system of course. So what we would like to do in fact is to split them out into different files so for each outlet a separate subcatchment raster and that's what we're going to produce next. So I'm going to remove or comment out this part because we are just going to read the flow direction map, takes too much time to calculate every time especially when we are debugging, we'd map ldd.map and I'm also going to comment out this part of the subcatchments, we're going to do that differently, we can also use shortcuts for that. So now we are going to first determine the maximum amount of outlets, so how many outlets are there because we need that value for looping over the different outlets. And we can use the map maximum function from PC Raster to determine that, the problem is that it will generate a PC Raster map and we need the value. So we can calculate the value, but the result will be a tuple with the cell value function. So use your cell value, as an input I use maximum outlets and that will give me back a tuple and I need to give it also the coordinates and because the map maximum is homogeneous for the whole map I can simply use 0 comma 0, the corner coordinate and now element 0 of the tuple gives me the value. So I define maximum outlets value as maximum outlets tuple, element 0. And let's see if that's indeed the case, we're printing it to the screen, save and run it, it goes much faster and we see that there are 56 subcatchments here and now we can make our for loop for each outlet in the range from 1 to maximum outlets value, but we need to do plus 1 because the range function includes the lower value but excludes the high value, so we add one to have them all and we want to give the user some feedback where we are, so processing subcatchments and then the number which is outlet and it would be nice if that would be not on a separate line all the time so what we can do is use the return carriage, the slash R and we have to add at the end some more code to make that happen and here you see what's needed and then we can calculate the subcatchment of that specific outlet, so subcatchment equals catchment function, so very similar as we did before but now with only one point, so we define the flow direction there and then we need to get a map with that outlet number, so we use if then else, if outlets that map equals the value of outlet then give me boolean 1 and else boolean 0, in this case it will generate a catchment with boolean 1 for inside the catchment and boolean 0 for outside the catchment, I also want to write this to disk, so I use report, subcatchment variable where the map is stored and then I need to build a string of the file name and let's go with subcatchment but add there the string of outlet, so make the outlet number a string and add it to the file name and add there dot map, it should have been dot map and visualize it and then we print to the screen and it's done, save it and run, I made it a bit faster but here you see all the subcatchments popping up, hope you've enjoyed the video, more free materials are available at GIS OpenCourseWare