Making maps in R:

Using the sf and tmap Packages

Suggested Citation:

Mendez C. (2020). Making maps in R: Using the sf and tmap Packages. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/tutorial-maps-in-r

This work is licensed under the Creative Commons Attribution-Share Alike 4.0 International License.

1 Libraries

knitr::opts_chunk$set(echo = TRUE, warning=FALSE)

library(tidyverse)  # Modern data science workflow
library(sf)         # Simple features for R
library(tmap)       # Thematic Maps
library(tmaptools)  # Thematic Maps Tools
library(RColorBrewer) # ColorBrewer Palettes
library(leaflet)    # Interactive web maps
library(rgdal)      # Bindings for the Geospatial Data Abstraction Library
library(rgeos)      # Interface to Geometry Engine - Open Source 


# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=4, scipen=999)

2 Tutorial objectives

• Load spatial data files into R
• Join non-spaital data to spatial data files
• Create simple choropleth maps

3 Preliminary material

• Creating thematic maps in R
  • Github repository
• Making maps with R

4 Replication files

• If you are a member of the QuaRCS lab, you can run this tutorial in R Studio Cloud.

4.1 Data

• Non-spatial
• Spatial
• Combined

5 Orginal data sources

The non-spatial datafile is from:

• Öcal, N., & Yildirim, J. (2010). Regional effects of terrorism on economic growth in Turkey: A geographically weighted regression approach. Journal of Peace Research, 47(4), 477-489.

The spatial (shapefile) is from:

• Matthew C. Ingram

6 Import the data

6.1 Non-spatial data

dat <-read.csv("DATASET.csv")

6.1.1 Explore the data

glimpse(dat) 

## Observations: 67
## Variables: 7
## $ id      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 2…
## $ growth  <dbl> 0.0254, -0.0254, 0.0773, 0.2583, 0.0709, 0.0247, 0.1059, 0.03…
## $ base    <dbl> 6.157, 5.900, 5.896, 5.315, 5.923, 6.198, 6.075, 6.173, 6.150…
## $ T       <dbl> 3.585, 2.896, 2.763, 3.140, 2.785, 4.118, 3.296, 3.163, 2.740…
## $ E       <dbl> 5.216, 3.458, 4.995, 3.478, 4.805, 5.492, 4.852, 5.131, 5.146…
## $ G       <dbl> 6.192, 6.048, 5.950, 5.354, 5.940, 6.200, 6.113, 6.219, 6.155…
## $ Coastal <int> 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0…

6.1.1.1 Definitions of variables

• growth: growth rate of per capita real income 1987-2001
• base: lograrithm of the per apita income in the base year 1987
• T: average terrorism index
• E: average years of schooling
• G: real per capita government expenditures in 1987
• Coastal: dummy variable which takes the value of one if the province is a coastal province.

6.2 Spatial data

mapData <- read_sf("MAP.shp")

6.2.1 Explore the data

glimpse(mapData)

## Observations: 67
## Variables: 21
## $ ObjectID   <int> 1946, 1949, 1952, 1957, 2116, 2150, 2186, 2236, 2278, 2371…
## $ NAME       <chr> "Adana", "Adiyaman", "Afyon", "Agri", "Amasya", "Antalya",…
## $ COUNTRY    <chr> "Turkey", "Turkey", "Turkey", "Turkey", "Turkey", "Turkey"…
## $ ISO_CODE   <chr> "TR01", "TR02", "TR03", "TR04", "TR05", "TR07", "TR08", "T…
## $ ISO_CC     <chr> "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR"…
## $ ISO_SUB    <chr> "01", "02", "03", "04", "05", "07", "08", "09", "10", "11"…
## $ ADMINTYPE  <chr> "Province", "Province", "Province", "Province", "Province"…
## $ DISPUTED   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ NOTES      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ AUTONOMOUS <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ COUNTRYAFF <chr> "Turkey", "Turkey", "Turkey", "Turkey", "Turkey", "Turkey"…
## $ CONTINENT  <chr> "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "A…
## $ Land_Type  <chr> "Primary land", "Primary land", "Primary land", "Primary l…
## $ Land_Rank  <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
## $ Shape_Leng <dbl> 8.589, 5.355, 6.698, 6.553, 4.831, 10.889, 4.740, 6.039, 9…
## $ Shape_Area <dbl> 1.6029, 0.8622, 1.4795, 1.1537, 0.6232, 2.1533, 0.7467, 0.…
## $ idnum      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19…
## $ province   <chr> "Adana", "Adiyaman", "Afyon", "Agri", "Amasya", "Antalya",…
## $ id         <dbl> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19…
## $ name_esri  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ geometry   <POLYGON [°]> POLYGON ((36.44 38.22, 36.4..., POLYGON ((39.15 38…

