# reproduction code from JSS article ‘tmap: Thematic Maps in R’

### 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) * 100

### Section 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")
    m1

### Section 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")