Mapping Travel Times with malariaAtlas and Friction Surfaces
Let’s say I’m somewhere in Colorado, and I want to go up a mountain. A big one, over 14,000 feet. About how long would it take, using the fastest overland travel methods available, for me to get to the top of one of those peaks?
An exciting new set of tools from geospatial modelers at the Malaria Atlas Project (MAP) have made it easy for me to find out. Early this year, the MAP team came out with a global map of travel time to cities — that is, for each 1km square chunk of land on the Earth’s surface, how many minutes on average it would take to get to the most-accessible city. (Notice that throughout this post I use the term “most accessible” instead of “nearest”. One of the key points of this research is that places a longer distance away may take less time to reach than places which are physically closer — for example, using public transit, it takes much less time to get from JFK to the Central Park Zoo than to the Queens Botanical Garden, even though the latter is closer as the crow flies.)
The map of travel times to cities is based on a versatile dataset called a friction surface. This is also a map of every 1km-square grid cell of the globe, but what it represents is a relative measure of how “hard” it is to cross that grid cell, based on whether the cell contains good roads, worse roads, railroads, rivers, various water bodies, or terrain with some slope. After generating this surface, the paper authors use it to generate a map of travel times to cities, but you could use it to generate travel times to anywhere you want, from anywhere you want. In this post, I’ll walk you through the process by calculating travel time from any pixel in Colorado to the most-accessible tall mountain peak in the state.
This task is extra-easy thanks to the malariaAtlas package, an R library recently released by MAP which provides downloadable versions of every publicly available MAP surface, along with tools to clip those surfaces to national or subnational shapefiles and make quick, attractive plots. The code in this tutorial is modified from the friction-mapping script made available with the paper release. The code and data needed to reproduce these examples are available from my GitHub.
Travel Time Tutorial
## Packages
library(gdistance)
library(abind)
library(rje)
library(ggplot2)
library(malariaAtlas)## Plot defaults
theme_set(theme_minimal(base_size=14))
After loading the necessary libraries, the first step is to use malariaAtlas::getShp to retrieve a state-level shapefile for the US, subset it to just Colorado, and plot it to make sure we’ve done it right. (Note: due to a change in malariaAtlas the below lines of code are not working. A bug fix is underway, but in the meantime you can find a workaround on GitHub):
USA.shp <- malariaAtlas::getShp(ISO = "USA", admin_level = "admin1")
analysis.shp <- USA.shp[USA.shp@data$name=="Colorado",]
plot(analysis.shp, main="Shape for Clipping")
Next we use the malariaAtlas::getRaster function to retrieve the friction surface from MAP’s database. The optional shp argument speeds this process up considerably by retrieving only the pixels contained within our shapefile.
friction <- malariaAtlas::getRaster(
surface = "A global friction surface enumerating land-based travel speed for a nominal year 2015",
shp = analysis.shp)
malariaAtlas::autoplot_MAPraster(friction)
Next, we have to convert the friction surface to a “transition matrix”. We’re trying to run an algorithm that finds the fastest path from one spot on this map to another, but in order to do that it needs to understand how different points are connected. Specifically, it needs a spreadsheet that states, for every pixel in our area of interest, how easy it is to jump to any other pixel in a single step. There are about 405,000 pixels in our area, so this spreadsheet (aka the transition matrix) has dimensions 405,000 x 405,000. Since we’re only allowing travel to a pixel’s eight neighbors in any one “step”, one row of the matrix will only have eight nonzero values. After we generate this matrix, we have to adjust the values a bit to account for the fact that our area of interest is stretched on a 3D sphere (the globe) so it’s not an evenly-spaced grid:
T <- gdistance::transition(friction, function(x) 1/mean(x), 8)
T.GC <- gdistance::geoCorrection(T) We’re almost there! To generate the “travel times” map, you need a list of coordinates for the places you’re trying to travel to. I pulled a list of the tallest mountains in the United States from Wikipedia, subset it down to those at least 14,000 feet high, and saved it as a csv called peaks.csv. We read in this dataset, and use the sp::over function to keep only those peaks inside the bounds of our shapefile:
## Point locations
point.locations <- read.csv(file = points.filename)
names(point.locations) <- c("X_COORD", "Y_COORD", "name")# Keep only point coordinates within the shapefile bounds
coordinates(point.locations) <- ~ X_COORD + Y_COORD
proj4string(point.locations) <- proj4string(analysis.shp)
overlap <- over(point.locations, analysis.shp)
point.locations <- point.locations[!is.na(overlap$gid),]
point.locations is now an object of class SpatialPointsDataFrame, but the accessibility algorithm needs a matrix of lat-longs:
points <- as.matrix(point.locations@coords)Finally, we’re ready to make our travel times map. The function accCost will run an “accumulated cost surface” algorithm to efficiently calculate travel time to the nearest peak:
access.raster <- gdistance::accCost(T.GC, points)The autoplot_MAPraster function will plot this surface as well — I’ve added a few extra lines of code to overlay our point locations and change the color scheme:
p <- malariaAtlas::autoplot_MAPraster(access.raster,
shp_df=analysis.shp, printed=F)full_plot <- p[[1]] + geom_point(data=data.frame(point.locations@coords),
aes(x=X_COORD, y=Y_COORD)) +
scale_fill_gradientn(colors = rev(rje::cubeHelix(gamma=1.0,
start=1.5,
r=-1.0,
hue=1.5,
n=16)),
name="Minutes \n of Travel") +
ggtitle("Travel Time to Most Accessible Peak") +
theme(axis.text=element_blank(),
panel.border=element_rect(fill=NA, color="white"))print(full_plot)
The final product has a lot of super interesting features. As one would expect, it’s easier to access a tall mountain if you’re already close to one, but there are a number of highly inaccessible locations very near peaks. For example, it’s estimated to take 500 minutes to get from the Pagosa Springs area to the most accessible peak, even though Windom Peak is less then 50km away:
Colorado’s road networks are clearly visible on the map (quite square to the east, more meandering to the west), as are a number of areas likely accessible only on foot on the west side of the mountains.
Discussion and Caveats
Of course, under no circumstances should you use these travel time estimates for mountaineering. The friction surface has no knowledge of trails, or trees, or ice and snow. While 1 kilometer square is an extremely fine grid size for the global analysis this surface was designed for, it may mask considerable variation in terrain within a single grid cell. And it’s directionless — so far as this analysis is concerned, you’re always going uphill. But this general methodology has limitless applications, and this tutorial should give you the tools to start exploring them. Be aware that these calculations take a very long time to run for large areas — the original script gives some tips for speedup.
To give a few other illustrative examples, the plots below show accessibility maps to a particular point on the Peruvian-Ecuadorian border, and to one of ten health facilities in a district of Mozambique:
Thanks to Daniel Weiss and Rebecca Stubbs for valuable edits and suggestions.
