Multiple-step Least Cost Path Analysis with QGIS and Python
Using Digital Elevation Models to estimate walking routes
In a previous post I wrote about finding and geocoding locations using textual analysis packages from R. In this series of posts, I am focusing on the use of these tools for Digital Humanities, but the methods are broadly applicable and useful in other disciplines, such as archaeology. In this post, I will provide a GIS work-flow for for examining topography in literature, following a methodology that uses cost-surface-analysis and least-cost-path analysis to examine most likely routes taken by writers of historic travelogues .
This tutorial requires some basic Python skills, and familiarity with QGIS.
Cost-Surface-Analysis (CSA)
CSA is a raster based analysis, which entails the creation of a cost grid (raster), where each cell value represents the cost of moving across the cell in space. Cost rasters are user-defined and can be anything from distance, to type of vegetation. In this case, we will use a cost surface that indicates the steepness of a slope derived from a digital elevation model (DEM).
Least-Cost-Path (LCP) Analysis
LCP analysis is a GIS tool that calculates the easiest route between two points based on a cost surface raster. The result will provide a line vector of the route, which can then be used for visualization.
Workflow
There are four steps to this process. Step 0. Generate a vector of points (see previous post). Step 1. Download a DEM. Step 2. Generate a cost-surface raster (in our case, using the slope) Step 3. Calculate the LCP between points.
Step 1: Download DEM
In order to generate a LCP vector, we must first have a DEM, which is a matrix that contains elevations for each grid. There are different sources for DEM downloads, but, in QGIS, DEMs can be downloaded using the SRTM downloader plug-in. This is convenient in that it will load the DEM directly into QGIS as a raster layer. It simply requires an extent that matches your points vector.
If your query results in multiple rasters, you will need to merge them into one layer. The resulting DEM can be quite large. If this is the case, you can reduce the resolution of the DEM by using the Build Pyramids function from the layer properties menu. This will reduce the size and processing requirements, but it will also be less spatially accurate, as your new DEM will have a larger grid size.
Step 2: Align the layers
After downloading the DEM, we will need to calculate our cost raster. In this case, we will use the slope of the terrain to calculate our cost value. Essentially, slope is calculated by finding the difference between adjacent matrix cells- the larger the difference, the steeper the slope. Luckily for us, there is a built in slope calculate in QGIS.
However, before calculating the slope, we need to be sure that our raster is in the proper projection- particularly, we need to be sure it is projected to the correct UTM zone. If you don’t know what UTM zone your location is in, there are many tools online to find UTM zones based on coordinates.
Once you have your zone, you will need to reproject the raster using the reproject tool found in Raster > Projection > Warp (Reproject). In the box, at the Target CRS option, select the icon at the right of the box. Search for the projection using the search bar. There will be a few different options, but select the option that matches the CRS of the coordinates vector (most likely WGS84 / UTM##). It is important that the coordinates vector and the raster layers are in the same projections, however you don’t need to project the coordinates vector into the correct UTM zone.
If you try to generate your slope layer from a DEM that is not projected into the correct UTM zone, the result will be unusable.
Step 3: Calculate the slope
After aligning all of the layers, we need to calculate the slope using the slope algorithm in QGIS. If we have done the previous step correctly, the result should look something like this:
The slope calculation gives cells a value between 0 and 90, which signifies the degree of the slope between the cells. We need to reassign these values to classes in order to run the LCP algorithm. To do this, right click on the slope layer and select the Properties > Symbology menu.
Here we can render new values from the single gray band. At the top (1) we need to change the render type to singleband pseudocolor. This will generate more options, where we can change the interpolation from linear to discrete (2), and finally change the mode to quantile (3), you may also choose to increase the number of classes, which will make the LCP more accurate. The result will group values and assign them a new color (something that you can change if you prefer a more topographical color palette).
Step 4: Run the function
We are finally ready to generate the LCP vector. While there are different ways of doing this, I have chosen to use the Cost Distance Analysis: Least Cost Path tool- a plugin that is also available from the plugin menu. While this function has an option to iterate, it does not have an option to do so sequentially. Instead, I have written a simple function that will break your coordinates vector into pairs to pass to the algorithm programatically. As I am not very expereienced with Python, and even less so with QtPy, this function may be a little bulky, however, it works.
You can copy and paste this function directly into a new Python script within the QGIS Python console (Ctrl + Alt+ P) and then provide it with your own variables. The result will be an individual vector for each point pair, which can then be joined within QGIS for further analysis or visualization.
Conclusion
In this post I have documented a relatively simple work-flow for calculating LCP from a slope raster. All tools and data used here can be accessed for free from within QGIS. I have also provided a custom Python function for step-wise iteration between vector points. In the final post, I will put everything together to show how to visualize the results of a LCP analysis.
Originally published at https://greggsaldutti.netlify.app on February 22, 2021.
