Create a groundwater level map from borehole data and a DEM

okay so um for today uh the first tutorial is a bit more about uh processing the data and uh you will see some more things there also be some uh repetition of what we did before so um go to the to the gis training on the open courseware platform and there um there's this uh first uh tutorial of today about creating a groundwater level map from borehole data and the dem and you can download the data before doing the tutorial and it is a polygon with a boundary it's an imaginary boundary because there are not many boundaries in the area that i wanted to use uh you could use the boundary of a municipality or of an aquifer here i just took a random extent like you learned last time how to make a polygon of the extent i'm going to download this folder because it's a shape file and you know from good practice that the shape file is not one file so i'm going to download the zip file and save it to a new folder and i'm saving the zip file there and then i go to that folder tutorial 3 and there's the zip file and i have 7zip that's the tool that i recommend it's open source for zipping and unzipping and if i click right i can choose here 7-zip and i can say extract here what's always good practice is first to open the zip file with 7-zip open archive and then you can see if there's folders in it or files if there's a folder in it then you will extract it at the folder level in this case there are files so we'll just extract the files in the subfolder which will be good enough so i'm going to extract all these files here to tutorial three and there they are i can remove the zip file not to get any confusion and there is and then i can start with the tutorial and we will load again a data set from a geonode we will calculate the density of boreholes within that study area we will download the srtm one arc second digital elevation model click and reproject it and then we style the dem and we sample the elevation from the borehole locations and we compare it with the elevation attribute which was registered at the boreholes we make corrections to the attribute table and some calculations and then we're going to end with an interpolation of the groundwater levels at the boreholes to a raster file so i'm going to start the tutorial um it starts with getting the data set from a geonode it's uh we're going to work on the stump reit transboundary aquifer and we get that data on the boreholes from the orange senku river basin gis server which is a geonode and here it is and we're going to use data from the area which is called stamp read which is a place and there's some data about stamp rate it's a trans boundary so have a look at it so the geo node you can always look at the metadata you see it's in botswana namibia south africa it's a real trans boundary and we're going to focus on the stump reit which is uh approximately over there to look at the the hydrological hydrogeological situation there and you can check the attributes and it comes with quite some attributes that are registered with this layer in the geo node as we learned we can make a connection with qgis so we're not going to use the download option and get a zipped shape file but we're going to make the connection so i'm gonna move to uh qgis and follow the steps of the the tutorial so i'm gonna make uh well first of all again the tip you can put your layers panel over the browser panels we have the two tabs and i'm going to start with making a new geonote connection to this server click this button and go to geonode make a new connection and this site is called erasicon and the url you can also copy it so here put the url don't make typos better to copy it and i'm gonna test the connection and again we are depending on the internet and it says the connection was successful it's a valid geonote instance so that's great so i click ok and i'm going to connect to load the layers it will give an error but this error doesn't affect what we do so if you click ok you will still get all those layers and what i'm going to do now is look for this stump read data so i'm going to filter on stump read and there you see it and remember that there are different web services and if we really want to work with the vector data we need a wfs web service and the one that we're going to use was the one that we were looking at that's this one and then the wfs so i'm going to click add and then it will take a bit because it's a very big data set so um let it load i'll do close and we'll just wait for it to uh to load many points but we will make it more useful by cutting it later to a specific study area you see the progress and there it's coming many many points great that these data sets are available on sdis so you can use them so we're going to the next step and that is that we want to export this to a local file so you remember maybe how to do that so click right and i do export save features as and there i'm going to save this to a new geo package that i'm going to make for this project i'm going to give it a new name that i call stamp read data that's the name of our deal package where we're going to store all the layers and the name of this layer we'll just call this one boreholes and i also want to change the projection to the projection that we use in this area which is utm zone 34 south on wgs 84 so i click this button to set the projection and i use the eps g code so 3 2 is for utm 7 is for south and then you get the zone number 34 and then don't forget to select it here so you can see that it covers our study area i click ok it's filled in here and that's all to save this boreholes layer with the projection in our new geo package and to add it to the map so do okay that will take a little bit to convert it and there it is in a different color qgis just assigns random colors to this you can always use the styling panel to style this and we don't need the geonote layer the online layer anymore so i'm going to remove it and you see that if i hover my mouse over it that it will say that it's the layer in the dual package with the utm coordinate but you can also notice here in the lower right that on the fly projection of this project is in 4326 which means latitude longitude so i need to change this another way of changing it is to go here and go to set crs and then we can say set project crs from layer if i click this the project will have the same projection as this layer so that's a nice way to change the projection of the project and let's save the project at this point i'll just save it as a normal project let's call it stumpreit there it is and let's continue the next step is yeah this is a huge data set the attribute table is also very large and will take a lot of time to display so i'm going to make it smaller and i'm gonna cut it to the boundary that was provided so i go here to the browser panel and i go to our tutorial folder tutorial 3. i click right and i make it a favorite so i don't have to browse to this whole list and here we see the geo package that we just created with the boreholes and here we have the shape file that we downloaded from the tutorial portal and i simply drag it there i can check the projection it has the same projection so that is okay so it means we can do all kinds of overlays with the borehole data i click right i choose zoom to layer and yeah it covers the points so let's style it a little bit go to the styling panel and i change the simple fill to a simple line make it red and make it a bit larger so this will be our imaginary study area so remember that you could choose there anything that is the boundary of your study area it can be the boundary of an administrative unit a national park or an aquifer or geological area that's up to you but here i just have this imaginary boundary and now we're going to clip the boreholes to this boundary and save it to our geo package so we've done these things before you go to vector geoprocessing tools and then to clip and there i choose as an input layer the boreholes and as an overlay layer the stumproot boundary and here you can check again in the brackets if they have the same projection that is very important and then i'm going to save it to the geo package so always use this to browse and then save to geo package and i'm going to choose here tutorial 3 stamp lead data save it will ask us how do you want to call the layer and i'm going to call this boreholes clipped and click ok and now it's going to clip it close and i can remove this layer and now you see we end up with only the boreholes inside our red rectangle now the next step is that we're going to calculate the density of these boreholes in the study area and therefore we use a processing algorithm so i'm going to close the styling panel i go to processing to toolbox and there you can find under vector analysis count points in polygons we need to first well mathematics if we want to know the density it's the amount of points and it's the area that it's covered and then we divide those so let's start with counting the points in the polygon and here for the polygon we choose the stamp read boundary for the points we choose our pore holes that are clipped to the boundary and here you can change the name of the field of account because what it will do it will make a copy of your polygon so it will create a new layer but it will it will copy all the attributes but it will add an attribute with the number of points so number of boreholes and then let's save it to our geo package and because it's still the polygon and it has everything that we need let's call it the stamp read study area so it will be clear from now that we will use that layer as the study area boundary stamp with study area then i run it close and you see it's a copy of our rectangle we can copy the style paste it and then we don't need this one anymore and uh let's have a look at the attribute table so click right open attribute table and there we see that it's one feature that's this red polygon and the number of boreholes in the rectangle is 64. that's great but now we need the area of the polygon so what i'm going to do is i'm toggling on the editing and i'm going to add a field which has the calculation of the area and i want it in square kilometers so i'm going to open the field calculator and create a new field which i call area the units so the user doesn't get confused later when they see the field and change it to decimal number wheel and you can calculate geometries from your uh vectors then you go here to the geometry function and you can choose their uh area so for a line you can choose dollar length and it will calculate the length of a line here we choose dollar area which works for polygons and it will give us the area here we see a preview that's the area in the map units our map units are square meters so i want to change it to square kilometers i'm going to divide it by one million oops and then the preview shows us the area in square kilometers and do okay it adds the area field so we're almost there we have all the ingredients to calculate now the number of boreholes per square kilometer which is the density so again we do that in this calculator and the output field name will now be borehole density per square kilometer and that will be also a decimal number and there we are going to do the calculation with the fields and values and we are going to create equation number of boreholes divided by the area and then we see the preview looks okay we do okay so the borehole density is 0.2 boreholes per square kilometer so this is just a demonstration on how you do calculations in the attribute table with these expressions i'm going to toggle off the editing mode and i'm going to save it now in the next section we are going to add a digital elevation model because we need the surface elevation also for the areas where it's not reported in the boreholes data set and we need to compare those values close the attribute table and we're going to use the srtm dataset that's the shuttle radar topography mission the space shuttle uh mapped almost the whole world at the one arc second spatial resolution for elevation using radar technology radar interferometry and the spatial resolution approximately at the equator is 30 meters and is freely available and covers the whole static region and much more there are different ways to download it so you can download it from the usgs website if you go to earth explorer then you can upload here the polygon you can go to data sets and you will find it here under digital elevation srtm the one arc second and with this information button you can get more information about this data always important as you have learned to check metadata and this gives us all the data information on what this data is and how we could use that so that's important but we are going to use another way we are going to use a plugin to download it into qgis so go to plugins manage and install plugins and there i'm going to choose the srtm downloader that's a plugin to download the tiles of srtm because it's a global data set which is cut into different tiles and it will download the tiles of your map canvas or in a specific study area you can choose that so i'm going to install it click close and it has added this button so i click the srtm downloader icon and i'm going to set the canvas extend so it will use these boundary coordinates of the map canvas and uses it in latitude longitude so it automatically converts it and i'm going to define the output path where i want to save the file so i go to the tutorial 3 folder and i want it to be saved there and then i click download it will add the username and password that's also on the earthexplorer website it's the same username and password so you can click this link and make a new username and password i already have one so i'm going to give it here and you can choose to save the credentials so it will not ask it next time once you have submitted your credentials it will start downloading in this case it only has one tile because i chose a small area for this tutorial and then it says download completed that's great so here we have our srtm one arc second elevation model and let's put the study area on top of it and we see that it covers much more than our study area so you can guess what the the next step is so after downloading we are going to clip and reproject the dm to our study area to do that well first of all let's check the projection of the global data sets are often in latitude longitude and you see here that it's in epsg-4326 and that is latitude longitude well we are working in utm so we want to change that and we want to clip the dm to the boundary of the stamp read study area so i click right go to export and choose save as and i'm going to save it directly to our geo package so change this to geo package and then we browse to our geo package so i have all the data together which is good practice i'm going to change now the projection and now i can choose it from the list because it's the projection of this project and i want to clip it to the boundaries of the stomp root study area layer so choose that here and it fits them to the boundary coordinates of our red rectangle i'm going to change the spatial resolution to 30 meters make them square and there's one more thing because we reproject it's important to define the node data value i check the box and i click the plus and i'm going to add minus 999 because there's no elevation that is -999 that's an of range value and with reprojection uh it and clipping it might introduce some no data areas at the border and this will be assigned the value minus 9999 and the software that knows that it should not visualize those pixels that's all so i'm going to click ok and it's now added to the map you see a bit of contrast difference that's because the standard the default visualization of a layer is take the minimum and maximum value and scale the grayscale between those values and because we clipped it we have a different minimum and maximum value to scale the gray values i'm going to remove this layer so we end up with our clipped digital elevation model with 30 meter resolution now the next step is to style this one because it doesn't look nice in grayscale there are some nice tricks to uh to style it so i'm going to open the layer styling panel this is continuous data so for continuous rasters we always use the single band pseudo color and then we can choose a color ramp and i'm going to create a new color ramp and i'm going to choose here a catalog cpt city and that's a very useful catalog because you can find a lot of predefined ramps as you will see there it is and in all kinds of categories so i'm going here to topography and i'm going to choose elevation and i'm going to use this a ramp and if you want to reuse this all the time you can click save as standard greedy gradient and then it'll always be available in in the list now i need to click classify and now it has our pixels in the colors of the color ramp that's one thing we can do i'm going to duplicate this layer that's a nice trick to style dems so i'm going to uncheck this one and check this one because it's duplicated it will be exactly in the same style but i'm going to use another render let me for good practice rename this one see i didn't change the layer name i guess in the geopackage so you should call it dem i left it as default so correct that when you go through the steps and i'm going to create here a hillshade so it refers to the same geo package but just in the layers list i give it another name so it's we can recognize what it is and then here i can choose the hill shade and for the renderer there's a hillshade renderer that automatically calculates on the fly the hill shades and we can see the shading here you see it's quite blocky therefore we need to change this setting of resampling and when it's zoomed in we choose the linear and then you see that it's smoothed and we can also put this one on average when we're zoomed out and that improves a bit the visualization and when we blend this with the dm so i go back to the dm in the layer styling panel and i change the blending mode here to multiply and then we have the heel shades with the dm and that gives a better visualization effect now this area is not very spectacular in real leaf but if you try to do that in other areas you will see that that gives a very nice effect that was the styling of the dm then it's time to to look more at the data and what we're going to do is uh sample the dm values at the locations of the boreholes and we did that before in another tutorial we need a plugin for that that is called the point sampling tool here it is install the plugin click close and it adds this button so i'm going to add it and what i now need is i want a copy of this borehole's clips but then with the elevation values from the dm added to it so i'm going to select here all the fields from the original layer the boreholes clip layer i use shift to select them all and i want the dem so i keep control pressed and i click on dm so it's also blue so we can see here that every field that i want in the output is now selected in blue and i'm going to create an output layer so it will make a copy go to tutorial 3 i i'm not saving it in the geopackage because this tool will not add a layer but will overwrite it so i'm going to choose a shape file and i'm going to call this one or holes that so i will recognize that this is the one with z values and i click save before clicking ok i go to the field step and there you can rename the fields if needed and i'm going to do that for the dm value i'm simply going to call this one z underscore dem so i know that's elevation from the dem when i click ok now it processes the data set and yeah we can now work with the sample data set so that's the next step we are going to look at the attribute table of this new layer i'm going to remove this one and click right open the attribute table 64 records as we have counted before and we see different fields in the attribute table and for us what is important for this uh tutorial is that we have the elevation that is mapped at the boreholes but you see that there's a lot of missing data these minus 9999 are missing data values so some of these borals have missing records for the elevation at the borehole mapped before we have some other useful fields we have our elevation from the srtm dm here added which we just did with the point sampling tool and there is a depth in meters which also has a lot of missing data but it also has some with data that we'll later progress with so first i want to compare the elevation from srtm with the reported elevation at the wells to see if they have a bit of a relation with each other so i'm going to select first the no data values because i don't want them in a scatter plot that we are going to create so i click this button select features using an expression and i want all the elevation features i go here to fields and values all the elevation values that are less than zero i want them to be selected so here elevation in double quotes that refers to the field name and then less than zero and then i click select features and i click close and i see that there are ten features selected we see them here in bright yellow also on the map and those are the ones that have missing data for the elevation here's the elevation field but for our scatter plot that we're going to make with a plugin we need to invert the selection we want to make the scatter plot for everything except those missing values so i'm going to use this button to invert the selection so now 54 records are selected those are the other ones and i can close the attribute table for now they remain selected and what we're going to do now is to make a scatter plot and we're going to use the data plot leap plugin which is a very nice plugin to make scatter plots or other types of graphs from your attribute tables so go to plugins manage and install plugins and i choose here data plotly there it is and i'm going to install the plugin plugin is successfully installed i click close and i'm just going to remove the layer styling panel because when i click this new button here it will open the data plotly panel make it a bit bigger so we see the whole panel and here you can choose the type of plot that you want and i'm choosing here the scatter plot and i'm choosing the layer here boreholes set that's important and what we need to set is the correct x field and the y field but first i should not forget that i only want to use the selected features from the attribute table so i check this box and i want on the x axis the elevation and on the y i want the set from the dm that we have sampled and i keep the styling as a default you can play around with that and i go to this tab to configure some more things i don't need a legend i'm going to change the title to borehole versus dm elevation always important to define these labels on the axis with the units so good practice borehole elevation in meters and this one we call it dm elevation meters and that's basically all the settings that we need now and then i can click this button here to create the plot and there we see our scatter plot and we see that they are related there are some outliers but generally it gives us a good correlation so we might in this imaginary case conclude okay if we have those missing data in the borehole data set for elevation we can replace it with the elevation values of srtm as a substitute you could also go back to the field with an altimeter and then measure then the elevation again but here i'm going to show you how how this will work to replace those no data values so that's the next step i'm going to close this data plotting panel and i go back to the attribute table there it is and i'm going to deselect all i don't want this selection anymore because i'm going to work on this attribute table now and i toggle to editing mode and what you see here in the editing mode on the just above your field names is something that you can read as an equation so this is the output field equals and then this button to set the equation or you can type the equation here and i'm going to make an equation with a condition so i'm going to choose here elevation and what i want is a condition that if the elevation is no data so less than zero then give me the values from the dem otherwise use the values that were recorded at elevation so i go to this button to build the expression and here under conditions conditionals you can have different functions to make conditions we are going to use the if function and you can see the syntax here on the right if i'm making a bit bigger if the condition is true do something that is here and then a comma what to do if it is false so that's the structure if i double click i can start building the expression expression it is written here and first i'm going to add the elevation field so go to fields and values and double click on elevation so it will be here in double cl quotes make it a bit bigger so you can read it better so if the elevation is less than zero then that's a comma it reminds me here about the syntax when it's true i want the values from the dm so i go here to zdm and i double click and then i write a comma and i say what it has to do when it's false so when it's not less than zero i wanted to add the original elevation values back so i double click and then i closed this function with the brackets and i see that the preview is okay so to remind you it's a condition if elevation pixels are less than zero if that is true then write for those fields the elevation from the dm if that is not true use the elevation values that were already there so i click ok and now the equation is written here and then you verify if the equation is well so elevation equals if elevation etc and only when that is the case you click update all and now we will not find any no data in the elevation field so we have every elevation every feature every borehole now has an elevation feature with data so that's great i'm going to toggle it off and save these results in this attribute table now we need to calculate the groundwater level in the boreholes because we have here the depth in meters not for all the boreholes but for some of them but we want to have the absolute elevation so i'm first going to filter out those no data values so i go here to select features using an expression as we learned before fields and values depth in meters and i'm going to say if that is less than zero so all those -999 values then i want to select them so choose select features and it will select a lot of them we're going to lose a lot of these points in this case and i'm going to simply remove those so i go to toggle to editing mode use the trash bin icon and we are left over now with elf point 11 points that we are going to use for our further analysis so um i'm going to save this by toggling off the editing and click save and then i can close this and we see here those points now the next step is that we need to go back to the attribute table because i now want to subtract the elevate the depth from the elevation of the surface to have the absolute elevation of those bore holes which is the groundwater level so toggle the editing mode and i'm going to the field calculator and there i create a new field which i call groundwater level meters it's a decimal a real keep the other defaults and i go here to fields and values double click on elevation and i do minus and then the depth in meters and the preview looks okay makes it a bit bigger so you can see it and this will give us then the groundwater level in absolute height so i click ok and what we have now here is the groundwater levels for these 11 points toggle off the editing and click save and the next step is that i want to interpolate these groundwater levels to a raster close the attribute table and we're going to use some tools to interpolate the first interpolation method that i'm going to demonstrate is the inverse distance weighing method and there are different ways to calculate it but i'm going to use a processing tool the advantage is that it will calculate it for the full extent instead of using the function that is in here here you also have inverse distance to a power but that will restrict the calculation to the convex hill or boundary of those borehole points so i'm going here to interpolation and i find here idw interpolation i double click i want my borehole set layer as an input and i want there the groundwater level in meters as the layer and then i need to use the plus button to add it here to this list i keep the default quadratic weight quadratic function exponential function there for the the decay function of the weights and i can choose here the extent and that's the advantage of using using this tool instead of the other one so i want to interpolate it to the boundary of the study area click ok and then it will use those coordinates then i'm going to save it here i don't have the option to save it to a geopackage so later if you want it in the geopackage you can drag it from the browser panel here i save it as a tiff file and call it groundwater level idw because later i'm going to show you another method and it will be adjective and then i run oh sorry i forget one more thing to set the pixel size and i wanted at the 30 meter resolution not 0.1 that would be really uh too many pixels for this uh very low amount of points so put it at 30 meters in this case and then i run it and there we see a smooth interpolation that's what the idw does uh it's always good to style this with a color ramp so i'm going to open the layer styling panel remember for continuous data we use single band pseudo color and i'm going to use a ramp and go now to all color ramps and i would like to have something like if it's very deep then it's not so good so that's red and if it's closer to the surface it's blue so i want something from from red to blue so i could use this one or i could use uh this one's in the tutorial or i could use red via yellow to blue so let's use this one in this case and uh there we have it so this is one interpolator um but there's uh yeah it looks maybe nice and smooth but it's very much biased to the points and here we have very sparse distribution of points let me put the points on top of it and uh you could also label those points with the values go to single labels and i choose then the field of the groundwater depth and sorry that's the last field here and then then you can see it just a quick and dirty labeling you can change a lot of these settings um so you see that we only have few points and the distribution is not ideal for this interpolator and i can make that even more visual what works and what what happens here by visualizing this in 3d so i go here to view new 3d map view and that will open the 3d map view and what i'm going to do is change some settings here need to indicate that we use a dem and the dem that we use in this case is the groundwater level idw and i'm going to exaggerate it because you've seen that the elevation values are very little so i go to 10 and then i click ok and then i can see it in 3d and what you see is typical for these kind of interpolators is that it will create bumps and valleys around the points so another interpolator might be better than using this one just to visualize the result idw works very well when we have many points a bit regularly spread over our study area i'm going to show another interpolator that's the decent polygons and for that purpose i can use their the tool from this raster menu go to analysis and i choose grid nearest neighbor that's the decent polygon it will simply assign to each pixel the value of the nearest point so choose here the boreholes set and it's important that here under advanced you choose the z field which will be groundwater level in meters and that's all i go here to save the interpolated raster to groundwater level and i call this one tsun i click save and you see that it uses a g dot tool and i click run there it is now it's limited to the boundaries of our points so that's a result of this algorithm and i'm going to style it to compare it single band pseudo color the yellow to blue blue to sorry red to blue and so it's at the same color skill as the other one and yeah normally when i then discuss which interpolator is uh is better uh students normally go for the smooth one but you've seen that those this one also has some strange effects if we look at it in 3d if we go to these and it simply says that the groundwater elevation is the same as the closest point so that's uh what happens with this interpolator so that's basically what i wanted to show you with this tutorial

Channel: Hans van der Kwast

Views: 8,472

Keywords: Groundwater, QGIS, GIS, Hydrogeology, IGRAC, IHE Delft, SRTM, DEM, sampling, interpolation, IDW

Id: CqyK_QsasSI

Length: 46min 56sec (2816 seconds)

Published: Thu Sep 10 2020