 Hello, in this video you'll learn how to do spatial interpolation of Borgl data with PCraster Python. This video follows the tutorial on GISOpenCourseWare.org. You can find the link in the description of this video. For this video I assume that you have installed Miniconda and created the environment and downloaded the tutorial as explained in the previous video. A nice trick is that you can pin the terminal to the taskbar. You can pin anything from the start menu to the taskbar by clicking right. So you can launch it quicker and of course in this course we will use the prompt a lot. So let's first activate the tutorials environment that we've created before and change the directory to the PCraster tutorials in my case on the Z drive. And then we launch JupyterLab by typing JupyterLab. Now JupyterLab will start in your browser. And let's set it up for our project here on groundwater interpolation. So I go to the subfolder, groundwater interpolation, where I find the data and the tutorial. The tutorial is in the Jupyter notebook. If I double click it will open in a tab. It's important to link the Jupyter notebook to a console. So click right and choose new console for notebook. This will add a tab with the console, a Python console. And it's also useful to have a terminal. You can drag these tabs to the position where you want them. So we can have several items here on the screen that we will use, the notebook, the Python console and the terminal where we can use the command line interface commands that you learned in module one. So now I'm set up to run a code from the Jupyter notebook. But I can run the same code also from the console. And because they are linked, the line numbers are following each other. So if I run first something in the Jupyter notebook and then on the console, the line numbers will follow up each other because it uses exactly the same kernel. So I've imported PCRuster and OS package. And first we need to convert the CSV file to the PCRuster map format. And in fact, that is a command line command called to map. If I type call to map, I get an explanation on how to use it. We can also find this in the documentation of PCRuster. And if I list the directory, I see I need to go to the data directory. And there I can have a look at the CSV file by using the type command. And here I see that the fourth column is the groundwater depth. And then we have an X and a Y coordinate in a projection. If you receive a CSV file from somebody, then the projection also need to be given to you because this is not stored in the file. Here we are not caring about that. So I can use call to map and the data scaler. So use dash S, the X column is column number two, the Y column is number three, and the value column is number one. Then the CSV file is the input. And I give an output file name, which is a PCRuster map. I call it borehole depth dot map. And it needs to use the properties of our clone map. Clone map is a mask, which defines the extent, pixel size, and coordinate system. So now we've imported the boreholes. It says it read 220 records. And they're all inside, but there are some cells with more than one record. In the Jupyter notebook, you can see how to execute a command line command from Python using os.system and then the command string. Now let's check the result. There's the Aguila command that we can use to visualize PCRuster maps. Note that the map will pop up behind your web browser. So you need to go to the taskbar and to find the map. Now for the Python part, you can see that we are not yet in the data folder. So I need to change the directory. And these are a few Python commands that you can use to run command line commands. So change the directory to data. You can then use again this print os getcwd, current working directory. And I see that I'm now with Python also in the data folder. Because the Jupyter notebook is linked, also the Jupyter notebook is now in that folder. So now I can read the map that we created from call to map. That's with borehole depth. If I use tap, I get a list of files that correspond to the pattern that I was typing. So that makes it faster to select the right file name and avoid typos. And then I can use the Agila command to visualize this borehole depth map. And note that also here it's behind the browser. And so you need to be aware that it doesn't pop up in front. So now we can interpolate using the inverse distance interpolation. There's a PCRuster operation called inverse distance. And it needs a mask, which is the extent that will be interpolated. So we need the scalar raster layer with the points. And we need to have the power, normally we use a power of two, a radius and a maximum number of points that will be used in interpolation. We'll set that to zero. So first we read the clone map that we will use as a mask. And then we can execute this line. So we can also execute this from the Jupyter notebook. And you see here that we are at line number nine. So number eight was typed in the Python console. And number nine was executed. So they are linked and follow each other. So if you define variables in the notebook or in the console, then they are recognized. So let's visualize. Note that there's a difference between using Agila from the Python console or from the terminal. From the Python console, you visualize variables that are in memory. And also the title bar of the map is a bit cryptic. From the terminal, you can visualize a .map file that is saved on your disk. So to show that, I can report this variable to a .map file. So with the report function, we can write maps to disk. Warhol depth IDW is the variable. And then in quotes here, I have the string of the file name warhol depth IDW.map. And then I can simply run Agila from the terminal. And also there it pops up in the back. And I will see here the result. But now the title bar shows the file name, which is much clearer in practice. So use Agila from the Python console just for quick checks and for final results. Save it to the disk and use Agila from the command prompt. So IDW is one way of interpolating points. Now we'll have a look at using tcent polygons. And the steps are assign a unique ID to each warhol, then assign to each pixel the closest warhol ID and then the corresponding values of the depth. So I run this defined operation to get a boolean layer which is true at the warhol depth and false at the others. And this is then the result. You see I now mix the notebook executable lines with the console ones, which is perfectly okay because they're linked. So now we have a boolean with true for the warhols and false for the rest. We can give each warhol a unique ID and that will be a scalar result and therefore we also convert it to a nominal. And then this is the result, a nominal layer with a unique value for each warhol. If you zoom in you can find there the warhols and you can query the values. Then we need to use the spread zone operation to create these tcent polygons with the unique numbers that we have assigned. So it needs the warhol ID, initial friction of zero and friction distance for each cell at one, which means basically that the distance is one pixel, that's the unit. And this is then the result, very nice polygon map. But now we need to assign the values of the warhol depth to each tcent polygons and we can use that with area minimum or area maximum because in the area there should be only one warhol and if I execute it and visualize that then this will be the result. If we compare IDW with tcent polygons, very different results so it really depends on the amount of points that you have and the assumptions for the interpolation. I can write these results to disk using report. The variable warhol depth tcent will be written to a string with the file name, warhol depth tcent dot map, shift enter to execute and then in the terminal I can visualize this. You can also use the prompt in the same way, it's exactly the same and here we visualize both layers and in the title bar you can see which is which, therefore it's more useful to use the terminal for that and we can compare the results because these viewers are linked. You can use file exit to exit all the linked windows. Now for IDW there was a nice operation, it would be really nice to have a Python function that does the same for tcent in an easy way. We can write that with the knowledge that we have. So you can write it here in the notebook but we can also create a Python file. If we choose here Python file, it opens an editor and if we drag it there to the right then we can have on the left the notebook and next to it we can write the script based on what we learned. So first we need to import PCRuster and then we define a function called tcent and it needs one argument, the borehole depth and first we need to get a Boolean layer which gives true for the boreholes and false for the rest of the roster. Then we need to give unique IDs using unique ID but that results in a scalar and we need nominal so we also convert to nominal and then we can get the tcent IDs so that's the polygon map but then with the unique IDs with the spread zone function with the cell size as the unit so therefore the last argument is a one and then the result borehole depth tcent is the area maximum of the borehole depth and the tcent ID so for each polygon it will find the maximum value in the polygon and this will be returned after executing the function. So we have created a function to calculate tcent polygons from pixels and here we can run it so we read the points from the map so the borehole depth PCRuster map and then we create a new variable depth tcent and we use our tcent function and it needs as an input the borehole depth so that's depth and then we want to see the result so we can use Aguila and we can also save the result to disk using the report. So this is our script we can save it to a Python file so we use savePythonFileS and save it under the data folder so we don't have to deal with the paths because our file names here are defined in the current folder which is the data folder and I call it tcent.py. I'm going to save it I can find it there. I can create a console for this editor so now it's linked and then I can run the code I can run line by line or I can run all the code and here you find the result and here you can also see that depth tcent.map was created and you can run this Python script from any terminal that has access to Python so also your miniconda prompt and then you just type python tcent.py and it will run. You see that the prompt now hangs that is because the Aguila map is in the background and it will give the prompt only back if you have closed Aguila because then I will go to the next line of the script which is nothing so it closes the script.