 In this video I'm going to demonstrate the PCRaster map algebra tutorial. The tutorial is in Jupyter notebook but in this demonstration I'm going to show you how to run the commands in a setup where you have the Anaconda distribution on your computer and you use a spider and you can use the tools from PCRaster from the command port. So first we're going to activate the environment, Conda activate tutorials, now we are in the correct environment, I go to the folder where the tutorials are, go to the map algebra folder there we find the Jupyter notebooks and we see the data folder. So I'm going to run spider but I still want to use this command prompt. So I do start flash B spider then it will start spider and I can still use the command prompt. So we have these different tip files and we need to convert it to the PCRaster format to proceed with this tutorial. So we need to do some preparations and we're going to use Gdoll, previously you have used Gdoll from the command prompt and now you're going to see the equivalent in Python. So I'm going to use this iPython prompt here in spider and import the library, you need Gdoll and you need Gdoll const, then I'm going to open the tip file and I can set here the working folder, so I'm going to change that so I don't have to change that all the time or so I don't have to explicitly mention it every time. So we are going to import the first tip file, the build g.tiff and we call the variable src underscore ds because Gdoll open, that's where we start with to open a data set, we can use single or double quotes to define the string of the file name, build g.tiff and we enter and now the tip file is opened with Gdoll and then we can translate it to the PCRaster map file, on the command prompt we would use Gdoll translate and remember that then we had to specify the output format and the value scale of PCRaster, here I will show the equivalent in Python, we define the output variable equals Gdoll.translate, that's the same as Gdoll translate, but a little different syntax, first we define the output data set name, that is build g.map, then the source data set, that is src, ds that we have just opened, the gdoll tiff, then we can give several options and what we need to give is the format, so that will be the same as minus of on the command line and the format is PCRaster, the output format and we need to give the output data type and the comma here output type and that would be Gdoll const dot and then you give GdT underscore int 32, so that's the way of specifying the integer 32 output type which is needed because this is a nominal map and what we need to also specify in the PCRaster data type and we do that with the metadata options and there we use value scale nominal, so the output will be a PCRaster nominal map and we close the bracket, so that's done, the file is being converted but not written to disk yet, therefore we need to flush the data set to disk by closing these variables and you do it simply by dst, ds equals none, then it will be written to disk and it's also nice to close the source data set in the same way. There it is now, let's check our disk, so I'm in the data folder, type there and I see there that we have now build g.map which is the PCRaster format and I can even visualize it using Agila, that's the PCRaster tool for visualization and there it is, it's nominal, we see the different classes, five classes and this is the map, so that all works, so we can now write a script that can translate all the maps at once based on these lines of code that we have used, I'm going to paste it here, we miss the GDAL const, I'm going to add that here and then this function needs as an input the source file name, the destination file name, the output data type and the value scale for PCRaster, then you recognize the commands that we've used and here let me put this on the next line by using backslash, so here we have the destination file name, the source, the format, the output type and the metadata options and we flush it and then here we can convert it with the different formats, you can see it completely and we don't need the first one because we already did that one, I'm going to hash it out and then we can save the script and I'll put it in a map algebra, so all the code is then relative to the map algebra folder if we run it and I'm going to call this one tiff to PCRaster and then I can run it, now it has run and let's check the result and now we see all the maps are here, you can use Aguila to view them all, Aguila, we have build g.map, we had dtm.map, gw-level.map and the roads, gw-map and now we can open them all and here we see the elevation in the area, land use, buildings, groundwater level and two different road types, you can close all the windows at once by doing file exit, so let's inspect the data before we proceed because for map algebra the properties of all raster layers used in calculations need to be the same, let's find out if they have the same number of rows, columns, coordinates and pixel size, let's first open one of these PCRaster files with gdoll, call it raster layer equals and then we learned gdoll open and then note that we are back in map algebra, so I'm going to change it to map algebra data, so that's our working directory, then I don't have to specify the folder and I'm going to just open dtm.map that we just created, now it's open and to get the numbers of rows and columns we can use different gdoll functions the raster x size and raster y size and we can use raster layer dot get description to get the relative path of the layer, so let's start with getting the path raster layer dot get description, maybe typo raster layer and then the number of columns is raster layer dot raster x size and number of rows equals raster layer dot raster y size and now let's print that and you can use this formatting option to substitute what's in brackets later with the variables as that's much easier than substituting inside the string and then you use format and then the first one is the path, the second one is number of columns and the third one number of rows, we need to close with two brackets and there we see it, there is dtm.map has 60 columns and 60 rows, to know the coordinates and pixel size we can use getGeotransform from gdoll, it returns a tuple with the following information, index 0 of the tuple gives the top left x coordinate, index 1 the west east pixel resolution, index 2 rotation 0 for north up, index 3 top left y coordinate, index 4 rotation 0 for north up and index 5 north south pixel resolution and then that should be negative, so let's use this raster layer dot getGeotransform and here we see the tuple with the information and let's make that also more readable, so again we do print and then we can say origin equals and then we fill it in later with the correct values which are rasterLayer dot getGeotransform index 0 of the tuple rasterLayer dot getGeotransform position 3 of the tuple and close the brackets, so now we have a readable origin coordinate and we can do the same for the pixel size print pixel size equals and substitute this later with format and rasterLayer dot getGeotransform position 1 and rasterLayer position 5 and 2 brackets to close it, made a typo, it's a dot here and here you see that the pixel size is 50, so the pixel size is in the x and the y direction and the y direction has a negative 1 because they are equal here it means these are square pixels of 50 units and the units depend on the projection, we can also use geot to get the projection information, the pcraster format doesn't store the projection information but when we converted the rasters it saves an xml file with the projection info in the well-known text format, the gdel function to use for that is getProjection rasterLayer dot getProjection and this will print the well-known text format for this projection, that's not very readable, so we are going to parse this using another library, the pycrs library that we already installed when we installed anaconda, so I'm going to import this library, import pycrs and then I'm going to read those projection properties in a variable rasterLayerProjection equals rasterLayer dot getProjection, so now the well-known text is stored in this variable and then I can parse it, so the coordinate reference system equals pycrs dot parse and then we need to define the format from ogcwkt, so it knows that it uses the ogc well-known text format and that is read from rasterLayerProjection which has the whole format, so that's done and stored in CRS and now we can check if it is a projected coordinate system, so we can then use isInstanceCRS, that's the variable that we made and then pycrs dot porch cs and it says true which means that it is a projected coordinate system and we are not dealing with a geographical coordinate system, but it will be more helpful to have the name of the projection, so for that we can use projectionName as a variable equals CRS dot name and we print it print projectionName and it says it is WGS84 and UTMZone30North, we can also get the units of the projection, so let's make a variable projectionUnits equals CRS dot unit dot unitName dot ogc underscore wkt, so from the projection use the unit, use the name of the unit and it's in the well-known text and then we can print it print projectionUnits and it says meters, finally it would be useful to get some statistics from the rasterLayer, let's calculate the minimum and maximum value because rasterLayers can have multiple bands, we need to select the band, in our case we have a single band rasterLayer, so we have to choose band1, it starts numbering from 1, not from 0, so let's first select the band, rasterLayerBand equals rasterLayer dot getRasterBand and then as we said we have one layer, one band, so it's always 1 and then we can calculate the minimum and maximum value of the band using getMinimum and getMaximum, so rasterLayer Minimum equals rasterLayerBand dot getMinimum and for the maximum, let's change it a bit and let's print them Minimum, rasterLayer Minimum 201, that's the minimum value and let's do the same for maximum, rasterLayerMaximum 298, so that's all very useful but if we want to compare multiple rasters it's better to write this in a script, I'm going to make a new script and this is the script that we need for that which contains all these code that we've already tested on one and it just needs the rasterLayer from GDAL that we have opened, so in the main script here you see the DTM layer, we open it and then we call the function with the DTM layer and it will run all these lines to get the information and print them in a proper way, so let's save this one and let's call this metadata.py and let's run it, there it is and it prints here in the console the results, so DTM 60 by 60, that one is for one band, these are rows and columns, the projection, units, the origin, pixel size and the minimum and maximum value, so we've made a script that can easily read metadata from the file and we can conclude here that all the projections and spatial properties are the same and that this is suitable for doing map algebra, so we can continue. So the first condition to find the inaccessible wells is that the wells need to be within 150 meters from houses or roads, so we need to create a Boolean layer with true for houses and false for other buildings and we need to create zones of 150 meters around the houses and then we repeat the same for the roads, so let's start with creating a Boolean layer with true for houses and false for other buildings and we're going to import PCRuster now, so from PCRuster import star, PCRuster is imported now and our data is stored in the data folder so I need to change it back here, when we have run this script it changed back to where the script was run so always keep an eye on where you are and then we can read with PCRuster the map file that we created in previous steps to a variable, so buildings equals and we use the readmap operator and here we read build g.map from the disk and it's stored in the variable buildings and we can easily visualize this with Agila that comes with PCRuster so in the notebook that is not possible because it is a separate program but here it is possible so we can run Agila from the command prompt there is an Agila Python operator here which will give a pop-up with the map but if we have multiple maps that's not so useful because it gives a temporary file so therefore it is recommended to not use this Agila from Python only if you want to quickly check the result but use it from the command prompt so I can do here Agila build g.map and that refers then to the file but if it's for the variables in the script we use Agila and then the name of the variable and it will give a pop-up so that doesn't work in the Jupyter notebook so let's have a look at it again and we see it has five classes class 0 is none class 1 is house class 2 is public building class 3 is landfill class 4 is industry and class 5 is mine that's given in the tutorial and we first need to create a boolean layer with true for houses and false for other buildings with not algebra that is really easy you can say houses equals buildings equals exactly one so one equal sign means added to the variable and two equal signs means it needs to be exactly equal so they have a different meaning so we have now created the houses variable which is a map and which we can visualize so let's do that and there it is so in green true we have all the houses and in red false we have all the other classes that's good that's what we need to proceed with and then we want to create zones of 150 meters around the houses and therefore we use the pc raster spread max zone function so remember you can find all the pc raster functions on the pc raster documentation website so let's apply this houses 150 meter equals and then the name of the operator spread max zone and then the input is the boolean map with the houses then we have to define an initial friction which is zero there's no initial friction but each cell that it crosses has a dimension of one friction which is the distance and then the maximum distance is 150 meters close the bracket i can run it you see that's quite quick so i can visualize it houses 150 and this is the result a buffer around the houses of 150 meters so we can do the same around the roads so let's have a look at the roads map let's first read the roads map and then have a look at it there are this and there are two types of roads and we want to have the distance to any type of road so one and two are valid so i'm going to make a new layer um is road equals roads and road should not be zero so if it's zero it is false and if it's one and two it's true that works we can check it if we can now see all the roads there it is we have all the roads and now i need the maximum distance of 150 meters so with the arrow up i can get back the one that we used for houses and we replace it with roads 150m and here we use this road to press enter and then we can visualize it and there we see the 150 meter buffer around the roads now we can save these files to disk using report so report and then the name of the variable so houses 150 meter and we save it as houses 150 meter dot map and we do the same for roads there it is and we can then check it and their houses 150 meter and roads 150 meter are added so we can with Aguila visualize them there they are the second condition is that we don't want industry mine or landfill within 300 meters from the wells for the second condition we first need to reclassify the buildings layer in such a way that the result is a boolean map with true for industry mine and landfill and false for the other classes then we need to calculate the distance to the true pixels and finally calculate the pixels that are further than 300 meter from industry mine and landfill so let's create a boolean raster with true for industry mine and landfill and false for other buildings for that purpose we can use so-called lookup tables and we use the lookup operator and it depends on the output data type how we write that operator so if it's a boolean we use lookup boolean it's nominal lookup nominal etc and we need the table the lookup table let's have a look at the lookup table so the table that we are going to use here is industry dot tbl and here we see it and it has on the first column the values of the raster and in the second column the recoded values and here we see that classes 0 1 and 2 will be 0 so boolean false and 3 4 and 5 will be boolean 1 so let's use that so we can say industry equals lookup because the output is boolean it will be lookup boolean and then we use the table industry dot tbl and the buildings layer that we have previously imported there it is so let's visualize the result so there it is in green are all the industry areas and in red are all the other areas so now we need to create zones of 300 meters around industry mine and landfill so around these green pixels and previously we used the spread max zone operation but there doesn't exist a spread min zone to get the minimum distance so there are different ways to solve this and i'm going to show that so the first method is a fastest method i think is we can use spread max zone and then invert the result by using the not operator so industry max 300 meter that's our output variable spread max zone as an input we use our boolean industry map initial friction zero then one cell for each that we cross and then a maximum of 300 there it is and then if we want the minimum we can say industry min 300 meter is the opposite is not there we use this symbol industry max 300 meters so not the tilde inverts to and false let's have a look at both maps and there they are indeed each other's inverse we can also do it in a different way we can calculate all the distances to industry and let's call that industry distance equals and there we use the spread function for pc raster and the input is the boolean industry map initial friction zero and then one for each cell that it crosses there it is let's have a look so here for each cell we have the distance in meters to an industry pixel note that you can change the scale by clicking right and then add a drop properties and you can choose another palette closer to industry would be bad so let's make it red to green now with this result we can create a boolean condition industry min 300 let's make it version 2 so we can later compare if they are exactly the same industry distance that's the distance map larger equal 300 so this means if industry distance has pixels which are larger equal to 300 meters then it's true else it's false you can visualize the result and there we see that it looks the same and in the third method we can use a lookup table and in lookup tables we can also use ranges so i'm going to show you that file this is the file and this means that any value less than 300 has zero and values larger than or equal to 300 get boolean true so let's use this industry min 300 equals lookup boolean the first argument is the table the map which is industry distance let's look at this one too it looks also the same now to be sure we can plot all of them and we see that they are all exactly the same so let's save this to the disk we can just choose one of them report industry min 300 m i'm going to save it as ind 300 m dot map and now it's on our disc for the last condition we need to identify the wells that are less than 40 meters deep the file groundwater level dot map gives the absolute elevation of the groundwater level in the well in meters above sea level in order to calculate the depth to the groundwater we need to subtract this from the surface elevation given in the digital terrain model dtm dot map so let's read the two files and with map algebra it's very easy to subtract them so the well depth equals the surface level from the dtm minus the absolute groundwater level let's look at the result there we see it continuous scalar map which gives the distances to the water the depth now we can calculate a boolean raster where pixels with a depth less than 40 meters are true and larger equal to 40 meters are false so let's call it not deep equals well depth less than 40 so that will give two for the wells that are less than 40 meters and false for larger equal to 40 meters let's look at the result and there it is so let's save this to the disc report not deep not deep dot map so we have calculated different boolean maps for all these criteria so let's combine them and the result should be a boolean map where the accessible wells are located at the cells where all the conditions are true therefore we need to use the ant operator so we can write the equation here accessible wells equals and then all the different boolean maps with the ampersand symbol which means end for pc raster and then we can visualize the results with only these three wells being true so accessible and the other ones being not accessible so what's the advantage of scripting map algebra because the same procedure can also be done in qgis as you can see in another tutorial so what's the advantage of using this instead of a graphical user interface well first you have more control over the functions because you determine each argument and operator that you use the process becomes reproducible anyone else can run the script with the input data to obtain exactly the same results the code is transparent so open source and you can easily modify the criteria and see how it affects the result for that last purpose we're going to put all this code together into one script and we make it also clear by writing certain processes functions and then we can play around with these different criteria and see what works best so it's a bit of a scenario analysis i'm going to make a new file let's run it you see it's fast this is the result but we can now play around with the results we can say what if we increase the first condition to 250 meters save it and run it and we see that there are two more accessible wells here in the north so with this script we can play around with these different conditions and see the effect of that which is very useful for planning