Spatial Analysis of Accessibility of Wells using Map Algebra
| Site: | OpenCourseWare for GIS |
| Course: | GIS training for Hydrogeological Applications |
| Book: | Spatial Analysis of Accessibility of Wells using Map Algebra |
| Printed by: | Guest user |
| Date: | Thursday, 21 November 2024, 7:15 PM |
Description
In this tutorial we'll apply map algebra to determine the accessibilty of wells. This tutorial was developed together with Kurt Menke (Bird's Eye View GIS).
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
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:
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.
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...
The Layer Properties window will open, showing the metadata of the layer.
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.
Now the Processing Toolbox panel shows up.
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
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 
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.
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).
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).
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.
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.
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.
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.
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 
.
4. Fill in the lookup table as shown in the figure below.
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 
and
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.
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.
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.
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.
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.
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.
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 
.
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.
5. Click the Fields tab.
6. Edit the name of the attribute that will be given to the output field if needed.
7. Click OK and Close.
8. Open the attribute table of wells_final and check the result.
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 
to 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 
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.
5. Scroll down in the Layer styling panel and find the Layer Rendering section. Expand it and select Draw effects. The Customize effects 
icon 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 
button to return to the main styling panel.
6. Switch to the Labels tab 
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 
button 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 
operator and the the New line 
operator 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 
tab 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.
15. Finally click the 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.
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