• Check the Coordinate Reference System

st_crs(mapData)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

7 Transform the data

No need to transform any data because both datasets share a common variable id

8 Merge the data

dat_map <- inner_join(
  dat,
  mapData,
  by = "id"
)

8.1 Keep data as sf object

• Keep the data as sf class, so we will not lose the coodinate system

dat_map <- st_as_sf(dat_map)
st_crs(mapData)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

glimpse(dat_map)

## Observations: 67
## Variables: 27
## $ id         <dbl> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19…
## $ growth     <dbl> 0.0254, -0.0254, 0.0773, 0.2583, 0.0709, 0.0247, 0.1059, 0…
## $ base       <dbl> 6.157, 5.900, 5.896, 5.315, 5.923, 6.198, 6.075, 6.173, 6.…
## $ T          <dbl> 3.585, 2.896, 2.763, 3.140, 2.785, 4.118, 3.296, 3.163, 2.…
## $ E          <dbl> 5.216, 3.458, 4.995, 3.478, 4.805, 5.492, 4.852, 5.131, 5.…
## $ G          <dbl> 6.192, 6.048, 5.950, 5.354, 5.940, 6.200, 6.113, 6.219, 6.…
## $ Coastal    <int> 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0…
## $ ObjectID   <int> 1946, 1949, 1952, 1957, 2116, 2150, 2186, 2236, 2278, 2371…
## $ NAME       <chr> "Adana", "Adiyaman", "Afyon", "Agri", "Amasya", "Antalya",…
## $ COUNTRY    <chr> "Turkey", "Turkey", "Turkey", "Turkey", "Turkey", "Turkey"…
## $ ISO_CODE   <chr> "TR01", "TR02", "TR03", "TR04", "TR05", "TR07", "TR08", "T…
## $ ISO_CC     <chr> "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR", "TR"…
## $ ISO_SUB    <chr> "01", "02", "03", "04", "05", "07", "08", "09", "10", "11"…
## $ ADMINTYPE  <chr> "Province", "Province", "Province", "Province", "Province"…
## $ DISPUTED   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ NOTES      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ AUTONOMOUS <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ COUNTRYAFF <chr> "Turkey", "Turkey", "Turkey", "Turkey", "Turkey", "Turkey"…
## $ CONTINENT  <chr> "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "A…
## $ Land_Type  <chr> "Primary land", "Primary land", "Primary land", "Primary l…
## $ Land_Rank  <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
## $ Shape_Leng <dbl> 8.589, 5.355, 6.698, 6.553, 4.831, 10.889, 4.740, 6.039, 9…
## $ Shape_Area <dbl> 1.6029, 0.8622, 1.4795, 1.1537, 0.6232, 2.1533, 0.7467, 0.…
## $ idnum      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19…
## $ province   <chr> "Adana", "Adiyaman", "Afyon", "Agri", "Amasya", "Antalya",…
## $ name_esri  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ geometry   <POLYGON [°]> POLYGON ((36.44 38.22, 36.4..., POLYGON ((39.15 38…

9 Plot Thematic Maps

9.1 Quick Map

• A quick map of the Terrorism variable

qtm(dat_map, fill = "T")

9.2 Simple Map

tm_shape(dat_map) + tm_fill("T") 

