Tutorial: Spatial Analysis of Accessibility of Wells using Map Algebra

Site: IHE DELFT OPEN COURSEWARE
Course: GIS training for Hydrogeological Applications
Book: Tutorial: Spatial Analysis of Accessibility of Wells using Map Algebra
Printed by: Guest user
Date: Friday, 3 December 2021, 11:00 AM

1. Introduction

With map algebra we can do calculations with raster layers. This is useful for spatial analysis. For example, when we need to evaluate different criteria to find suitable or unsuitable locations we can use map algebra.

After this lesson you will be able to:

  • apply map algebra for raster analysis
  • distinguish Boolean, discrete, and continuous rasters
  • make legends for Boolean, discrete, and continuous maps
  • understand the use of Nodata
  • use logical operators
  • calculate distances from rasters
  • reclassify rasters
  • convert raster to vector points
  • sample raster values with vector points

In this lesson we'll address the following case.
The municipality of the (imaginary) oasis Aïn Kju Dzjis has hired you to analyse which wells are unsuitable for its inhabitants based on the following conditions:

Condition 1:
    The wells should be within 150 meters of houses or roads.

Condition 2:
    No industry, mine, or landfill within 300 meters of the wells.

Condition 3:
    The wells should be less than 40 meters deep.

You will use map algebra to perform the required analysis.

This lesson will follow this workflow:

Flowchart map algebra

The raster layers for this lesson can be downloaded from the main lesson page.


2. Preparation

For this task you are provided with the following raster layers: buildg.tif, roads.tif, dtm.tif, and gwlevel.tif.

We'll first inspect the metadata of these layers. Then we'll setup the QGIS Processing Toolbox.


2.1. Checking the Metadata

1. In the Browser panel right-click on Favorites and choose Add a directory.
add to favorites


2. Choose the folder where you have stored the layers.

By adding a folder to Favorites you can quicker access you project files.

3. Click the + to expand the contents of the folder.

4. Preview the maps and metadata of these raster layers by right-clicking on the layer and choosing Layer Properties...

browser layer properties

The Layer Properties window will open, showing the metadata of the layer.

layer properties

5. Check and interpret the metadata.

2.2. Using the Processing Toolbox

The Processing Toolbox in QGIS provides a lot of tools for processing GIS data. Besides QGIS tools, it also has tools from GDAL, GRASS, and SAGA that are very useful.

1. First activate the Processing Toolbox by choosing Processing | Toolbox from the main menu.

processing toolbox menu

Now the Processing Toolbox panel shows up.

Next, we are going to change a default setting of QGIS. By default, the Processing Toolbox doesn’t use the file name of the output as a layer name, which can be confusing.

2. In the main menu choose Settings | Options | Processing.

3. Expand the General menu by clicking on the plus sign. Then check the box at Use filename as layer name.

use filename as layer name

4. Click OK to close the dialogue.

3. Condition 1: Wells within 150 Meters of Houses or Roads

Now that we've looked at the metadata of our input layers and we've set up the Processing Toolbox, it's time to proceed with analysing the three conditions:

  • Condition 1: Wells within 150 meters of houses or roads
  • Condition 2: No industry, mine, or landfill within 300 meters from wells
  • Condition 3: Wells less than 40 meters deep

For Condition 1 we'll first look at the houses and we'll do the following:

  • Style the buildg raster layer
  • Create a boolean layer with True for houses and False for other buildings
  • Create zones of 150 meters around the houses
Then we'll repeat the steps for roads.


3.1. Style a discrete raster layer

Let's first look at the houses. The houses are a class in the buildg layer. Rasters normally only store values, so we need to assign a proper legend to raster layers.

1. In the Layers panel click on the Open layer styling panel button to open the Layer styling panel.

2. Choose Paletted/Unique values as render type. Keep Random as Color ramp type and click Classify.

3. Use the screenshot below as a reference to add the labels to the class numbers by double-clicking on each label name. Double-click on the color to make it more intuitive as shown in the example.

