Introduction
The article Tennekes, M. tmap: Thematic Maps in R (2018), Journal
of Statistical Software 84(6), 1-39, http://dx.doi.org/10.18637/jss.v084.i06 describes
version 1.x of tmap and tmaptools. Since tmap and tmaptools version 2.0
are based on sf objects instead of sp objects,
there are some changes in the code, especially regarding the necessary
processing step. The code for the plotting methods is fully backward
compatible, although some messages or warnings may be given for
deprecated functions or arguments.
This code replaces the script v84i06.R for tmap and tmaptools version 2.0. The produced figures are identical to those included in the article.
See vignette("tmap-changes")
for changes in version 2.x/3.x, which is recommended for people who have
read the JSS article. For people who are new to tmap, see vignette("tmap-getstarted").
###########################################################################
## This script will replicate the figures (except the screenshots)
## and write them to files.
## The working directory should be set the the parent folder of this script.
###########################################################################
## install.packages(c("tmap", "tmaptools"))
library("tmap") # required version 2.0 or later
library("tmaptools") # required version 2.0 or later
data("World", "metro", package = "tmap")
metro$growth <- (metro$pop2020 - metro$pop2010) / (metro$pop2010 * 10) * 100Section 1.1
#############################
## Figure 1
#############################
m1 <- tm_shape(World) +
tm_polygons("income_grp", palette = "-Blues",
title = "Income class", contrast = 0.7, border.col = "grey30", id = "name") +
tm_text("iso_a3", size = "AREA", col = "grey30", root = 3) +
tm_shape(metro) +
tm_bubbles("pop2010", col = "growth", border.col = "black",
border.alpha = 0.5,
breaks = c(-Inf, 0, 2, 4, 6, Inf) ,
palette = "-RdYlGn",
title.size = "Metro population (2010)",
title.col = "Annual growth rate (%)",
id = "name",
popup.vars = c("pop2010", "pop2020", "growth")) +
tm_style("gray") +
tm_format("World", frame.lwd = 2)
tmap_save(m1, "bubble.png", width = 6.125, height = 3, scale = .75, dpi = 300, asp = 0, outer.margins = 0)Section 3.1
#############################
## Figure 2
#############################
m0 <- tm_shape(metro) +
tm_bubbles(size = "pop2030") +
tm_style("cobalt") +
tm_format("World")
tmap_save(m0, "metro2030.png", width = 6.125, scale = .5, dpi = 300, outer.margins = 0)Section 3.2
#############################
## Figure 3
#############################
m21 <- tm_shape(World) + tm_polygons(c("blue", "red")) + tm_layout(frame.lwd = 1.5)
tmap_save(m21, "facets1.png", width = 6.125, height = 1.54, scale = .75, dpi = 300, outer.margins = 0)
#############################
## Figure 4
#############################
m22 <- tm_shape(World) + tm_polygons("red") + tm_facets(by = "continent", free.coords = FALSE)
tmap_save(m22, "facets2.png", width = 6.125, height = 1.8, scale = .75, dpi = 300, outer.margins = 0)
#############################
## Figure 5
#############################
tm1 <- tm_shape(World) + tm_polygons()
tm2 <- tm_shape(metro) + tm_dots()
png("facets3.png", width = 6.125, height = 1.54, units = "in", res = 300)
tmap_arrange(tm1, tm2, outer.margins = .01)
dev.off()Section 4.1
#############################
## Figure 6
#############################
tmap_mode("view")
m1Section 4.2
#############################
## Figure 7
#############################
data("land", "rivers", package = "tmap")
m2 <- tm_shape(land) +
tm_raster("elevation", breaks = c(-Inf, 250, 500, 1000, 1500, 2000, 2500, 3000, 4000, Inf),
palette = terrain.colors(9), title = "Elevation (m)") +
tm_shape(rivers) +
tm_lines("lightblue", lwd = "strokelwd", scale = 1.5, legend.lwd.show = FALSE) +
tm_shape(World, is.master = TRUE) +
tm_borders("grey20", lwd = .5) +
tm_grid(projection = "longlat", labels.size = 0.4, lwd = 0.25) +
tm_text("name", size = "AREA") +
tm_compass(position = c(0.08, 0.45), color.light = "grey90", size = 3) +
tm_credits("Eckert IV projection", position = c("RIGHT", "BOTTOM")) +
tm_style("classic",
bg.color = "lightblue",
space.color = "grey90",
inner.margins = c(0.04, 0.04, 0.03, 0.02),
earth.boundary = TRUE) +
tm_legend(position = c("left", "bottom"),
frame = TRUE,
bg.color = "lightblue")
tmap_save(m2, "classic.png", width = 6.125, scale = .7, dpi = 300, outer.margins = 0)Section 4.3
#############################
## Figure 8
#############################
m3 <- tm_shape(World, projection = "robin") +
tm_polygons(c("HPI", "gdp_cap_est"),
palette = list("RdYlGn", "Purples"),
style = c("pretty", "fixed"), n = 7,
breaks = list(NULL, c(0, 500, 2000, 5000, 10000, 25000, 50000, Inf)),
title = c("Happy Planet Index", "GDP per capita")) +
tm_style("natural", earth.boundary = c(-180, -87, 180, 87)) +
tm_format("World", inner.margins = 0.02, frame = FALSE) +
tm_legend(position = c("left", "bottom"), bg.color = "gray95", frame = TRUE) +
tm_credits(c("", "Robinson projection"), position = c("RIGHT", "BOTTOM"))
tmap_save(m3, "world_facets2.png", width = 5, scale = .7, dpi = 300, outer.margins = 0)Section 4.4
#############################
## Figure 9
#############################
library("readxl")
library("grid")
# function to obtain Food Environment Atlas data (2014)
get_food_envir_data <- function() {
dir <- tempdir()
if (!file.exists(file.path(dir, "DataDownload.xls"))) {
download.file("https://www.ers.usda.gov/webdocs/DataFiles/48731/February2014.xls?v=41688", destfile = file.path(dir, "DataDownload.xls"), mode = "wb")
}
res <- tryCatch({
read_excel(file.path(dir, "DataDownload.xls"), sheet = "HEALTH")
}, error = function(e) {
stop("The excel file cannot be read. Please open it, and remove all sheets except HEALTH. The location of the file is: ", normalizePath(file.path(dir, "DataDownload.xls")))
})
}
# function to obtain US county shape
get_US_county_2010_shape <- function() {
dir <- tempdir()
download.file("http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip", destfile = file.path(dir, "gz_2010_us_050_00_20m.zip"))
unzip(file.path(dir, "gz_2010_us_050_00_20m.zip"), exdir = dir)
US <- sf::read_sf(file.path(dir, "gz_2010_us_050_00_20m.shp"))
levels(US$NAME) <- iconv(levels(US$NAME), from = "latin1", to = "utf8")
US
}
# obtain Food Environment Atlas data
FEA <- get_food_envir_data()
# obtain US county shape
US <- get_US_county_2010_shape()
us1 <- qtm(US)
tmap_save(us1, "US1.png", scale = .5, width = 6.125, asp = 0, outer.margins = 0)
#############################
## Figure 10
#############################
library(dplyr)
library(sf)
US$FIPS <- paste0(US$STATE, US$COUNTY)
# append data to shape
#US <- append_data(US, FEA, key.shp = "FIPS", key.data = "FIPS", ignore.duplicates = TRUE)
US <- left_join(US, FEA, by = c("FIPS", "FIPS"))
unmatched_data <- FEA %>% filter(!(FIPS %in% US$FIPS))
tmap_mode("view")
qtm(US, fill = "PCT_OBESE_ADULTS10")
#############################
## Figure 11
#############################
US_cont <- US %>%
subset(!STATE %in% c("02", "15", "72")) %>%
simplify_shape(0.2)
US_AK <- US %>%
subset(STATE == "02") %>%
simplify_shape(0.2)
US_HI <- US %>%
subset(STATE == "15") %>%
simplify_shape(0.2)
# create state boundaries
US_states <- US_cont %>%
dplyr::select(geometry) %>%
aggregate(by = list(US_cont$STATE), FUN = mean)
# contiguous US
m_cont <- tm_shape(US_cont, projection = 2163) +
tm_polygons("PCT_OBESE_ADULTS10", border.col = "gray50", border.alpha = .5, title = "", showNA = TRUE) +
tm_shape(US_states) +
tm_borders(lwd = 1, col = "black", alpha = .5) +
tm_credits("Data @ Unites States Department of Agriculture\nShape @ Unites States Census Bureau", position = c("right", "bottom")) +
tm_layout(title = "2010 Adult Obesity by County, percent",
title.position = c("center", "top"),
legend.position = c("right", "bottom"),
frame = FALSE,
inner.margins = c(0.1, 0.1, 0.05, 0.05))
# Alaska inset
m_AK <- tm_shape(US_AK, projection = 3338) +
tm_polygons("PCT_OBESE_ADULTS10", border.col = "gray50", border.alpha = .5, breaks = seq(10, 50, by = 5)) +
tm_layout("Alaska", legend.show = FALSE, bg.color = NA, title.size = 0.8, frame = FALSE)
# Hawaii inset
m_HI <- tm_shape(US_HI, projection = 3759) +
tm_polygons("PCT_OBESE_ADULTS10", border.col = "gray50", border.alpha = .5, breaks = seq(10, 50, by = 5)) +
tm_layout("Hawaii", legend.show = FALSE, bg.color = NA, title.position = c("LEFT", "BOTTOM"), title.size = 0.8, frame = FALSE)
# specify viewports for Alaska and Hawaii
vp_AK <- viewport(x = 0.15, y = 0.15, width = 0.3, height = 0.3)
vp_HI <- viewport(x = 0.4, y = 0.1, width = 0.2, height = 0.1)
# save map
tmap_mode("plot")
tmap_save(m_cont, "USchoro.png", scale = 0.7, width = 6.125, outer.margins = 0,
insets_tm = list(m_AK, m_HI),
insets_vp = list(vp_AK, vp_HI))Section 4.5
#############################
## Figure 12a
#############################
library("sf")
library("rnaturalearth")
# functions to obtain crimes data
get_crimes_data <- function(path) {
stopifnot(file.exists(path), ("crimes_in_Greater_London_2015-10.zip" %in% list.files(path)))
tmp_dir <- tempdir()
unzip(file.path(path, "crimes_in_Greater_London_2015-10.zip"), exdir = tmp_dir)
rbind(read.csv(file.path(tmp_dir, "2015-10-city-of-london-street.csv")),
read.csv(file.path(tmp_dir, "2015-10-metropolitan-street.csv")))
}
# please download the file "crimes_in_Greater_London_2015-10.zip" (available on https://www.jstatsoft.org as a supplement of this paper), and change the path argument below to the location of the downloaded file:
crimes <- get_crimes_data(path = "./")
# create sf of known locations
crimes <- crimes[!is.na(crimes$Longitude) & !is.na(crimes$Latitude), ]
crimes <- st_as_sf(crimes, coords = c("Longitude", "Latitude"), crs = 4326)
# set map projection to British National Grid
crimes <- st_transform(crimes, crs = 27700)
c1 <- qtm(crimes)
tmap_save(c1, "crimes1.png", scale = .6, width = 3, units = "in", outer.margins = 0)
#############################
## Figure 12b
#############################
crimes_osm <- read_osm(crimes)
c2 <- qtm(crimes_osm) + qtm(crimes, symbols.col = "red", symbols.size = 0.5)
tmap_save(c2, "crimes2.jpg", scale = .6, width = 3, units = "in", outer.margins = 0)
#############################
## Figure 13
#############################
c3 <- qtm(crimes_osm, raster.saturation = 0, raster.alpha = 1) +
qtm(crimes, symbols.col = "Crime.type", symbols.size = 0.5) +
tm_legend(outside = TRUE)
tmap_save(c3, "crimes3.png", scale = .8, width = 5, height = 4, units = "in", outer.margins = 0)
#############################
## Figure 14
#############################
regions <- ne_download(scale = "large", type = "states", category = "cultural", returnclass = "sf")
london <- regions[which(regions$region == "Greater London"),]
london <- st_transform(london, crs = 27700)
# remove crimes outside Greater London
crimes_london <- crop_shape(crimes, london, polygon = TRUE)
c3b <- qtm(crimes_london, dots.alpha = 0.1) +
tm_shape(london) +
tm_borders()
tmap_save(c3b, "crimes3b.png", scale = .7, width = 6.125, units = "in", outer.margins = 0)
#############################
## Figure 15
#############################
library("devtools")
install_github("mtennekes/oldtmaptools")
library(oldtmaptools)
crime_densities <- smooth_map(crimes_london, bandwidth = 0.5, breaks = c(0, 50, 100, 250, 500, 1000), cover = london)
# download rivers, and get Thames shape
rivers <- ne_download(scale = "large", type = "rivers_lake_centerlines", category = "physical")
thames <- crop_shape(rivers, london)
c4 <- tm_shape(crime_densities$polygons) +
tm_fill(col = "level", palette = "YlOrRd", title = expression("Crimes per " * km^2)) +
tm_shape(london) + tm_borders() +
tm_shape(thames) + tm_lines(col = "steelblue", lwd = 4) +
tm_compass(position = c("left", "bottom")) +
tm_scale_bar(position = c("left", "bottom")) +
tm_style("gray", title = "Crimes in Greater London\nOctober 2015")
tmap_save(c4, "crimes4.png", scale = .7, width = 6.125, units = "in", outer.margins = 0)
#############################
## Figure 16
#############################
london_city <- london[london$name == "City",]
crimes_city <- crop_shape(crimes_london, london_city, polygon = TRUE)
london_osm <- read_osm(london_city, type = "stamen-watercolor", zoom = 13)
c5 <- qtm(london_osm) +
tm_shape(crimes_city) +
tm_dots(size = .2) +
tm_facets("Crime.type", free.coords = FALSE)
tmap_save(c5, "crimes5.png", scale = 1, width = 6.125, asp = NA, outer.margins = 0)
#############################
## Figure 17
#############################
crime_lookup <- c("Anti-social behaviour" = 2,
"Bicycle theft" = 1,
"Burglary" = 1,
"Criminal damage and arson" = 2,
"Drugs" = 6,
"Other crime" = 7,
"Other theft" = 1,
"Possession of weapons" = 3,
"Public order" = 2,
"Robbery" = 1,
"Shoplifting" = 1,
"Theft from the person" = 1,
"Vehicle crime" = 4,
"Violence and sexual offences" = 5)
crime_categories <- c("Property Crime",
"Criminal damage and anti-social behaviour",
"Possession of weapons",
"Vehicle crime",
"Violence and sexual offences",
"Drugs",
"Other crime")
crimes_city$Crime.group <- factor(crime_lookup[crimes_city$Crime.type], labels = crime_categories)
tmap_mode("view")
tm_shape(crimes_city) +
tm_dots(jitter = 0.2, col = "Crime.group", palette = "Dark2", popup.vars = TRUE) +
tm_view(alpha = 1,
basemaps = "Esri.WorldTopoMap")