9.3 Remove frame

tm_shape(dat_map) + 
  tm_fill("T") +
  tm_layout(frame = FALSE)

9.4 Set legend outside

tm_shape(dat_map) + tm_fill("T") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

9.5 Set legend size

tm_shape(dat_map) + tm_fill("T") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_legend(text.size = 0.3) 

9.6 Set color palette

See color reference here

tm_shape(dat_map) + 
  tm_fill("T", palette = "Greens") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

tm_shape(dat_map) + 
  tm_fill("T", palette = "-Greens") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

9.7 Add polygon names

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_text("province", size = "Shape_Area",  auto.placement = F, legend.size.show = FALSE) 

Documentation

9.8 Add scale bar

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_text("province", size = "Shape_Area",  auto.placement = F, legend.size.show = FALSE)  + 
  tm_scale_bar(position = c("RIGHT", "BOTTOM"))

9.9 Set color intervals

Enter “style =” followed by one of the options below.

• equal: divides the range of the variable into n parts.
• pretty: chooses a number of breaks to fit a sequence of equality spaced ‘round’ values.
• quantile: equal number of cases in each group
• jenks: looks for natural breaks in the data
• cat: when the variable is categorical

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "Reds"
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

Change the number of intervals in the color scheme and how the intervals are spaced. Changing the number of intervals n = 7. So, we have 7 shades instead of the default 5.

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          n = 7,
          palette = "Reds"
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

9.10 Add histogram

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          n = 5,
          palette = "YlOrBr",
          legend.hist = TRUE
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

9.11 Add borders

You can edit the borders of the shapefile with the tm_borders() function which has many arguments. alpha denotes the level of transparency on a scale from 0 to 1 where 0 is completely transparent.

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "Blues"
    ) + 
  tm_borders(alpha=.4) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

9.12 Add compass

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "viridis"
    ) + 
  tm_borders(alpha=.4) + 
  tm_compass(type = "8star", 
              position = c("RIGHT", "BOTTOM"),
              show.labels = 2,
              text.size = 0.35) +
  tm_layout(legend.outside = TRUE,
            frame = FALSE)

Documentation

9.13 Edit the layout

tm_shape(dat_map) + 
  tm_fill("T",
          palette = "viridis",
          style = "quantile",
          title = "Terrorist Attacks"
          ) + 
tm_borders(alpha=.4) + 
tm_layout(
  legend.text.size = 0.7,
  legend.title.size = 1,
  legend.position = c("right", "bottom"),
  legend.outside = TRUE,
  frame = FALSE
  ) 

Documentation

9.14 Interactive map

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(dat_map) + 
  tm_fill("T",
          palette = "viridis",
          id="province",
          popup.vars=c("T", "growth", "base", "E", "G")
          ) + 
    tm_legend(outside=TRUE) +
  tm_layout(frame = FALSE) 

tmap_mode("plot")

## tmap mode set to plotting

Documentation

10 Save a new shapefile

st_write(dat_map, "dat_map.shp")

11 Other tutorials

• ECS530: (IV) Visualization of spatial data

• Practical 5: Making maps in R

• Making maps with R

• Chapter 8 Solutions

12 References

• https://rdrr.io/cran/tmap/

• Tennekes, M., 2018, tmap: Thematic Maps in R, Journal of Statistical Software, 84(6), 1-39

• https://marcinstepniak.eu/post/interactive-choropleth-maps-with-r-and-tmap-part-i/

• Modern Geospatial Data Analysis with R

• Webinar on Spatial Data Analysis with R

• https://tlorusso.github.io/geodata_workshop/tmap_package

• https://rstudio-pubs-static.s3.amazonaws.com/229482_b5b182f650434a36ace4eeef2d2bcd83.html#17

• https://mran.microsoft.com/snapshot/2017-12-11/web/packages/tmap/vignettes/tmap-nutshell.html

• Chapter 8 Solutions

• Geocomputation with R website

• https://bookdown.org/lexcomber/brunsdoncomber2e/

END