buildg styling

Because buildg is a discrete raster we used the Paletted/Unique values render type. This renderer classifies each unique integer value found in the raster layer.

3.2. Create a Boolean Layer with True for Houses and False for Other Buildings

If we want to create a Boolean layer with True (1) for houses and False (0) for the other classes in the buildg layer we can use the Raster Calculator.

1. In the main menu go to Raster | Raster Calculator.

2. In the Raster Calculator dialogue double-click on buildg@1, click the = button and type 1.

Now the equation reads buildg@1 = 1, which means: if the buildg@1 layer equals 1 (which is the houses class), then the output layer has True (value 1), else it has False (value 0). @1 means band 1. In our case we only use single band raster layers (multiple bands are used in other applications such as remote sensing).

house raster calculator

3. Call the output layer houses.tif and click OK to perform the calculation.

Following good practice, we are going to style this Boolean layer. For Boolean layers we also use the Paletted/Unique values renderer.

4.  Repeat steps 1 and 2 of lesson 4.1 to style the houses layer and choose intuitive colors for True and False.


3.3. Create Zones of 150 Meters Around the Houses

Now that we have a layer with only houses, we can proceed with the calculation of the zones of 150 m around them. We're going to create a Boolean layer which is True (1) within a zone of 150 m around houses and False (0) further than 150 m from houses.

1. In the main menu go to Raster | Analysis | Proximity (Raster Distance).

proximity menu

2. In the Proximity (Raster Distance) dialogue make sure the houses layer is selected as Input layer. Set the Distance units to Georeferenced coordinates. Type for the Maximum distance to be generated 150 meters and type 1 for Value to be applied to all pixels that are within the -maxdist of target pixels. Set the Output data type to Byte (because we use only 0 and 1) and call the output Proximity map houses150m.tif. Leave the other settings at default.

proximity dialog

3. Click Run, then Close when the calculation is done.

4. Style the Boolean map. Make the True pixels green and the False pixels red.

house zone styled


3.4. Create Zones of 150 Meters Around the Roads

In a similar way we can now calculate the 150 m buffer around the roads.

1. Style the roads layer. Use these classes:

0: no roads
1: dirt road
2: tarmac

2. Repeat the steps used to calculate the zones of 150 m around the houses, but use the roads layer as the Input layer and name the output roads150m.tif.

3. Style the Boolean map. Make the True pixels green and the False pixels red.


4. Condition 2: No Industry, Mine, or Landfill within 300 Meters from Wells

For the second condition we first need to reclassify the buildg layer in such a way that the result is a Boolean map with True (1) for industry, mine, and landfill and False (0) 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 m from industry, mine, and landfill.


4.1. Create a Boolean Raster with True for Industry, Mine, and Landfill, and False for Other Buildings

1. Choose in the Processing Toolbox Raster analysis | Reclassify by table.

reclassify by table menu

The Reclassify table dialogue appears. We are going to reclassify the buildg raster using a lookup table.

2. Fill in the dialogue exactly as in the figure below.

reclassify by table

Here you can identify a nodata value for the output layer for values that are excluded from the lookup table. The Range boundaries define if values are included or excluded from the ranges in a row of the lookup table. Here we don't use ranges, but reclassify each value. Therefore we choose min <= value <= max. For the Output data type we use Int16 because the output values are whole numbers.

3. Then go to Lookup table and click Browse button.

4. Fill in the lookup table as shown in the figure below.
lookup table

5. Click OK and Run. Close the dialogue when the processing is finished.

6. Check the result: 1 for mines, industry, and landfills, 0 for the other classes. Use the Identify tool identify tool buttonand click on the map. On the lower right panel you can find the identify results. It displays the value of the pixel of the selected layer in the Layers panel. You might have to resize the columns to see the pixel values.

As you can see the industry is a Boolean layer, so we need to style it with Paletted/Unique Values renderer as we did before.

