Creating a figure of map layers in R
I like figures of map layers to illustrate the many different types of data sets
we combine to do urban and transport modeling. And oftentimes I get obsessed with like making maps that are reproducible with code in R. In this post I’ll be sharing a reproducible example showing how to create a figure of stacked maps like this one below.
Quick background:
In 2014, I was trying to find a way to create map layers in R. This was before the sf library was created. Most of us were using the sp library for handling spatial data and
Barry Rowlingson was super helpful, as usual. I used Barry’s suggestion to create a
reproducible example so I could use it latter, but then sf was created and it completely changed how we do spatial analysis in R. Since then,
Lauren O’brien
proposed a simple way to
tilt and stack sf objects and
Stefan Jünger created
an elegant function to do this. I’ll be using Stefan’s function in my example below.
Load libraries
library(easypackages)
easypackages::packages("sf",
"raster",
"stars",
"r5r",
"geobr",
"aopdata",
"gtfs2gps",
"ggplot2",
"osmdata",
"h3jsr",
"viridisLite",
"ggnewscale",
"dplyr",
"magrittr",
prompt = FALSE
)
Functions to tilt sf
Original function created by Stefan Jünger.
rotate_data <- function(data, x_add = 0, y_add = 0) {
shear_matrix <- function(){ matrix(c(2, 1.2, 0, 1), 2, 2) }
rotate_matrix <- function(x){
matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2)
}
data %>%
dplyr::mutate(
geometry = .$geometry * shear_matrix() * rotate_matrix(pi/20) + c(x_add, y_add)
)
}
rotate_data_geom <- function(data, x_add = 0, y_add = 0) {
shear_matrix <- function(){ matrix(c(2, 1.2, 0, 1), 2, 2) }
rotate_matrix <- function(x) {
matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2)
}
data %>%
dplyr::mutate(
geom = .$geom * shear_matrix() * rotate_matrix(pi/20) + c(x_add, y_add)
)
}
Load data
We’ll be using a few data sets available from the packages used here. The first thing we need to do is to load the data and crop them to make sure they have the same extent.
### get terrain data ----------------
# read terrain raster and calculate hill Shade
dem <- stars::read_stars(system.file("extdata/poa/poa_elevation.tif", package = "r5r"))
dem <- st_as_sf(dem)
# crop
bbox <- st_bbox(dem)
### get public transport network data ----------------
gtfs <- gtfs2gps::read_gtfs( system.file("extdata/poa/poa.zip", package = "r5r") )
gtfs <- gtfs2gps::gtfs_shapes_as_sf(gtfs)
# crop
gtfs <- gtfs[bbox,]
gtfs <- st_crop(gtfs, bbox)
plot(gtfs['shape_id'])
### get OSM data ----------------
# roads from OSM
roads <- opq('porto alegre') %>%
add_osm_feature(key = 'highway',
value = c("motorway", "primary","secondary")) %>% osmdata_sf()
roads <- roads$osm_lines
# crop
roads2 <- roads[bbox,]
roads2 <- st_crop(roads2, bbox)
plot(roads2['osm_id'])
### get H3 hexagonal grid ----------------
# get poa muni and hex ids
poa <- read_municipality(code_muni = 4314902 )
hex_ids <- h3jsr::polyfill(poa, res = 7, simple = TRUE)
# pass h3 ids to return the hexagonal grid
hex_grid <- h3jsr::h3_to_polygon(hex_ids, simple = FALSE)
plot(hex_grid)
# crop
hex_grid <- hex_grid[bbox,]
hex <- st_crop(hex_grid, bbox)
plot(hex)
### get land use data from AOP project ----------------
#' more info at https://www.ipea.gov.br/acessooportunidades/en/
landuse <- aopdata::read_access(city = 'poa', geometry = T, mode='public_transport')
# crop
landuse <- landuse[bbox,]
landuse <- st_crop(landuse, bbox)
plot(landuse['CMATT30'])
# hospitals
# generate one point per hospital in corresponding hex cells
df_temp <- subset(landuse, S004>0)
hospitals <- st_sample(x = df_temp, df_temp$S004, by_polygon = T)
hospitals <- st_sf(hospitals)
hospitals$geometry <- st_geometry(hospitals)
hospitals$hospitals <- NULL
hospitals <- st_sf(hospitals)
plot(hospitals)
# schools
# generate one point per schools in corresponding hex cells
df_temp <- subset(landuse, E001>0)
schools <- st_sample(x = df_temp, df_temp$E001, by_polygon = T)
schools <- st_sf(schools)
schools$geometry <- st_geometry(schools)
schools$schools <- NULL
schools <- st_sf(schools)
plot(schools)
Plot
### plot ----------------
# annotate parameters
x = -141.25
color = 'gray40'
temp1 <- ggplot() +
# terrain
geom_sf(data = dem %>% rotate_data(), aes(fill=poa_elevation.tif), color=NA, show.legend = FALSE) +
scale_fill_distiller(palette = "YlOrRd", direction = 1) +
annotate("text", label='Terrain', x=x, y= -8.0, hjust = 0, color=color) +
labs(caption = "image by @UrbanDemog")
temp2 <- temp1 +
# pop income
new_scale_fill() +
new_scale_color() +
geom_sf(data = subset(landuse,P001>0) %>% rotate_data(y_add = .1), aes(fill=R001), color=NA, show.legend = FALSE) +
scale_fill_viridis_c(option = 'E') +
annotate("text", label='Population', x=x, y= -7.9, hjust = 0, color=color) +
# schools
geom_sf(data = hex %>% rotate_data(y_add = .2), color='gray50', fill=NA, size=.1) +
geom_sf(data = schools %>% rotate_data(y_add = .2), color='#0f3c53', size=.1, alpha=.8) +
annotate("text", label='Schools', x=x, y= -7.8, hjust = 0, color=color) +
# hospitals
geom_sf(data = hex %>% rotate_data(y_add = .3), color='gray50', fill=NA, size=.1) +
geom_sf(data = hospitals %>% rotate_data(y_add = .3), color='#d5303e', size=.1, alpha=.5) +
annotate("text", label='Hospitals', x=x, y= -7.7, hjust = 0, color=color) +
# OSM
geom_sf(data = roads2 %>% rotate_data(y_add = .4), color='#019a98', size=.2) +
annotate("text", label='Roads', x=x, y= -7.6, hjust = 0, color=color) +
# public transport
geom_sf(data = gtfs %>% rotate_data(y_add = .5), color='#0f3c53', size=.2) +
annotate("text", label='Public transport', x=x, y= -7.5, hjust = 0, color=color) +
# accessibility
new_scale_fill() +
new_scale_color() +
geom_sf(data = subset(landuse, P001>0) %>% rotate_data(y_add = .6), aes(fill=CMATT30), color=NA, show.legend = FALSE) +
scale_fill_viridis_c(direction = 1, option = 'viridis' ) +
theme(legend.position = "none") +
annotate("text", label='Accessibility', x=x, y= -7.4, hjust = 0, color=color) +
theme_void() +
scale_x_continuous(limits = c(-141.65, -141.1))
# save plot
ggsave(plot = temp2, filename = 'map_layers.png',
dpi=200, width = 15, height = 16, units='cm')
