1 - Principles
the mapextrud R package allows to make extruded maps from sf objects (polygons or multipolygons). For the moment, it is only available on github.
Until a more advanced version is available on CRAN, you can install the package by writing this…
… and load the package like that.
2 - Basemap setup
The first thing to do before making an extruded map is to deform the basemap to give it a perspective effect. To do this, you can use the deform() function.
We Load an example background map on the US.
See below the original basemap
And here, the deformed basemap.
To better perceive the deformation, you can add a frame with the getframe() function
Notice that you can configure the deformation. The parameter flatten indicates the scaling in Y. The parameter rescale adjust the deformation of the top and bottom of the figure (c(1,1) = no deformation).
basemap <- deform(us, flatten=0.5, rescale = c(0.5, 3))
frame <- deform(getframe(us), flatten=0.5, rescale = c(0.5, 3))
plot(st_geometry(frame), col="#e1e5eb", border = "#99aed1", lwd=0.4)
plot(st_geometry(basemap), col="#99aed1", border = "white", lwd=0.4, add = TRUE)
It can also be useful to rotate the basemap with the rotate() function.
3 - Build an extruded map
Now we’re ready to start extruding… To do this, we use the function - unsustainable suspensions - extrud(). Beware, this is still an experimental function. It has bugs and works very slowly. Sorry for that!
Find below a simple example with the default settings.
Of course, you can custom the map by changing colors (col) and heights (k). On the reprentation below, you will probably notice some imperfections. I hope to be able to fix them in the next version of the package.
4 - Another example with smoothed data
Find here an example of the number of hospital beds in Paris provided in the SpatialPosition package.
library(SpatialPosition)
library(cartography)
pot <- quickStewart(x = hospital,
var = "capacity",
span = 800,
beta = 2, mask = paris,
returnclass = "sf")
breaks <- sort(c(unique(pot$min), max(pot$max)), decreasing = FALSE)
choroLayer(x = pot,
var = "center", breaks = breaks,
legend.pos = "topleft",
legend.title.txt = "Nb. of Beds")
And here’s the same map after extrusion.
library(SpatialPosition)
data("hospital")
# Compute potentials
pot <- quickStewart(x = hospital,
var = "capacity",
span = 800,
beta = 2, mask = paris,
returnclass = "sf")
basemap <- deform(pot)
extrude(basemap, "center", k = 4, col = "white", lwd = 0.7, regular = FALSE)
Be informed that, instead of defining a single color, you can use a field defining the color of each object.
library(SpatialPosition)
data("hospital")
# Compute potentials
pot <- quickStewart(x = hospital,
var = "capacity",
span = 800,
beta = 2, mask = paris,
returnclass = "sf")
pot$colors <- cartography::carto.pal(pal1 = "blue.pal" ,n1 = length(pot$id))
basemap <- deform(pot)
extrude(basemap, "center", k = 4, col = "colors", lwd = 0.2, regular = FALSE)
6 - What about regular grids?
For regular geometries, you can use the parameter regular = TRUE to speed up the calculations. Tip: on irregular basemaps, you can also use this option to test and then remove it to finalize your map.
Here, an example on world population based on on a sample extracted form gpw-v4 data available here
library(raster)
library(rnaturalearth)
# Map Projection
crs <- "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Countries
countries <- ne_countries(scale = 110, type = "countries", returnclass = "sf")
countries <- st_transform(countries, crs)
countries <- countries[countries$adm0_a3 != "ATA",]
countries <- st_sf(st_buffer(st_union(countries),1))
# World population
r <- raster(system.file("worldpop.tif", package="mapextrud"))
dots <- as(r, 'SpatialPointsDataFrame')
dots <- st_as_sf(dots)
dots <- st_transform(dots, crs)
colnames(dots) <- c("pop2020","geometry")
grid <- st_make_grid(countries, cellsize = 300000, square = FALSE)
grid <- aggregate(dots["pop2020"], grid, sum)
# Cartography
basemap <- deform(grid, flatten = 1.2)
extrude(x = basemap, var = "pop2020", k = 4, col = "white", lwd = 0.2, regular = TRUE, add = F)
As everything is based on sf, you can simply combine different layers to finalize the map. Here, we only keep the cells where the population is greater than 5 million people. Thus, 77 % of the world’s population is represented on the map.
library(raster)
library(rnaturalearth)
# Map Projection
crs <- "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Countries
countries <- ne_countries(scale = 110, type = "countries", returnclass = "sf")
countries <- st_transform(countries, crs)
countries <- countries[countries$adm0_a3 != "ATA",]
countries <- st_sf(st_buffer(st_union(countries),1))
# World population
r <- raster(system.file("worldpop.tif", package="mapextrud"))
dots <- as(r, 'SpatialPointsDataFrame')
dots <- st_as_sf(dots)
dots <- st_transform(dots, crs)
colnames(dots) <- c("pop2020","geometry")
grid <- st_make_grid(countries, cellsize = 200000, square = FALSE)
grid <- aggregate(dots["pop2020"], grid, sum)
# Deform
basemap <- deform(grid, flatten = 1.2)
frame <- deform(getframe(countries), flatten = 1.2)
world <- deform(countries, flatten = 1.2)
# Selection
total <- sum(basemap$pop2020, na.rm = TRUE)
basemap <- basemap[basemap$pop2020 > 5000000 & !is.na(basemap$pop2020),]
pct <- round((sum(basemap$pop2020) / total) * 100,0)
# Cartography
plot(st_geometry(frame), col = "#c5dbed", lwd = 0.2)
plot(st_geometry(world), col = "#e3dbca", border = "#7f8aad", lwd = 0.3, add = TRUE)
extrude(x = basemap, var = "pop2020", k = 5, col = "#e65c5c", lwd = 0.2, regular = TRUE, add = TRUE)
# Layout
text(x= -28000000, y =11500000, "An Urban World (population in 2020)", cex = 0.9, col = "#5e698a", pos = 4)
text(x= 42000000, y = -11500000, "Made with the R package mapextrud", cex = 0.6, col = "#414245", pos = 2)
