Making Maps with R

Richard A. Lent

Maps, R Notebooks, and Code Chunks

Making maps on a computer has traditionally required the use of desktop geographic information system (GIS) software such as ArcGIS or QGIS. An alternative is to use R, a free software environment for statistical computing and graphics. R has many features that allow it to read GIS data and produce both static and interactive maps. This document (which is an R Notebook) shows how to make maps with R and RStudio, using R base graphics and the maps and mapdata packages, in addition to the leaflet and tmap packages. We will also map U.S. Census data with help from the tigris and acs packages.

It is assumed that the reader has some familiarity with R and RStudio. For additional information see An Introduction to R and RStudio as a Research and Writing Platform. See also Why I love R Notebooks. Also please see Habitat structure and phenotypic variation in the invading butterfly Coenonympha tullia, an R Notebook that illustrates how to use R and RStudio to write a complete academic manuscript.

To install R, go to the R download site. You must install R before installing RStudio. After installing R, visit the RStudio download site to download the RStudio installer for your computer platform. Choose the RStudio Desktop Open Source License (where it says FREE). When you run the RStudio installer, answer “yes” if it asks to install the command-line tools. Some R packages require the command-line tools so that they can be compiled and installed.

This R Notebook may be viewed online, where you can download the R Notebook document (click the first button labeled Code, then click Download Rmd). The notebook is a plain text file named maps.Rmd, which can be opened in RStudio. The notebook contains text written using R Markdown, a simple formatting language, interspersed with blocks of R source code, called “chunks.” The code chunks can be executed independently and interactively from inside the notebook using controls contained in each chunk. Output is visible immediately beneath the code chunk. Within a code chunk, the # symbol is used to denote comments, which are ignored by the R interpreter and serve to document the code.

A code chunk begins and ends with three backticks (found to the left of the numeral 1 on your keyboard). Following the initial three backticks is a pair of curly braces containing the letter r to signify that it is an R code chunk. The curly braces can contain additional R options, and all R code must be contained inside of the backticks. You can use the Insert menu of the RStudio text editor to insert an empty R code chunk ready for coding.

The following chunk establishes global options that apply to all of the other chunks in the document. Each chunk has a short title, following the r after the first curly brace, making it easier to navigate the chunks in the RStudio text editor.

knitr::opts_chunk$set(
    echo = FALSE,
    message = FALSE,
    warning = FALSE,
    comment = NA
)

Next we install a series of R packages, and other packages that those packages require, which are needed for the examples. The code first makes a list of the packages we want to install, then checks if they are already installed. The packages and their dependencies are then installed and loaded into the R environment. If you are running this for the first time after you have installed R and RStudio, be patient, as the packages may take a while to download and install.

# List of required packages.
.packages = c("maps", "mapdata", "Hmisc", "leaflet", "htmlwidgets", "rgdal", "tmap", "tigris", "acs", "stringr", "stringi", "sp", "dplyr", "mapview", "knitr", "rmarkdown")
# Install CRAN packages if not already installed.
.inst <- .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])
# Load packages without printing messages. 
suppressMessages(invisible(lapply(.packages, require, character.only=TRUE)))

In the various code chunks, we use the R functions invisible and suppressMessages to eliminate unwanted console output and automatic printing of values returned from functions. This is handy when the return value is an intermediate map layer that we want to use later but do not want to be printed to the screen or to the R Notebook. We also use function arguments, such as verbose = FALSE, to suppress additional console output that would otherwise clutter the R Notebook.

Each of the following four R code chunks produces a map using different technologies in R. The code can be executed interactively from inside the R Notebook. The code can also be copy-pasted directly into the R console and run from there, or copy-pasted into an empty R script file and run as a standalone program. The first comment line of each code chunk contains a suggested name for the script file.

Map 1: Base R graphics with the maps and mapdata packages

Our first map is created using the base R graphics package along with the maps and mapdata packages. Base R graphics is the standard graphics package that comes with R. It works by writing a series of graphics commands to a graphics window. The maps package enables display of geographic maps using map databases contained in the mapdata package.

We start by reading an external data file called sites.csv into an R data frame, a spreadsheet-like structure in which rows represent individual items of interest and columns are variables measured for each item. sites.csv is a real dataset and is described in Habitat structure and phenotypic variation in the invading butterfly Coenonympha tullia. The maps package is then used to draw the boundaries of Vermont, New Hampshire, and Massachusetts. Then, the sites data frame is used to map the locations of 11 study sites using their latitude-longitude coordinates, and to plot site labels next to each location. Next, the base graphics text function is used to place text labels for each of the three states. Finally, a box is drawn around the map and a scale added.

