9.4. Lesson: Supplementary Exercise
In this lesson, you will be guided through a complete GIS analysis in QGIS.
Nota
Lesson developed by Linfiniti and S Motala (Cape Peninsula University
of Technology)
9.4.1. Problem Statement
You are tasked with finding areas in and around the Cape Peninsula that are a
suitable habitat for a rare fynbos plant species. The extent of your area of
investigation in the Cape Peninsula is: south of Melkbosstrand, west of Strand.
Botanists have provided you with the following preferences exhibited by the
species in question:
- It grows on east facing slopes.
- It grows on slopes with a gradient between 15% and 60%.
- It grows in areas that have a total annual rainfall of > 1200 mm.
- It will only be found at least 250 m away from any human settlement.
- The area of vegetation in which it occurs should be at least 6000m2 in area.
As a volunteer for Cape Nature, you have agreed to search for the plant on the
closest suitable piece of land to your house. Use your GIS skills to determine
where you should go to look.
9.4.2. Solution Outline
In order to solve this problem, you will have to download data from
https://www.dropbox.com/s/q5evvkizuunrcs0/more_analysis.zip?dl=0 and place it
under exercise_data/more_analysis folder.
You will use it to find the candidate area
that is closest to your house. If you don’t live in Cape Town (where this
problem is based) you can choose any house in the Cape Town region. The
solution will involve:
- analysing the DEM to find the east facing slopes and the slopes with the
correct gradients;
- analysing the rainfall raster to find the areas with the correct amount of
rainfall;
- analysing the Zoning vector layer to find areas that are away from human
settlement and are of the correct size.
9.4.3. Setting up the Map
- Click on the “CRS status” button in the extreme lower right corner of the
screen. Under the CRS tab of the screen that appears, you will
see the box Coordinate reference systems of the world.
- In this box, navigate to .
- Select the entry WGS 84 / UTM zone 33S (with the EPSG code
32733).
- Click OK. The map is now in the UTM33S coordinate
reference system.
- Save the map by clicking on the Save Project As toolbar button,
or use the menu item.
- Save the map in a directory called Rasterprac that you should create
somewhere on your computer. You will save whatever layers you create in this
directory as well.
9.4.4. Loading Data into the Map
In order to process the data, you will need to load the necessary layers
(street names, zones, rainfall, DEM) into the map canvas.
9.4.4.1. For vectors ...
- Click on the Add Vector Layer button, or use the
menu item.
- In the dialog that appears, ensure that the File radio button is
selected.
- Click on the Browse button.
- In the dialog that appears, open the
exercise_data/more_analysis/streets directory.
- Select the file Street_Names_UTM33S.shp.
- Click Open.
The dialog closes and shows the original dialog, with the file path specified
in the text field next to the Browse button. This allows you to
ensure that the correct file is selected. It is also possible to enter the file
path in this field manually, should you wish to do so.
- Click Open. The vector layer will load in your map. Its color is
automatically assigned. It will be changed later.
- Rename the layer to Streets.
- Right-click on it in the Layers list (by default, the pane along
the left-hand side of the screen).
- Click Rename in the dialog that appears and rename it, pressing
the Enter key when done.
- Repeat the vector adding process, but this time select the
Generalised_Zoning_Dissolve_UTM33S.shp file in the
Zoning directory.
- Rename it to Zoning.
9.4.4.2. For rasters ...
- Click on the Add Raster Layer button, or use the
menu item.
- Navigate to the appropriate file, select it, and click Open.
- Do this for each of the two raster files. The files you want are
DEM/reproject/DEM and
Rainfall/reprojected/rainfall.tif.
- Rename the rainfall raster to Rainfall (with an initial capital).
Initially when you load them, the images will be gray rectangles. Don’t
worry, this will be changed later.
- Save the map.
In order to properly see what’s going on, the symbology for the layers needs to
be changed.
9.4.5. Changing the symbology of vector layers
- In the Layers list, right-click on the Streets layer.
- Select Properties from the menu that appears.
- Switch to the Style tab in the dialog that appears.
- Click on the button labelled Change, with a square showing the
current color of the Streets layer.
- Select a new color in the dialog that appears.
- Click OK.
- Click OK again in the Layer Properties dialog. This
will change the color of the Streets layer.
- Follow a similar process for the Zoning layer and choose an
appropriate color for it.
9.4.6. Changing the symbology of raster layers
Raster layer symbology is somewhat different.
- Open the Properties dialog for the Rainfall raster.
- Switch to the Style tab. You’ll notice that this style dialog is
very different from the version used for vector layers.
- Ensure that the button Use standard deviation is selected.
- Change the value in the associated box to 2.00 (it should be set to
0.00 by default).
- Under the heading Contrast enhancement, change the value of the
Current dropdown list to Stretch to MinMax.
- Click OK. The “Rainfall” raster, if visible, should change
colors, allowing you to see different brightness values for each pixel.
- Repeat this process for the DEM, but set the standard deviations used for
stretching to 4.00.
9.4.7. Changing the layer order
- In the Layers list, click and drag layers up and down to change
the order they appear in on the map.
- Newer versions of QGIS may have a Control rendering order
checkbox beneath the Layers list. Ensure that it is checked.
Now that all the data is loaded and properly visible, the analysis can begin.
It is best if the clipping operation is done first. This is so that no
processing power is wasted on computing values in areas that aren’t going to be
used anyway.
9.4.8. Find the Correct Districts
- Load the vector layer admin_boundaries/Western_Cape_UTM33S.shp into
your map.
- Rename it to Districts.
- Right-click on the Districts layer in the Layers list.
- In the menu that appears, select the Query... menu item. The
Query Builder dialog appears.
You will now build a query to select only the following list of districts:
- Bellville,
- Cape,
- Goodwood,
- Kuils River,
- Mitchells Plain,
- Simons Town, and
- Wynberg.
- In the Fields list, double-click on the NAME_2
field. It appears in the SQL where clause text field below.
- Click the = button; an = sign is added to the SQL query.
- Click the All button below the (currently empty)
Values list. After a short delay, this will populate the
Values list with the values of the selected field
(NAME_2).
- Double-click the value Bellville in the Values
list. As before, this will be added to the SQL query.
In order to select more than one district, you’ll need to use the OR
boolean operator.
Click the OR button and it will be added to the SQL query.
Using a process similar to the above, add the following to the existing SQL
query:
Add another OR operator, then work your way through the list of
districts above in a similar fashion.
The final query should be
"NAME_2" = 'Bellville' OR "NAME_2" = 'Cape' OR "NAME_2" = 'Goodwood' OR
"NAME_2" = 'Kuils River' OR "NAME_2" = 'Mitchells Plain' OR "NAME_2" =
'Simons Town' OR "NAME_2" = 'Wynberg'
Click OK. The districts shown in your map are now limited to
those in the list above.
9.4.9. Clip the Rasters
Now that you have an area of interest, you can clip the rasters to this area.
- Ensure that the only layers that are visible are the DEM,
Rainfall and Districts layers.
- Districts must be on top so that they are visible.
- Open the clipping dialog by selecting the menu item .
- In the Input file (raster) dropdown list, select the
DEM layer.
- Specify an output location in the Output file text field by
clicking the Select... button.
- Navigate to your Rasterprac directory.
- Enter a file name.
- Save the file. Leave the No data value checkbox unchecked.
- Use the Extent clipping mode by ensuring the correct radio button
is selected.
- Click and drag an area in the canvas, so that the area which includes the
districts is selected.
- Check the Load into canvas when finished box.
- Click OK.
- After the clipping operation is completed, DO NOT CLOSE the
Clipper dialog. (Doing so would cause you to lose the clipping
area that you have already defined.)
- Select the Rainfall raster in the Input file (raster)
dropdown list and choose a different output file name.
- Do not change any other options. Do not alter the existing clipping area
which you drew previously. Leave everything the same and click
OK.
- After the second clipping operation has completed, you may close the
Clipper dialog.
- Save the map.
9.4.10. Clean up the map
- Remove the original Rainfall and DEM layers from the
Layers list:
- Right-click on these layers and select Remove.
- This will not remove the data from your storage device, it will merely take
it out of your map.
- Deactivate the labels on the Streets layer:
- Click the Labeling button.
- Uncheck the Label this layer with box.
- Click OK.
- Show all the Streets again:
- Right-click on the layer in the Layers list.
- Select Query.
- In the Query dialog that appears, click the Clear
button, then click OK.
- Wait while the data is loaded. All the streets will now be visible.
- Change the raster symbology as before (see Changing the symbology of raster
layers).
- Save the map.
- You can now hide the vector layers by unchecking the box next to them in the
Layers list. This will make the map render faster and will save
you some time.
In order to create the hillshade, you will need to use a plugin that was
written for this purpose.
9.4.11. Activating the Raster Terrain Analysis plugin
This plugin is included by default in QGIS 1.8. However, it may not be
immediately visible. To check if it is accessible on your system:
- Click on the menu item .
- Ensure that the box next to Raster Terrain Analysis plugin is
selected.
- Click OK.
You will now have access to this plugin via the menu item.
Remember that plugins may sometimes depend on certain Python modules being
installed on your system. Should a plugin refuse to work while complaining of
missing dependencies, please ask your tutor or lecturer for assistance.
9.4.12. Create the hillshade
- In the Layers list, ensure that the DEM is the active
layer (i.e., it is highlighted by having been clicked on).
- Click on the menu
item to open the Hillshade dialog.
- Specify an appropriate location for the output layer and call it
hillshade.
- Check the Add result to project box.
- Click OK.
- Wait for it to finish processing.
The new hillshade layer has appeared in your Layers
list.
- Right-click on the hillshade layer in your Layers
list and bring up the Properties dialog.
- Click on the Transparency tab and set the transparency slider to
80%.
- Click OK on the dialog.
- Note the effect when the transparent hillshade is superimposed over the
clipped DEM.
9.4.13. Slope
- Click on the menu item .
- Select the Slope analysis type, with the clipped DEM as the input
layer.
- Specify an appropriate file name and location for output purposes.
- Check the Add result to project box.
- Click OK.
The slope image has been calculated and added to the map. However, as usual it
is just a gray rectangle. To properly see what’s going on, change the symbology
as follows.
- Open the layer Properties dialog (as usual, via the right-click
menu of the layer).
- Click on the Style tab.
- Where it says Grayscale (in the Color map dropdown
menu), change it to Pseudocolor.
- Ensure that the Use standard deviation radio button is selected.
9.4.14. Aspect
- Use the same approach as for calculating the slope, but select
Aspect in the initial dialog box.
Remember to save the map periodically.
9.4.15. Reclassifying rasters
- Click the menu item .
- Specify your Rasterprac directory as the location for the output
layer.
- Ensure that the Add result to project box is selected.
In the Raster bands list on the left, you will see all the raster
layers in your Layers list. If your Slope layer is called
slope, it will be listed as slope@1.
The slope needs to be between 15 and 60 degrees. Everything less
than 15 or greater than 60 must therefore be excluded.
Using the list items and buttons in the interface, build the following
expression:
((slope@1 < 15) OR (slope@1 > 60)) = 0
Set the Output layer field to an appropriate location and file
name.
Click OK.
Now find the correct aspect (east-facing: between 45 and 135
degrees) using the same approach.
Having reclassified all the rasters, you will now see them displayed as gray
rectangles in your map (assuming that they have been added to the map
correctly). To properly display raster data with only two classes (1 and
0, meaning true or false), you will need to change their symbology.
9.4.16. Setting the style for the reclassified layers
- Open the Style tab in the layer’s Properties dialog
as usual.
- Under the heading Load min / max values from band, select the
Actual (slower) radio button.
- Click the Load button.
The Custom min / max values fields should now populate with
0 and 1, respectively. (If they do not, then there was a mistake
with your reclassification of the data, and you will need to go over that part
again.)
- Under the heading Contrast enhancement, set the
Current dropdown list to Stretch To MinMax.
- Click OK.
- Do this for all three reclassified rasters, and remember to save your work!
The only criterion that remains is that the area must be 250m away from
urban areas. We will satisfy this requirement by ensuring that the areas we
compute are 250m or more from the edge of a rural area. Hence, we need
to find all rural areas first.
9.4.17. Finding rural areas
Hide all layers in the Layers list.
Unhide the Zoning vector layer.
Right-click on it and bring up the Query dialog.
Build the following query:
See the earlier instructions for building the Streets query if
you get stuck.
When you’re done, close the Query dialog.
You should see a collection of polygons from the Zoning layer. You
will need to save these to a new layer file.
- On the right-click menu for Zoning, select Save
as....
- Save your layer under the Zoning directory.
- Name the output file rural.shp.
- Click OK.
- Add the layer to your map.
- Click the menu item .
- Select the rural layer as your input vector layer, while leaving
the Use only selected features box unchecked.
- Under Dissolve field, select — Dissolve all —.
- Save your layer under the Zoning directory.
- Click OK. A dialog will appear asking whether you want to add the
new layer to the TOC (“Table of Contents”, referring to the Layers
list).
- Click Yes.
- Close the Dissolve dialog.
- Remove the rural and Zoning layers.
- Save the map.
Now you need to exclude the areas that are within 250m from the edge of
the rural areas. Do this by creating a negative buffer, as explained below.
9.4.18. Creating a negative buffer
- Click the menu item .
- In the dialog that appears, select the rural_dissolve layer as
your input vector layer (Use only selected features should not be
checked).
- Select the Buffer distance button and enter the value -250
into the associated field; the negative value means that the buffer must be
an internal buffer.
- Check the Dissolve buffer results box.
- Set the output file to the same directory as the other rural vector files.
- Name the output file rural_buffer.shp.
- Click Save.
- Click OK and wait for the processing to complete.
- Select Yes on the dialog that appears.
- Close the Buffer dialog.
- Remove the rural_dissolve layer.
- Save the map.
In order to incorporate the rural zones into the same analysis with the three
existing rasters, it will need to be rasterized as well. But in order for the
rasters to be compatible for analysis, they will need to be the same size.
Therefore, before you can rasterize, you’ll need to clip the vector to the same
area as the three rasters. A vector can only be clipped by another vector, so
you will first need to create a bounding box polygon the same size as the
rasters.
9.4.19. Creating a bounding box vector
- Click on the menu item .
- Under the Type heading, select the Polygon button.
- Click Specify CRS and set the coordinate reference system
WGS 84 / UTM zone 33S : EPSG:32733.
- Click OK.
- Click OK on the New Vector Layer dialog as well.
- Save the vector in the Zoning directory.
- Name the output file bbox.shp.
- Hide all layers except the new bbox layer and one of the
reclassified rasters.
- Ensure that the bbox layer is highlighted in the
Layers list.
- Navigate to the menu item and ensure that
Digitizing is selected. You should then see a toolbar icon with a
pencil or koki on it. This is the Toggle editing button.
- Click the Toggle editing button to enter edit mode. This allows
you to edit a vector layer.
- Click the Add feature button, which should be nearby the
Toggle editing button. It may be hidden behind a double arrow
button; if so, click the double arrows to show the Digitizing
toolbar’s hidden buttons.
- With the Add feature tool activated, left-click on the corners of
the raster. You may need to zoom in with the mouse wheel to ensure that it is
accurate. To pan across the map in this mode, click and drag in the map with
the middle mouse button or mouse wheel.
- For the fourth and final point, right-click to finalize the shape.
- Enter any arbitrary number for the shape ID.
- Click OK.
- Click the Save edits button.
- Click the Toggle editing button to stop your editing session.
- Save the map.
Now that you have a bounding box, you can use it to clip the rural buffer
layer.
9.4.20. Clipping a vector layer
- Ensure that only the bbox and rural_buffer layers are
visible, with the latter on top.
- Click the menu item .
- In the dialog that appears, set the input vector layer to
rural_buffer and the clip layer to bbox, with both
Use only selected features boxes unchecked.
- Put the output file under the Zoning directory.
- Name the output file rural_clipped.
- Click OK.
- When prompted to add the layer to the TOC, click Yes.
- Close the dialog.
- Compare the three vectors and see the results for yourself.
- Remove the bbox and rural_buffer layers, then save
your map.
Now it’s ready to be rasterized.
9.4.21. Rasterizing a vector layer
You’ll need to specify a pixel size for a new raster that you create, so first
you’ll need to know the size of one of your existing rasters.
- Open the Properties dialog of any of the three existing rasters.
- Switch to the Metadata tab.
- Make a note of the X and Y values under the heading
Dimensions in the Metadata table.
- Close the Properties dialog.
- Click on the menu item.
You may receive a warning about a dataset being unsupported. Click it away
and ignore it.
- Select rural_clipped as your input layer.
- Set an output file location inside the Zoning directory.
- Name the output file rural_raster.tif.
- Check the New size box and enter the X and Y values
you made a note of earlier.
- Check the Load into canvas box.
- Click the pencil icon next to the text field which shows the command that
will be run. At the end of the existing text, add a space and then the text
-burn 1. This tells the Rasterize function to “burn” the existing
vector into the new raster and give the areas covered by the vector the new
value of 1 (as opposed to the rest of the image, which will
automatically be 0).
- Click OK.
- The new raster should show up in your map once it has been computed.
- The new raster will look like a grey rectangle – you may change the display
style as you did for the reclassified rasters.
- Save your map.
Now that you have all four criteria each in a separate raster, you need to
combine them to see which areas satisfy all the criteria. To do so, the rasters
will be multiplied with each other. When this happens, all overlapping pixels
with a value of 1 will retain the value of 1, but if a pixel has
the value of 0 in any of the four rasters, then it will be 0 in
the result. In this way, the result will contain only the overlapping areas.
9.4.22. Combining rasters
Click the menu item.
Build the following expression (with the appropriate names for your layers,
depending on what you called them):
[Rural raster] * [Reclassified aspect] * [Reclassified slope] *
[Reclassified rainfall]
Set the output location to the Rasterprac directory.
Name the output raster cross_product.tif.
Ensure that the Add result to project box is checked.
Click OK.
Change the symbology of the new raster in the same way as you set the style
for the other reclassified rasters. The new raster now properly displays the
areas where all the criteria are satisfied.
To get the final result, you need to select the areas that are greater than
6000m^2. However, computing these areas accurately is only possible for
a vector layer, so you will need to vectorize the raster.
9.4.23. Vectorizing the raster
- Click on the menu item .
- Select the cross_product raster.
- Set the output location to Rasterprac.
- Name the file candidate_areas.shp.
- Ensure that Load into canvas when finished is checked.
- Click OK.
- Close the dialog when processing is complete.
All areas of the raster have been vectorized, so you need to select only the
areas that have a value of 1.
Open the Query dialog for the new vector.
Build this query:
Click OK.
Create a new vector file from the results by saving the
candidate_areas vector after the query is complete (and only the
areas with a value of 1 are visible). Use the Save as...
function in the layer’s right-click menu for this.
Save the file in the Rasterprac directory.
Name the file candidate_areas_only.shp.
Save your map.
9.4.24. Calculating the area for each polygon
Open the new vector layer’s right-click menu.
Select Open attribute table.
Click the Toggle editing mode button along the bottom of the
table, or press Ctrl+E.
Click the Open field calculator button along the bottom of the
table, or press Ctrl+I.
Under the New field heading in the dialog that appears, enter the
field name area. The output field type should be an integer, and the
field width should be 10.
In Field calculator expression, type:
This means that the field calculator will calculate the area of each polygon
in the vector layer and will then populate a new integer column (called
area) with the computed value.
Click OK.
Do the same thing for another new field called id. In
Field calculator expression, type:
This ensures that each polygon has a unique ID for identification purposes.
Click Toggle editing mode again, and save your edits if prompted
to do so.
9.4.25. Selecting areas of a given size
Now that the areas are known:
You now have your solution areas, from which you will pick the one nearest to
your house.
9.4.26. Digitize your house
- Create a new vector layer as before, but this time, select the
Type value as being a Point.
- Ensure that it is in the correct CRS!
- Name the new layer house.shp.
- Finish creating the new layer.
- Enter edit mode (while the new layer is selected).
- Click the point where your house or other current place of residence is,
using the streets as a guide. You might have to open other layers to help you
find your house. If you don’t live anywhere nearby, just click somewhere
among the streets where a house could conceivably be.
- Enter any arbitrary number for the shape ID.
- Click OK.
- Save your edits and exit edit mode.
- Save the map.
You will need to find the centroids (“centers of mass”) for the solution area
polygons in order to decide which is closest to your house.
9.4.27. Calculate polygon centroids
- Click on the
menu item.
- Specify the input layer as solution.shp.
- Provide the output location as Rasterprac.
- Call the destination file solution_centroids.shp.
- Click OK and add the result to the TOC (Layers list),
then close the dialog.
- Drag the new layer to the top of the layer order so that you can see it.
9.4.28. Calculate which centroid is closest to your house
- Click on the menu item .
- The input layer should be your house, and the target layer
solution_centroids. Both of these should use the id field
as their unique ID field.
- The output matrix type should be linear.
- Set an appropriate output location and name.
- Click OK.
- Open the file in a text editor (or import it into a spreadsheet). Note which
target ID is associated with the shortest Distance. There may be
more than one at the same distance.
- Build a query in QGIS to select only the solution areas closest to your house
(selecting it using the id field).
This is the final answer to the research question.
For your submission, include the semi-transparent hillshade layer over an
appealing raster of your choice (such as the DEM or the
slope raster, for example). Also include the polygon of the closest
solution area(s), as well as your house. Follow all the best practices for
cartography in creating your output map.