Making maps with R
Paula Moraga, KAUST
Materials based on the book Geospatial Health Data (2019, Chapman & Hall/CRC)
http://www.paulamoraga.com/book-geospatial
Maps allow us to easily convey spatial information. Here, we show how to create both static and interactive maps by using several mapping packages including ggplot2, leaflet, mapview, and tmap. We create maps of areal data using several functions and parameters of the mapping packages. We also briefly describe how to plot point and raster data. Then, we show how to create maps of flows among locations with the flowmapblue package.
The areal data we map correspond to sudden infant deaths in in the counties of North Carolina, USA, in 1974 and 1979 which are in the sf package.
The path of the data can be obtained with system.file() specifying the directory of the data ("shape/nc.shp") and the package ("sf").
Then, the data can be read with the st_read() function of sf by specifying the path of the data.
We call this data d, and create variables vble and vble2 that we wish to map with the values of the sudden infant deaths in years 1974 and 1979, respectively.
library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
d$vble <- d$SID74
d$vble2 <- d$SID79
head(d)1 ggplot2
The ggplot2 package allows us to create graphics based on the grammar of graphics which defines the rules of structuring mathematic and aesthetic elements to build graphs layer-by-layer.
To create a plot with ggplot2, we call ggplot() specifying arguments data which is a data frame with the variables to plot, and mapping = aes() which are aesthetic mappings between variables in the data and visual properties of the objects in the graph such as position and color of points or lines.
Then, we use + to add layers of graphical components to the graph. Layers consist of geometries, stats, scales, coordinates, facets and themes. For example, we add objects to the graph with geom_*() functions (e.g, geom_point() for points, geom_line() for lines). We can also add color scales (e.g., scale_colour_brewer()), faceting specifications (e.g., facet_wrap() splits data into subsets to create multiple plots), and coordinate systems (e.g., coord_flip()).
We can create maps by using the geom_sf() function and providing a simple feature (sf) object.
For example, we create a map of sudden infant deaths in North Carolina in 1974 (vble) with viridis scale as follows.
library(ggplot2)
library(viridis)
ggplot(d) + geom_sf(aes(fill = vble)) +
scale_fill_viridis() + theme_bw()
Figure 1.1: Map of sudden infant deaths in North Carolina in 1974 created with ggplot2.
Plots created with ggplot2 can be saved with the ggsave() function.
Alternatively, we can specify a device driver (e.g., png, pdf), print the plot, and then shut down the device with dev.off():
png("plot.png")
ggplot(d) + geom_sf(aes(fill = vble)) +
scale_fill_viridis() + theme_bw()
dev.off()Moreover, the plotly package can be used in combination with ggplot2 to create an interactive plot.
Specifically, we can turn a static ggplot object to an interactive plotly object by calling the ggplotly() function of plotly by providing the ggplot object.
library(plotly)
g <- ggplot(d) + geom_sf(aes(fill = vble))
ggplotly(g)The gganimate package can also be used to create animated plots. The syntax of this package is similar to that of ggplot2 and has additional functions to define how data should change along the animation. Information on how to use the gganimate package can be seen at the package’s website.
2 leaflet
The leaflet package makes it easy to create maps using Leaflet which is a very popular open-source JavaScript library for interactive maps.
We can create a leaflet map by calling the leaflet() function passing the spatial object, and adding layers such as polygons and legends using a number of functions.
The sf object that we pass to leaflet() needs to have a geographic coordinate reference system (CRS) indicating latitude and longitude (EPSG code 4326).
Here, we use the st_transform() function of sf to transform the data d which has CRS given by EPSG code 4267 to CRS with EPSG code 4326.
st_crs(d)## Coordinate Reference System:
## User input: NAD27
## wkt:
## GEOGCRS["NAD27",
## DATUM["North American Datum 1927",
## ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4267]]d <- st_transform(d, 4326)We create a color palette using the colorNumeric() function of leaflet specifying the color palette "YlOrRd" of the RColorBrewer package, and the domain with the possible values that can be mapped (d$vble).
Then, we create the map using leaflet() and addTiles() to add a background map to put data into context.
Then, we use addPolygons() to add the polygons representing counties specifying the color of the border (color), interior (fillColor) and opacity (fillOpacity).
Finally, we add a legend with addLegend().
library(leaflet)
pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)
l <- leaflet(d) %>% addTiles() %>%
addPolygons(color = "white", fillColor = ~ pal(vble), fillOpacity = 0.8) %>%
addLegend(pal = pal, values = ~vble, opacity = 0.8)
lFigure 2.1: Map of sudden infant deaths in North Carolina in 1974 created with leaflet.
Note that the default background map added with addTiles() can be changed by another map with addProviderTiles() specifying another tile layer.
Examples of tile layers can be seen at the leaflet providers’ website.
We can also use the addMiniMap() function to add an inset map.
l %>% addMiniMap()The saveWidget() function of htmlwidgets allows us to save the map created to an HTML file.
If we wish to save an image file, we can first save the leaflet map as an HTML file with saveWidget(), and then capture a static version of the HTML using the webshot() function of the webshot package.
The use of webshot requires the installation of the external program PhantomJS which can be installed with webshot::install_phantomjs().
For example, to save the leaflet map l as .png we can proceed as follows:
# Saves map.html
library(htmlwidgets)
saveWidget(widget = l, file = "map.html")
# Takes a screenshot of the map.html created and saves it as map.png
library(webshot)
# webshot::install_phantomjs()
webshot(url = "map.html", file = "map.png")Note we can use getwd() and setwd() to get and set, respectively, the current working directory.
This allows us to see where the files were saved.
Moreover, if in saveWidget() we specify the path where to save map.html as file = "directory/map.html",
the same path needs to be specified in the argument url of webshot() as url = "directory/map.html".
The package webshot2 is meant to be a replacement for webshot.
The webshot2 package uses headless Chrome via the chromote package whereas webshot uses PhantomJS.
3 mapview
The mapview package allows us to very quickly create similar interactive maps as leaflet.
Specifically, we just need to use the mapview() function passing as arguments the spatial object and the variable we want to plot (zcol).
The map created is interactive and by clicking each of the areas we can see popups with the data information.
library(mapview)
mapview(d, zcol = "vble")Maps created with
mapview can also be customized to add elements such as legends and background maps.
For example, we can choose another background map by using the argument map.types, change the color palette with col.regions, and put a title for the legend with layer.name as follows:
library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
mapview(d, zcol = "vble", map.types = "CartoDB.DarkMatter", col.regions = pal, layer.name = "SDI")An inset map can also be added by using the addMiniMap() function of leaflet.
map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)We can save maps created with mapview by using the mapshot() function of mapview
as an HTML file or as a PNG, PDF, or JPEG image.
4 Side by side plots with mapview
We can create a side-by-side plot of maps obtained with mapview separated with a slider by using the | operator. To use the | operator, we need to install the leaflet.extras2 package. When creating the individual maps with mapview(), we use a common legend by specifying the same color palette (col.regions) and breaks (at) so in both maps colors correspond to the same intervals of values.
# common legend
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
at <- seq(min(c(d$vble, d$vble2)), max(c(d$vble, d$vble2)), length.out = 8)
m1 <- mapview(d, zcol = "vble", map.types = "CartoDB.Positron", col.regions = pal, at = at)
m2 <- mapview(d, zcol = "vble2", map.types = "CartoDB.Positron", col.regions = pal, at = at)
m1 | m2