require(tidyverse)
require(ggspatial)
require(ggnewscale)
require(terra)
require(sf)
require(raster)
require(tidyterra)
require(elevatr)
require(topoDistance)Raster Cost Distance in R
Overview
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 for these tasks, 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 the
topoDistancepackage
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-familar tidyverse, sf, and terra packages, we’ll be using ggspatial and tidyterra for plotting, and we’re also adding ggnewscale, which gives a function we’ll use to make some cool-looking terrain maps. For computing cost distances, we’ll be using the topoDistance package, and, as mentioned above, we’ll use the elevatr package to access our DEMs. Finally, because some of the packages we will be using still require raster, terra’s predecessor package, we will be using it, as well.
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 seal 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)Mosaicing & ProjectingNote: Elevation units are in meters.Okay, let’s have a look at what we’ve got here:
everestclass : RasterLayer
dimensions : 1414, 2134, 3017476 (nrow, ncol, ncell)
resolution : 0.000164731, 0.000164731 (x, y)
extent : 86.74805, 87.09958, 27.83905, 28.07198 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : file777872d5508d.tif
names : file777872d5508d This is a RasterLayer, which is the raster object type used in the raster package, terra’s predecessor. Over the next few years, the raster package will be discontinued, but in the meantime there are still some R packages that continue to use it. Fortunately, most of the differences between terra and raster are “under the hood,” so the intuitions you’ve built up working with terra will generally continue to serve you well.
For example, to plot a raster object:
plot(everest)
We can also use all the same map algebra tricks that we know from the terra package.
Plotting beautiful terrain maps with terra
Converting from a raster object to a terra object
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.
First things first, we need to convert our DEM from a raster package raster to our familiar terra SpatRaster format, because we’ll be using terra functions for most of the work here. Fortunately, it’s really easy to move back and forth between the two packages:
everest_spatrast <- rast(everest)That’s it! terra uses the same function, rast() to bringing in data, whether we are reading it from a file, or just converting it from the raster package.
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. If the raster is in latitude and longitude, as is the case here, the elevation must be in meters. As you can see from the output when we downloaded the DEM, the elevation is indeed in meters, and you might also have noticed when we looked at our raster object that its CRS is WGS 1984, so we should be good to go.
Here’s how we use the terrain() function:
everest_slope <- terra::terrain(everest_spatrast, v = "slope",
unit = "degrees", neighbors = 8)
everest_aspect <- terra::terrain(everest_spatrast, 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)
Task 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?
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.
Task 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. The output should look something like this:
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_spatrast) +
scale_fill_hypso_c()<SpatRaster> resampled to 500544 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_spatrast) +
scale_fill_hypso_c(palette = "gmt_globe_hypso")<SpatRaster> resampled to 500544 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_spatrast, 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 500544 cells.
<SpatRaster> resampled to 500544 cells.
Not bad, but we can probably do better.
Task 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.
Topographic distance in the topoDistance package
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 distance. Basically, topographic distance measures the shortest path between points, taking into consideration changes in elevation.
Because topographic distance only considers elevation, all you need to use it is a DEM and a set of locations. For our locations, we just need a matrix of latitude and longitude coordinates. We’ll use the location of Everest Base Camp and the mountain’s summit.
base_camp <- c(86.85, 28)
summit <- c(86.93, 27.99)
locations <- rbind(base_camp,
summit)
colnames(locations) <- c("longitude", "latitude")
locations longitude latitude
base_camp 86.85 28.00
summit 86.93 27.99topoDistnace really only has one function, topoDist(), that does the entire job of building a cost surface and finding the cost paths for you. All it needs are the DEM, the locations, whether or not you want it to return the actual least-cost paths, and how much it should weight the differences in height, as compared to horizontal distance (it defaults to equal weights).
topoDistance is built to use the raster package, so we’ll want to use our original everest raster object as the DEM. Be patient - this is a computationally intensive operation and may take a while:
base_camp_dist <- topoDist(DEM = everest, pts = locations,
paths = TRUE, zweight = 1)Let’s have a look at what we’ve got:
base_camp_dist[[1]]
base_camp summit
base_camp 0.000 9511.112
summit 9511.112 0.000
[[2]]
class : SpatialLines
features : 1
extent : 86.84993, 86.92999, 27.99003, 28.00008 (xmin, xmax, ymin, ymax)
crs : NA This gives us a list object with two components. First, we have a distance matrix that shows the topographic distance between the locations we fed the function. Then we have a SpatialLines object, which is a line vector object from the predecessor of the sf package. We can extract that line object from the list, convert it into an sf object, and then plot it on top of the DEM to see the path:
base_path <- base_camp_dist[[2]] %>% #Gets the 2nd object in a list
st_as_sf() #Converts to an sf object
st_crs(base_path) <- "epsg:4326" #Specifying that we are in WGS 1984
ggplot() +
geom_spatraster(data=everest_shade, show.legend = FALSE) +
scale_fill_distiller(palette = "Greys") +
new_scale_fill() +
geom_spatraster(data=everest_spatrast, alpha = 0.7) +
scale_fill_hypso_tint_c() +
geom_sf(data=base_path, color="yellow") +
labs(fill = "Elevation (masl)",
title = "Shortest Topographic Path from Everest Base Camp to Summit") <SpatRaster> resampled to 500544 cells.
<SpatRaster> resampled to 500544 cells.
Task 4: Compute two more least topographic distance paths between base camp and the Mount Everest summit. For each, set the relative weight of changes in height (zweight) to be greater than 1 (remember that this is a multiplier, so you probably don’t want to set the values too high). Make a map showing all three paths, each in a different color. Use the subtitle of the map to identify which color corresponds to which weighting.
Habitat-based cost distance calculation in topoDistance
topoDistance also allows you to add in other information about the difficulty of moving across a given pixel in the form of a conductance raster, which is really just the inverse of a cost raster (R tends to use conductance rasters because it is less computationally intensive to do so).
You’ll remember that a cost raster’s pixels measure the difficulty of moving across the surface. Because in this tutorial we’ll be building on the analysis conducted by Cisneros-Araujo, et al. (2020), and because they constructed their cost raster based on the types of land cover that might be easier or more difficult for elephants to traverse, we can use their land-cover raster, along with the costs that they assigned to the different land-cover classes, to create our conductance raster.
So let’s start by reading in their land-cover raster. The original data used the Sentinel satellite and was therefore at a 15-meter resolution. I have upscaled that to 250 meters. This will affect our results, but it will expedite our processing time:
land_cover <- rast("Kilombero_LC_250m.tif")
plot(land_cover, type="classes") #Setting type to classes tells R that
#the raster has discrete values.But what do these values represent? I’ve put together a table that we can use to translate the pixel values. Let’s have a look (the Raster_Code column is the pixel value):
lc_classes <- read.csv("Kilombero_LC_Codes.csv",
stringsAsFactors = FALSE)
lc_classes Habitat Resistance Raster_Code
1 Mosaic 1 4
2 Dry Forest 5 2
3 Moist Forest 5 1
4 Grassland 5 5
5 Wetlands 5 7
6 Croplands 20 6
7 Tree Plantation 20 3
8 Railways 20 11
9 Open Water 100 9
10 Irrigated Land 100 10
11 Urban Area 500 8
12 Roads 500 12The Resistance column is the cost per pixel for an elephant to traverse. Now, since we need to have a conductance raster, we have to compute the inverse of Resistance.
lc_classes <- lc_classes %>%
mutate(Conductance = 1/Resistance) #mutate() allows you to add new
#columns to a data frame or tibble
lc_classes Habitat Resistance Raster_Code Conductance
1 Mosaic 1 4 1.000
2 Dry Forest 5 2 0.200
3 Moist Forest 5 1 0.200
4 Grassland 5 5 0.200
5 Wetlands 5 7 0.200
6 Croplands 20 6 0.050
7 Tree Plantation 20 3 0.050
8 Railways 20 11 0.050
9 Open Water 100 9 0.010
10 Irrigated Land 100 10 0.010
11 Urban Area 500 8 0.002
12 Roads 500 12 0.002So the next thing we need to do is convert the pixel values in the land-cover raster to their corresponding Conductance values. terra has a handy function, classify(), that can take care of that for us. It’s a flexible function that can either recode continuous valued rasters into discrete bins or reclassify discrete value rasters like the one here, where the numbers represent different land-cover features.
To reclassify unique numbers into new ones, classify() needs a two-column matrix, where the first column is the value in the raster and the second is the new value we want to assign. That means we need to get just the Raster_Code and Conductance columns from lc_classes and then convert it into a matrix. Here we go:
reclass_mat <- lc_classes[,c("Raster_Code", "Conductance")] %>%
as.matrix()
reclass_mat Raster_Code Conductance
[1,] 4 1.000
[2,] 2 0.200
[3,] 1 0.200
[4,] 5 0.200
[5,] 7 0.200
[6,] 6 0.050
[7,] 3 0.050
[8,] 11 0.050
[9,] 9 0.010
[10,] 10 0.010
[11,] 8 0.002
[12,] 12 0.002Then, to reclassify the raster, we just need to tell classify() which raster object and which reclassification matrix to use:
land_cover_conductance <- classify(x = land_cover,
rcl = reclass_mat)
plot(land_cover_conductance)
Now that we have our conductance raster, we need the points that we’ll be examining. I’ve sort of eyeballed them, based on the article, because the points that they used for suitable elephant habitat are unfortunately not posted online. Let’s create the location matrix to feed into topoDistance:
lats <- c(-8, -8.59, -9.17, -9.43, -9.98)
lons <- c(36.46, 36.86, 35.76, 36.1, 35.82)
habitats <- cbind(lons, lats)
colnames(habitats) <- c("longitude", "latitude")
habitats longitude latitude
[1,] 36.46 -8.00
[2,] 36.86 -8.59
[3,] 35.76 -9.17
[4,] 36.10 -9.43
[5,] 35.82 -9.98Now, these points are in longitude and latitude, so we’ll need to convert them to an sf object and reproject them into the land_cover_conductance object’s CRS. Here’s how:
habitats <- habitats %>%
as_tibble() %>% #sf can't use a plain matrix
st_as_sf(coords = c("longitude", "latitude"),
crs = "epsg:4326") %>% #For data with point coordinates, this
#Tells sf where to find the locations and then specifies the CRS.
#Remember that the coordinate order must be in XY order!
st_transform(st_crs(land_cover_conductance)) %>%
st_coordinates() #Extract the reprojected coordinates
habitats X Y
[1,] 881513.7 9114097
[2,] 925044.3 9048335
[3,] 803331.8 8985188
[4,] 840479.7 8956097
[5,] 809196.4 8895479topoDistance also requires us to have a DEM, in addition to a conductance raster, so that means we need to download an DEM for our area. Luckily, because elevatr can take a SpatRaster input, we can use land_cover_conductance to identify the area we want to cover (notice that because we’re using a SpatRaster, we don’t need to specify the prj argument. Also, we can use the clip argument to clip the layer to the boundaries of the SpatRaster). Because we’re working with a much larger area than before, I’m lowering the zoom value from 12 to 9, which is roughly the resolution of the SpatRaster:
land_cover_elev <- get_elev_raster(locations = land_cover_conductance,
src = "aws", z = 9, clip = "bbox")Mosaicing & ProjectingClipping DEM to bboxNote: Elevation units are in meters.plot(land_cover_elev)
Next, because topoDistance still uses the raster package, we need to convert our SpatRaster to a raster:
land_cover_conductance_rast <- raster(land_cover_conductance)Finally, we need to resample our DEM to match our conductance raster’s pixels:
land_cover_elev <- raster::resample(land_cover_elev,
land_cover_conductance_rast)Now we have all we need to get the least-cost paths between our locations using topoDistance. We’ll use a new function, topoLCP, which allows us to input our conductance raster, in addition to our DEM and points of interest. For right now, we’re actually going to ignore topography to try to approximate the approach taken in the article upon which we’re building, so we’re setting zweight to zero. Here we go (be patient - even with upscaled rasters, this is computationally intensive!):
lcp_kilombero <- topoLCP(DEM = land_cover_elev,
costSurface = land_cover_conductance_rast,
pts = habitats,
paths = TRUE,
zweight = 0)We can quickly visualize the results using a function called topoPathMap():
topoPathMap(DEM = land_cover_elev, pts = habitats,
topoPaths = lcp_kilombero)
Task 5: Whew! That was a lot. Just to solidfy all the steps we just went through, go back through this section and create a flowchart showing all the steps that were required to produce this result.
Task 6: Convert the least-cost paths above into sf objects and create a polished map showing the paths over the land_cover_conductance SpatRaster. Be sure to include a title and scale bar. Compare your results with the map on Page 437 of Cisneros-Araujo, et al. (2020). Where do you see differences? What might account for them?
Task 7: Compute new least-cost paths using a higher zweight value and create a map of your results. Do they differ much from the results with zweight was set to zero? Why or why not? It may be helpful to compute the slope of the DEM and have a look at that when thinking about your answer.