7. Style the industry layer.


4.2. Create Zones of 300 Meters Around Industry, Mine, and Landfill

With the Boolean map with True for industry, mine, and landfill we can calculate the distance to these pixels. Because the Proximity (Raster Distance) tool does not allow us to assign values for pixels larger than a threshold, we have to calculate all distances and then use the Raster Calculator to calculate a Boolean map with True for pixels further than 300 meters from industry, mine, and landfill.

1. Open the Proximity (Raster Distance) tool.

2. Make sure that Industry is the Input layer. Keep all other settings as default. Name the Output proximity map inddist.tif.

proximity industry dialog


3. Check the result.

The inddist layer is a continuous raster, with real values representing distances.

4. Make a legend that is appropriate for this raster type using Singleband pseudocolor as render type. Use intuitive colors (e.g. a ramp from red to green).

5. Use the Raster Calculator to calculate inddist@1 >= 300. Call the output layer ind300m.tif.

This again results in a Boolean layer.


6. Style the ind300m layer.

ind300m result


5. Condition 3: Wells Less than 40 Meters Deep

For the last condition we need to identify the wells that are less than 40 m deep.

The gwlevel layer 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).

1. Open the Raster Calculator.

2. Subtract the absolute well depth from the DTM using this calculation: dtm@1 - gwlevel@1. Call the output layer welldepth.tif.

3. Style the welldepth layer with the appropriate renderer for this raster type.

4. Next, calculate in the Raster Calculator a Boolean map with wells less than 40 m deep. Call the output layer notdeep.tif.

5. Style the notdeep layer.

notdeep result

6. Combine the Three Conditions

After calculating the Boolean maps for the three conditions we need to combine them to come to the final result.

1. Use the Raster Calculator to combine the three conditions. Because all Boolean results for the conditions need to be true we have to use the AND operator. Call the output layer accessiblewells.tif.

calculate accessible wells

2. Check the resulting raster layer and style the layer.

7. Report the results

Now we have identified the accessible and inaccessible wells we need to prepare the result map for reporting to the contractor.

In the next steps, we are going to:

  • Convert raster cells of the wells to point vectors for a better visualisation
  • Sample raster values to add attributes from the calculated layers
  • Style the analysis results

7.1. Convert Raster Cells to Point Vectors

To present the end result we can style the accessiblewells layer. However, it is nicer to present the wells as point features on the map. Therefore we need to convert the well pixels to point vectors.

1. In the Processing Toolbox go to Vector creation | Raster pixels to points.

raster pixels to points

2. In the dialogue choose the accessiblewells layer as the input Raster layer. For Field name, enter Accessibility. The raster values will be saved under this field in the attribute table.
Choose wells.shp as the output Vector points layer.

raster pixels to points dialog

3. Click Run, then Close when the conversion is finished.

7.2. Sample Raster Values

The point vector layer wells now only contains the field Accessibility. It is however, more informative to also include other data in the attribute table. With the Point sampling tool plugin we can sample the raster layers in this project and add that information to the point attribute table.

1. Install the Point sampling tool plugin.

2. In the Layers panel only check the boxes before the layers you want to sample from and uncheck the others. Choose the following layers: dtm, welldepth, gwdepth, notdeep, ind300m, roads150m, and houses150m.
3. Click the Point sampling tool button point sampling tool button.

4. In the General tab of the Point Sampling Tool dialogue choose wells.shp as the layer containing sampling points, select all Layers with fields/bands to get values from and save the Output point vector layer to wells_final.shp.

point sampling tool general

5. Click the Fields tab.

6. Edit the name of the attribute that will be given to the output field if needed.

point sampling tool fields

7. Click OK and Close.

8. Open the attribute table of wells_final and check the result.
wells final attribute table


7.3. Style the Analysis Results

Now we can style the layers.

The first step in showing the results will be to style the wells_final data by whether the wells are accessible or not.

