require(tidyverse)
require(ggspatial)
require(ggnewscale)
require(terra)
require(sf)
require(tidyterra)
require(elevatr)Raster Cost Distance in R
Overview
Big picture: cost distance as raster-based reasoning
This tutorial extends earlier work on map algebra and raster alignment by treating space as an impedance surface rather than a neutral plane. Each pixel represents a cost of movement, and distance is accumulated across the raster rather than measured “as the crow flies.”
As you work, consider how assumptions about cost, raster resolution, and accumulation influence the conclusions one might draw from spatial models.
In this tutorial, we’ll be looking at raster cost distance - that is, functions that use quantitative measure of the difficulty of traversing terrain to characterize the connectivity (or lack thereof) between locations. We’ll be using a combination of land-cover and topographic data, which will also give us an opportunity to further develop our skills for creating high-quality maps in R
The Objectives
- Learn how to read in and process digital elevation model (DEM) data using the
elevatrpackage - Learn how to recode raster values
- Learn how to identify least-cost paths using
terra’scostDist()andshortestPath()functions
The Data
For our topographical data, we’ll be using the Amazon Web Services (AWS) Open Topography Global Datasets, which we will download using a package called elevatr. The land-cover data we will use come from one of our readings for this week: Cisneros-Araujo, et al. (2020). Remote sensing of wildlife connectivity networks and priority locations for conservation in the Southern Agricultural Growth Corridor (SAGCOT) in Tanzania. Remote Sensing in Ecology and Conservation, 7(3), 430-444.
We have two files that we will be reading in:
- Kilombero_LC_250m.tif is a lower-resolution version of the land-cover dataset used by Cisneror-Araujo, et al. (2020)
- Kilombero_LC_Codes.csv is a crosswalk file that we will use to identify the land-cover types designated by the pixel values in the Kilombero_LC_250m.tif file
These data can be accessed from this zip file.
Loading Packages
In addition to the now-familiar tidyverse, sf, and terra packages, we’ll be using ggspatial and tidyterra for plotting, and we’re also using ggnewscale, which gives a function we’ll use to make some cool-looking terrain maps. We’ll use elevatr to access digital elevation models. Unlike earlier approaches to cost distance in R, we’ll be computing cost distances directly using terra’s own built-in functions, so we don’t need any additional packages for that purpose.
Let’s load these packages in now:
Working with topographic data in R
For the first section of the tutorial, we’ll take some time to get familiar with digital elevation models (DEMs), a common and very useful form of geospatial data. As you probably remember, DEMs are simply rasters whose pixel values show the pixel’s elevation, usually in meters above sea level (masl).
There are several global DEM datasets at various resolutions, and there are a few packages that allow you to download tiles from those datasets within an R session. We’re using the elevatr package for that purpose here, but other packages are available.
Downloading DEMs with elevatr
We’re going to work with an area around Mount Everest, just so we have an interesting topography to play around with. elevatr can download topographic data based on the extent of sf or terra objects, as well as based on points.
To make things simpler, we’re just going to use points that define a rectangle around our area of interest. Here’s how:
latitudes <- c(27.9, 27.9, 28.05, 28.05)
longitudes <- c(86.8, 87.05, 86.8, 87.05)
locations <- cbind(longitudes, latitudes) %>%
as.data.frame()
names(locations) <- c("x", "y")
locations x y
1 86.80 27.90
2 87.05 27.90
3 86.80 28.05
4 87.05 28.0527.9 - 28.05; 86.8 - 87.05
So that gives us a data frame (the old version of a tibble) with the coordinates of our points. Now we can use that to get our topographic data.
We can download a DEM using the get_elev_raster() function. We just need to input values for locations, which could be an sf or terra object, or a data frame like the one we just created. If we are using a list of coordinates, we also need to specify the projection using the prj argument. src tells R where to get the data. Finally, there is an argument z that specifies the zoom (really, the resolution) of the desired raster. The levels of z correspond to the table at this link.
Okay, so let’s download our DEM. I’m setting z to 12, which will give us about a 30-meter resolution DEM:
everest <- get_elev_raster(locations = locations, src = "aws",
prj = "epsg:4326", z = 12) |>
rast() |> #Converting from raster to terra format
project(y = "epsg:32645")#We're setting the projection to the local UTM zone forMosaicing & ProjectingNote: Elevation units are in meters. #Mount Everest
plot(everest)
Plotting beautiful terrain maps with terra
While you don’t necessarily need to be able to make beautiful maps to conduct spatial analysis, it’s just so easy to create some really cool visualizations using the terra and tidyterra packages and a couple additions to ggplot2 that I want to take a brief detour to show you R’s geovisualization capabilities.
Slope, aspect, and hillshade in terra
You probably remember that we need to do some processing of DEMs in order to take advantage of them for sophisticated terrain modeling. Here, we’re going to use a few different functions that people commonly apply to DEMs:
- Slope: Compute a measure of the change in elevation in a pixel (could be as an angle, percentage, or ratio of horizontal to vertical change)
- Aspect: Compute the direction, in degrees or radians, that a pixel faces
- Hillshade: Simulate shadows cast over a landscape from a light source at a specified location
The first two of these can be computed using terra’s terrain() function. Because there is also a terrain function in the raster package, we have to use a special format to tell R which function we want: terra::terrain(). You can use this way of writing things for any function, and it is quite handy when you have a lot of packages loaded, some of which might have functions with the same name, making it easy for R to get confused about which one you want.
terra’s terrain() function can work on either projected or unprojected rasters. If the raster is projected, the elevation units must be in the same units as the projection. We already made sure that everthing is in meters, so we should be good to go.
Here’s how we use the terrain() function:
everest_slope <- terra::terrain(everest, v = "slope",
unit = "degrees", neighbors = 8)
everest_aspect <- terra::terrain(everest, v = "aspect",
unit = "degrees", neighbors = 8)Let’s unpack those functions. The first argument of each is just the DEM that you are using for the computation. v is the type of operation you want to perform. You can look at help("terrain") to see what kinds of options you have. unit only applies to slope or aspect and can be either “degrees” or “radians”.
neighbors specifies whether we should consider all the surrounding pixels, including the diagonals (8) or just the pixels whose sides touch (4). We call the first the “queen’s case” and the second the “rook’s case”:
Queen’s Case (8):
V1 V2 V3
1 * * *
2 * Pixel *
3 * * *Rook’s Case (4):
V1 V2 V3
1 *
2 * Pixel *
3 * Let’s have a look at what we’ve made:
plot(everest_slope)
plot(everest_aspect)
Deliverable 1
Describe the slope and aspect rasters you have created. What do they mean? How could you use the slope map to identify the steepest parts of Mount Everest? How can you use the aspect map to identify the ridgelines at the top of surrounding peaks? (Hint: remember that aspect shows you which way the pixel is facing)
Once we have computed the slope and aspect, we can use them to create a hillshade raster. terra’s shade() function allows us to do this:
everest_shade <- terra::shade(slope = everest_slope,
aspect = everest_aspect,
angle = 30, direction = 315,
normalize = TRUE)The two arguments here that aren’t immediately obvious are angle and direction. angle sets the vertical angle of the simulated light source relative to the horizon. At 30 degrees, this would be something like a late afternoon. Direction is the clockwise angle relative to north, so at 315 degrees, the light source would be coming from the northwest.
Here’s the result:
plot(everest_shade)
Huh . . . something’s definitely wrong! Turns out, shade() requires that the slope and aspect rasters be in radians! We’ll have to fix that.
Deliverable 2
Recreate the everest_slope and everest_aspect rasters, but this time, compute the results in radians. Then recompute the hillshade and plot your result. Be sure to include both the code and your output in your tutorial report. The output should look something like the following.
Ah, much better!
Making terrain maps with ggplot2 and tidyterra
So the hillshade looks kinda cool, but it’s not really clear what we’re supposed to do with it. Nevertheless, we now have all the tools we need to create a kicka– terrain map in ggplot2.
Cool-looking terrain maps are generally made up of two layers carefully stacked and blended. On the bottom is the hillshade, in grayscale. Above it, slightly transparent, is the DEM with a hyposometric color scale. Hypsometric colors are palettes designed for elevation maps, that give a sense of changing terrain with elevation. For example:
ggplot() +
geom_spatraster(data = everest) +
scale_fill_hypso_c()<SpatRaster> resampled to 500409 cells.
So you can see that the effect is to make it look like Mount Everest’s peak is snowcapped. You can look at help("scale_hypso") to see other palettes you can use. For example:
ggplot() +
geom_spatraster(data = everest) +
scale_fill_hypso_c(palette = "gmt_globe_hypso")<SpatRaster> resampled to 500409 cells.
Good start. Now we layer this on top of the hillshade. Remember that ggplot2 draws layers in the order they are run in the code, so we want to list the hillshade layer first. Also, a ggplot2 object usually can have only one scale for each visual variable - so usually only one fill palette. With the ggnewscale package, though, we can add distinct scales for the DEM and the hillshade layers:
ggplot() +
geom_spatraster(data=everest_shade, show.legend = FALSE) +
scale_fill_distiller(palette = "Greys") +
new_scale_fill() +
geom_spatraster(data=everest, alpha = 0.7) +
#alpha is the proportion opaque for the layer, so we're saying to make
#this 30% transparent
scale_fill_hypso_tint_c() +
labs(fill = "Elevation (masl)")<SpatRaster> resampled to 500409 cells.
<SpatRaster> resampled to 500409 cells.
Not bad, but we can probably do better.
Deliverable 3
Create your own more polished version of the terrain map for the area around Mount Everest. Be sure to include a scale bar and title. You might also experiment with using scale_fill_gradient(), instead of the hypsometric scale, to try to capture the fact that Everest and its environs will be even more snowy. Be sure to inlcude both your code and the map you create.
Topographic cost distance with terra
Now that we’ve completed our diversion into terrain mapping, let’s get back to cost distance modeling. One of the simplest examples, which we can use to get a feel for the approach, is topographic cost distance: the shortest path between two points, taking into account that steeper terrain is harder to traverse.
To find this path, we need two things: a set of locations and a cost surface — a raster where each pixel value represents how expensive it is to move through that pixel. Because we want to penalize steep terrain, we’ll use the slope raster we already computed as our cost surface. Pixels with gentle slopes are cheap to cross; pixels with steep slopes are costly. Note that everest_slope at this point is in radians (from Deliverable 2), so values range from just above 0 (flat) to roughly 1.5 (near-vertical cliff faces). Adding 1 gives us cost values ranging from 1 to about 2.5. Let’s assume that terrain difficulty literatally increases exponentially. We’ll add a small constant to the slope raster so that even perfectly flat terrain carries a non-zero traversal cost (representing the baseline cost of horizontal movement). These are probably terrible assumptions that we’re just using for illustrative purposes:
slope_cost <- (everest_slope^3) + 0.25Creating point locations with vect()
For our two locations, we’ll use Everest Base Camp and the mountain’s summit. So far in this class we’ve been using sf to represent point data. terra has its own equivalent vector object type, the SpatVector, which we create using vect(). Here we create one directly from a coordinate matrix:
base_camp <- c(86.85, 28)
summit <- c(86.93, 27.99)
locations <- rbind(base_camp, summit)
colnames(locations) <- c("longitude", "latitude")
locations_sv <- vect(locations,
crs = "epsg:4326") |>
terra::project(y="epsg:32645") # Reproject to match the cost raster
plot(everest)
plot(locations_sv, add=TRUE)
vect() takes several kinds of inputs. Here we pass it a matrix of coordinates. If there only two columns, it assumes the first is the x and the second the y coordinate. The result is a SpatVector of two point features, analogous to an sf points object.
Computing accumulated cost with costDist()
terra’s costDist() function computes the accumulated cost of reaching every cell in the raster from one or more target locations, traveling along paths that minimize cumulative cost. Think of it as asking: starting from Base Camp, what is the cheapest total cost to reach every other pixel on the map?
costDist() requires a specific type of input: a raster where all pixels are the estimated cost of traversing an area, with a special value given to pixels from which we want to compute distances, called the target pixels.
#Set the target pixels
base_camp_sv <- locations_sv[1, ] # Extract just the first point
base_camp_sv <- terra::rasterize(x= base_camp_sv, y = slope_cost) # Convert to a raster matching slope_cost
slope_cost[!is.na(base_camp_sv)] <- -9999 #Add a clear flag to show the target pixels in slope_cost
base_camp_costdist <- costDist(x = slope_cost, target = -9999)
plot(base_camp_costdist)
The output is a raster with the same extent and resolution as slope_cost. Each pixel’s value is the total accumulated cost of the cheapest possible route from Base Camp to that pixel. Notice that valleys and ridges create natural “corridors” of low accumulated cost — the algorithm finds routes that stay on gentler terrain wherever possible.