The code chunk also performs a simple linear regression of some butterfly morphological variables (see Habitat structure and phenotypic variation in the invading butterfly Coenonympha tullia) and saves the regression residuals to the sites data frame. There is an optional line of R code that uses the residuals to control the size of the plotting symbol. This illustrates the great advantage of using R to make maps, in that it is easy (if you know R) to map all manner of statistical results.

The appearance of the map can be changed by altering the parameters passed into the various graphics statements. For example, the col parameter controls the color used to draw a particular feature, pch changes the plotting symbol, and cex alters the size of a plotting symbol. At the R console, you can enter show.col and show.pch to display available colors and plotting symbols, respectively. Also, the colors() function will list the names of 657 available colors that can be used with the col parameter.

# sitesmap1.R
# Map the sites data using the maps package.
# Use the maps package to plot state boundaries 
# contained in the mapdata package
# and then overlay site locations using the points 
# function. Additional text is plotted
# using the text function to locate text using the 
# latitude-longitude coordinates of the map.
# We then draw a box around the map using the box 
# function and add a scale with the map.scale function.
# All of this works because the base graphics package in R 
# will keep drawing a sequence of graphic statements to 
# the same graphics window until the window is closed 
# or a new graphics window is created.
# The Hmisc package has functions show.col and show.pch to 
# display available color codes and plotting symbols, respectively.
# Use colors() to list names of 657 additional colors.
# See also the R color cheatsheet:
# https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf
library(maps)     # Draw geographical maps.
library(mapdata)  # Map databases.
library(Hmisc)    # Miscellaneous useful functions.
sites <- read.csv("http://college.holycross.edu/faculty/rlent/sites/sites.csv")
# Do a simple linear regression of butterfly forewing length vs thorax width, 
# and save the regression residuals to the sites data frame.
reg <- with(sites, lm(fwlength ~ thorax))
sites$resid <- reg$residuals
# Make the map.
map('state', c('Vermont','New Hampshire','Massachusetts'), 
    xlim=c(-74, -69), ylim=c(41, 45.5),  
    col='gray90', fill=TRUE, mar = c(0, 0, 1, 0))
with(sites, points(lon_dd, lat_dd, pch=16, col='black', cex=.7))

# Comment out the previous 'with' statement, and uncomment the following 
# 'with' statement, to size the site points by the magnitude of their regression residual.
# with(sites, points(lon_dd, lat_dd, pch=16, col='black', cex=abs(resid)+1))
with(sites, text(lon_dd-.15, lat_dd, cex=0.7, labels=code))
text(-72, 42.4, labels='MA', cex=1.5)

text(-71.5, 43.4, labels='NH', cex=1.5)
text(-72.9, 44.3, labels='VT', cex=1.5)

box()
map.scale(ratio = FALSE, metric = FALSE)

Map 2: The leaflet package

We next produce the same basic map as in the previous example but using the leaflet package, which produces a prettier, interactive map that can be saved as a web page, exported to an image file, or embedded in an R Markdown document. Leaflet is the leading open-source JavaScript library for mobile-friendly interactive maps. (See this example from the New York Times.) The Leaflet package is an R wrapper for the Leaflet JavaScript mapping library, making it easy to generate interactive maps based on spatial data you have in R.

The function leaflet() returns a Leaflet map widget object that can receive various map layers, such as base maps and point, line, and polygon features. (A list of free base map providers for Leaflet is here.) Leaflet can map spatial data from several sources in R, including data frames containing latitudes and longitudes in addition to maps from the maps package. An additional package called htmlwidgets is used to save the Leaflet widget to an HTML page. The mapshot function of the mapview package can save a Leaflet interactive map as a static image file, for inclusion in documents, presentations, etc. There is an optional line of code at the end of this chunk to perform this function.

# sitesmap2.R
# Map the sites data using the leaflet package.
# This creates a prettier, interactive map.
library(leaflet)
library(maps)
library(htmlwidgets) # To save the map as a web page.
# The data to map.
sites <- read.csv("http://college.holycross.edu/faculty/rlent/sites/sites.csv")
# State boundaries from the maps package. The fill option must be TRUE.
bounds <- map('state', c('Massachusetts', 'Vermont', 'New Hampshire'), fill=TRUE, plot=FALSE)
# A custom icon.
icons <- awesomeIcons(
    icon = 'disc',
    iconColor = 'black',
    library = 'ion', # Options are 'glyphicon', 'fa', 'ion'.
    markerColor = 'blue',
    squareMarker = TRUE
)
# Create the Leaflet map widget and add some map layers.
# We use the pipe operator %>% to streamline the adding of
# layers to the leaflet object. The pipe operator comes from 
# the magrittr package via the dplyr package.
map <- leaflet(data = sites) %>%
    # setView(-72.14600, 43.82977, zoom = 8) %>% 
    addProviderTiles("CartoDB.Positron", group = "Map") %>%
    addProviderTiles("Esri.WorldImagery", group = "Satellite") %>% 
    addProviderTiles("Esri.WorldShadedRelief", group = "Relief") %>%
    # Marker data are from the sites data frame. We need the ~ symbols
    # to indicate the columns of the data frame.
    addMarkers(~lon_dd, ~lat_dd, label = ~locality, group = "Sites") %>% 
    # addAwesomeMarkers(~lon_dd, ~lat_dd, label = ~locality, group = "Sites", icon=icons) %>%
    addPolygons(data=bounds, group="States", weight=2, fillOpacity = 0) %>%
    addScaleBar(position = "bottomleft") %>%
    addLayersControl(
        baseGroups = c("Map", "Satellite", "Relief"),
        overlayGroups = c("Sites", "States"),
        options = layersControlOptions(collapsed = FALSE)
    )
invisible(print(map))

# Save the interactive map to an HTML page.
saveWidget(map, file="sitesmap2.html", selfcontained=TRUE)
# Uncomment the following line to save the interactive map to a static image file.
# mapshot(map, file="sitesmap2.png")

Map 3: The tmap package

The tmap package “offers a flexible, layer-based, and easy to use approach to create thematic maps, such as choropleths and bubble maps.” See tmap in a nutshell for further details. Here we present a simple example of using tmap to create a thematic map from an ESRI shapefile. The shapefile is a common GIS vector data format that has been around for a long time, and is used to represent points, lines, and enclosed areas (polygons). A shapefile is actually composed of several files that function together to provide the geographic features (points, lines, or polygons) along with a data table containing attributes of each feature.

There are many online repositories of GIS data that provide downloadable map layers in shapefile format. Here we use a shapefile from the Massachusetts Office of Geographic Information (MassGIS) that contains polygons delineating areas of critical enviromental concern. The R package rgdal provides function readOGR to read a shapefile into a SpatialPolygonsDataFrame, which is an R data frame that is able to contain spatial data. We then use the qtm (“quick thematic map”) function of the tmap package to produce the map. The polygons are color-coded according to their names, which are contained in the field NAME from the attribute table of the original shapefile downloaded from MassGIS. The code also uses two functions from the R utils package. The download.file function downloads the shapefile from MassGIS, using a URL passed to it. The multiple files comprising the shapefile map layer are provided by MassGIS as a compressed zip file. We then use the unzip function to extract the shapefile components to our working directory. The shapefile layer, named acecs_poly, is then passed to the readOGR function so that it can be brought into R.

tmap can function in two modes: plot mode and view mode. Plot mode creates a non-interactive map using R graphical plotting, whereas view mode creates an interactive Leaflet map. The tmap_mode function can be used to switch between plot mode and view mode.

# shape.R
# Use the tmap package to make a thematic map from an ESRI shapefile.
# The shapefile is from MassGIS, Areas of Critical Environmental Concern. 
# See https://goo.gl/vQxakO.
library(rgdal)
library(tmap)
# Download and extract the data.
# The polygon shapefile will be extracted to the working directory (use getwd()).
url <- "http://wsgw.mass.gov/data/gispub/shape/state/acecs.zip"
download.file(url, destfile = "acecs.zip", quiet = TRUE)
unzip("acecs.zip")
suppressMessages(tmap_mode("view")) # Options are "plot", "view".
map <- readOGR("./","acecs_poly", stringsAsFactors = FALSE, verbose = FALSE)
invisible(print(qtm(map, "NAME", fill.title="Name", format = "Europe_wide")))

Map 4: Mapping U. S. Census Data with the acs, tigris, and leaflet packages

Our final example is a bit more complicated. We use the leaflet package to map some data downloaded from the U. S. Census Bureau, using the acs and tigris packages. The acs package downloads tabular census data, while the tigris package is used to load Census TIGER shapefiles into R. TIGER stands for Topologically Integrated Geographic Encoding and Referencing, and TIGER shapefiles contain features such as roads, railroads, rivers, as well as legal and statistical geographic areas used by the U. S. Census.

This example was compiled and modified from Manipulating and mapping US Census data in R using the acs, tigris and leaflet packages, using the easy(er) way described in that document. Basically, the census shapfiles are downloaded using tigris and read into R spatial data frames. The tablular census data are downloaded using acs and the merged into the shapefile’s attribute table. From there, the merged data can then be given to the leaflet function to produce the interactive map.

While appreciably more complex than our previous examples, this map serves to illustrate the degree of data access, processing, and spatial mapping that can be accomplished with R and RStudio.

# census.R
# Mapping U. S. Census Data with R using the acs, tigris and leaflet packages.
# R code compiled and modified from https://goo.gl/TvLCRo
# Load required R packages.
library(tigris)      # U. S. Census maps.
library(acs)         # U. S. Census data.
library(stringr)     # To pad fips codes.
library(rgdal)       # For readOGR and others.
library(sp)          # For spatial objects.
library(leaflet)     # For interactive maps.
library(dplyr)       # For working with data frames.
library(htmlwidgets) # For saving the map.
# Download and store XML tables used to look up variable codes, table names, 
# and other metadata associated with the acs package. (May not be needed.)
# acs.tables.install() 
# Download the spatial data (tigris).
# Using Federal Information Processing Standards (FIPS) codes for counties.
# You can view this layer at the R console with plot(tracts). 
counties <- c(5, 47, 61, 81, 85)
suppressMessages(tracts <- tracts(state = 'NY', county=c(5, 47, 61, 81, 85), cb=TRUE))
# We need an API key (Application Program Interface) to access U. S. Census data.
# You can get your own Census data API key from 
# http://api.census.gov/data/key_signup.html.
api.key.install(key="474205a8b83bc35e69282fe60f5ea05170a06161")
# Create a geographic set to grab tabular data (acs).
geo <- geo.make(state=c("NY"),
              county=c(5, 47, 61, 81, 85), tract="*")
# Retrieve data for the five year span ending in 2012.
# You can view this data table by typing 'income' at the R console.
# Use of col.names = "pretty" gives the full column names.
# If instead you want Census variable IDs use col.names="auto". 
income <- acs.fetch(endyear = 2012, span = 5, geography = geo,
                  table.number = "B19001", col.names = "pretty")
# Convert to a data.frame for merging
income_df <- data.frame(paste0(str_pad(income@geography$state, 2, "left", pad="0"), 
          str_pad(income@geography$county, 3, "left", pad="0"), 
          str_pad(income@geography$tract, 6, "left", pad="0")), 
          income@estimate[,c("Household Income: Total:",
          "Household Income: $200,000 or more")], 
          stringsAsFactors = FALSE)
income_df <- select(income_df, 1:3)
rownames(income_df) <- 1:nrow(income_df)
names(income_df) <- c("GEOID", "total", "over_200")
income_df$percent <- 100*(income_df$over_200/income_df$total)
income_merged <- geo_join(tracts, income_df, "GEOID", "GEOID")
# There are some tracts with no land that we should exclude
income_merged <- income_merged[income_merged$ALAND>0,]
popup <- paste0("GEOID: ", income_merged$GEOID, "<br>", "Percent of Households above $200k: ", round(income_merged$percent,2))
pal <- colorNumeric(
    palette = "YlGnBu",
    domain = income_merged$percent
)
censusmap <- leaflet() %>%
    addProviderTiles(providers$Esri.WorldTopoMap) %>%
    addPolygons(data = income_merged, 
                fillColor = ~pal(percent), 
                color = "#b2aeae", # Need to use hex color codes.
                fillOpacity = 0.7, 
                weight = 1, 
                smoothFactor = 0.2,
                popup = popup) %>%
    addLegend(pal = pal,
              values = income_merged$percent, 
              position = "bottomright", 
              title = "Percent of Households<br>above $200k",
              labFormat = labelFormat(suffix = "%")) %>%
    addScaleBar(position = "topright")
# Display and save the map.
invisible(print(censusmap))

saveWidget(censusmap, file="censusmap.html", selfcontained=TRUE)