1. In the Layers panel click on the Open layer styling panel buttonto open the Layer styling panel. Set wells_final as the target layer.

2. Change from a Single symbol to a Categorized renderer. Then set the Column to Accessible. Click Classify.

3. When using this renderer QGIS will create a category to capture any NULL values. Here there are no wells with NULL values in the Accessible column. Select the entry where the Value reads all other values and click the Delete minus button button to remove it. Give the wells with a value of 1 a green symbol with a size of 4 mm. Give the wells with a value of zero a red color with a size of 2 mm.

4. In the Legend column rename the entry with a Value of 1 to Accessible and the entry with a Value of zero to Inaccessible.

styling wells

5. Scroll down in the Layer styling panel and find the Layer Rendering section. Expand it and select Draw effects. The Customize effects customize effects iconicon becomes active. Click on it to open the Effects properties panel. Here you can add Inner and Outer Glows, Drop Shadows, and other effects. Click the box next to Drop Shadow. Select the Drop Shadow effect so that you are seeing the parameters of the Drop Shadow. Reduce the Offset distance to 1.0 mm. Click the Go back go back buttonbutton to return to the main styling panel.

styling dock labels

6. Switch to the Labels tab labeling settings of the Layer styling panel.

7. Change the setting from No labels to Single labels. Set the Label with column to welldepth.

You will now use a labeling expression very similar to that used at the end of the "Preparing Data from Hardcopy Maps" chapter where we labelled the peaks. You will set the expression up so the labels read something like: Well Depth: 76.4 m.

8. Click the Expression expression buttonbutton to open the Expression Dialog window.

9.To begin, add the string Well Depth in front of the well depth value by adding 'Well Depth: ' before the existing expression.

10. To accomplish the second component you will use the format_number function. Use the search box to find the format_number function. Insert it right before the welldepth field. Remember, this function requires a number and a number of decimal places. The number will be the  welldepth field and the number of places 2.

11. Use what you learned about the String Concatenation String concatenationoperator and the the New line new line operatoroperator to nicely format the label and add the units to the end.

12. Your expression should resemble 'Well Depth:'|| '\n' || format_number(welldepth,2) || ' m'.

13. Set the Font to a sans serif font with a Size of 9 points.

14. Switch to the Placement label placement tab icontab and set the Placement to Offset from point. Select the lower right Quadrant placement with an Offset X Y setting of 1.5 mm each.

well label placement

15. Finally click the Automated placement settings Automated placement settings button to open the Automated Placement Engine window. Uncheck the box for Allow truncated labels on edges of map option. This will prevent labels from being cut off.

16. To complete the styling you will work with the dtm layer. Make this layer the target layer in the Layer styling panel. Click the drop-down for the Color ramp and choose Create new color ramp.

17. The Color ramp type window opens. Choose Catalog:cpt-city as the type. Click OK.

18. The Cpt-city Color Ramp window will open. Select the Topography category.

19. Choose sd-a and click OK.

Did you know that you can save this color ramp to your style library? To do this access the color ramp context menu and select Save color ramp. The Save New Color Ramp window will open allowing you to give it a Name and provide Tag(s).

20. Scroll down to the Layer Rendering section and set the Blending mode back to Normal.

21. Finally right-click on the dtm layer and choose Duplicate from the context menu.

22. Drag the dtm copy layer above the dtm layer and turn it on.

23. Make this duplicated dtm layer the target layer in the Layer styling panel. Change the renderer to Hillshade. Scroll down to the Layer Rendering section and set the Blending mode back to Multiply.

wells final


8. Conclusions

In this lesson you have learned to:

  • apply map algebra for raster analysis
  • distinguish Boolean, discrete, and continuous rasters
  • make legends for Boolean, discrete, and continuous maps
  • understand the use of Nodata
  • use logical operators
  • calculate distances from rasters
  • reclassify rasters
  • convert raster to vector points
  • sample raster values with vector points
You can watch this webinar for the full procedure: