This tutorial demonstrates how to create interactive choropleth like the one below using the leaflet library, there are no alternative libraries currently considered in this guide.

It is important to note that such maps are generated from two datasets: points on the map and region shape files, the colours are then assigned according to the number of points in each region.

Shapefiles for the regions of interest as supported by readOGR (see below for further details)
A .csv file contain a set of geographic points with decimal formatted longitude and latitude coordinates

example <- structure(list(longitude = c(-104.990251, -104.990251, 13.960129, 
13.960129, 13.960129), latitude = c(39.739236, 39.739236, 50.97178, 
50.97178, 50.97178)), .Names = c("longitude", "latitude"), row.names = c(NA, 
5L), class = "data.frame")
## kable is loaded from the knitr library to convert the data.frame into a Markdown table for display purposes
knitr::kable(example)

longitude	latitude
-104.99025	39.73924
-104.99025	39.73924
13.96013	50.97178
13.96013	50.97178
13.96013	50.97178

Importing Shapefiles

Shapefiles will be imported with the readOGR function and the stored as a SpatialPolygonsDataFrame - this workflow requires what we might term ESRI or DBF formatted shapefiles. An informative wikipedia page on this format is available here, the first check for whether your shapefiles are correctly formatted is whether you have a directory containing the following files:

shapefile_name.dbf
shapefile_name.shx
shapefile_name.shp
shapefile_name.prj (if this file is not provided additional work is required to discover the projection of the shapefiles)

In this example the US state boundaries as provided by https://www.census.gov/geo/maps-data/data/tiger-cart-boundary.html are used. This is primarily because the US government provides an excellent repository for correctly formatted shapefiles and the US states are large regions that are easy to compare, unlike the UK administrative borders.

The shapefiles are imported using readOGR from the rgdal library; specifying the directory containing the shape files to dsn and the shapefile_name to layer (note the use of verbose = F to prevent uninteresting messages from the function). Details of the leaflet library are covered later, but for demonstrative purposes it is useful to see the US state boundaries are clearly visible:

library(rgdal)
full_us_shapefiles <- readOGR(dsn = "data/shapefiles", layer = "cb_2015_us_state_500k", verbose = F)
library(leaflet)
leaflet(data = full_us_shapefiles) %>% addTiles() %>% 
  addPolygons(
  weight = 1
)

Unfortunately, the shapefiles include all US terrorities. Often it is only the contiguous United States that are displayed in visualisations, and this does make a very useful example for this guide. Fortunately, it is fairly easy to extract this set of states through the use of the US FIPS Codes which are documented here and this data is stored in the full_us_shapefiles object but need to be converted to numeric format for use later:

full_us_shapefiles$STATEFP <- as.numeric(as.character(full_us_shapefiles$STATEFP))
full_us_shapefiles$STATEFP

##  [1] 31 53 35 46 21 13  5 42 28  8 49 47 56 18  2 32 17 50 30 19 45  4 60
## [24] 39  1 55 41 29 37 40 51 54 22 26 25 16 72 12 20 33 10 48 27 78  9 34
## [47] 38 24 23 15 66 69 44  6 36 11

The function make_us_contiguous uses the US-FIPS-Codes.csv file hosted on Figshare to extract only those polgyons in full_us_shapefiles that are in the contiguous USA:

fips_codes <- read.csv("https://ndownloader.figshare.com/files/5548022", stringsAsFactors = F)
make_us_contiguous <- function(spatial_polgyon= NA){
  contiguous_fips_codes <- fips_codes[fips_codes$Contiguous.United.States. == "Y",]$STATE
  contiguous <- full_us_shapefiles[spatial_polgyon$STATEFP %in% contiguous_fips_codes,]
  # Drop unnecessary levels
  contiguous@data <- droplevels(contiguous@data)
  contiguous
}
contiguous_us_shapefiles <- make_us_contiguous(full_us_shapefiles)

The fips_codes object also contains the names of the states, which may be added to the contiguous_us_shapefiles object and can be visualised with leaflet once more. Note that the function mapvalues from the plyr library is used as it is massively convenient for modifying a vector according to a set of key-value pairs.

library(plyr)
contiguous_fips_codes <- fips_codes[fips_codes$Contiguous.United.States. == "Y",]
contiguous_us_shapefiles$State_Name <- mapvalues(contiguous_us_shapefiles$STATEFP, 
                                                 from = contiguous_fips_codes$STATE, 
                                                 to = contiguous_fips_codes$STATE_NAME)
leaflet(data = contiguous_us_shapefiles) %>% addTiles() %>% 
  addPolygons(
  weight = 1
)

Count Points Per Location

In order to count the number of points within each state it is necessary to create a SpatialPointsDataFrame, in the code below the lat-long-points.csv file is first imported and then converted into such an object. It is necessary to specify the projection for the SpatialPointsDataFrame, which is simply extracted from the contiguous_us_shapefiles object.

locations_df <- read.csv("data/lat-long-points.csv", stringsAsFactors = F)
locations_spdf <- SpatialPointsDataFrame(coords = locations_df,
                                         data = locations_df,
                                         proj4string = contiguous_us_shapefiles@proj4string)

The GISTools library includes poly.counts which provides a simple way to count the number of points within each region and this data can easily be added to contiguous_us_shapefiles in the same manner as the State names:

contiguous_counts <- poly.counts(pts = locations_spdf, polys = contiguous_us_shapefiles)
contiguous_us_shapefiles$Count.of.locations <- as.numeric(contiguous_counts)

Basic Leaflet Map of Shapefiles

Leaflet maps are created in layers and very much favour the pipe “%” operator, the most basic map of the world can be constructed in one of two ways:

base_map <- leaflet()
addTiles(base_map)

## or
leaflet() %>% addTiles()

The shapefiles are added to the map using addPolygons, in the code below a number of optional arguments have been supplied to customise the appearance of the borders.

leaflet(data = contiguous_us_shapefiles) %>% 
  addTiles() %>% 
  addPolygons(
  stroke = TRUE,
  color = "#000000",
  smoothFactor = 0.2,
  weight = 1
)

Before creating a colour palette for the map, thereby makign it a choropleth, it is useful to note how popups can be added to the map. In the map below the State_Name will be shown when a state is clicked, the tilda (~) allows the namespace of the contiguous_us_shapefiles object to be accessed.

leaflet(data = contiguous_us_shapefiles) %>% 
  addTiles() %>% 
  addPolygons(
  stroke = TRUE,
  color = "#000000",
  smoothFactor = 0.2,
  weight = 1,
  popup = ~State_Name
)

Custom Choropleth Color Palette

The RColorBrewer library provides a number of versatile and accessible color palettes, see colorbrewer2.org for a breakdown of these, we can use in the choropleth via brewer.pal. In the code below a six colour palette is defined over the range 0 - 350 with colorBin, note the use of alpha to ensure that coloured regions are opaque to the map details behind them.

palette <- colorBin(c("#cccccc",brewer.pal(5, "YlGnBu")), 
                    bins = c(0,1,5,10,20,50,350),
                    alpha = TRUE
                    )

This fill colour of the polygons can now be specified with palette(Count.of.locations):

leaflet(data = contiguous_us_shapefiles) %>% 
  addTiles() %>% 
  addPolygons(
  stroke = TRUE,
  color = "#000000",
  smoothFactor = 0.2,
  fillOpacity = 0.6,
  fillColor = ~palette(Count.of.locations),
  weight = 1,
  popup = ~State_Name
)

The popup displayed at the moment in the map contains only the state name, the following utility function allows the state name and the number of data points in the location to be displayed:

region_labeller <- function(state_name = NA, number_of_points = NA){
  paste0(
    "<p>State: ", state_name, "</p>",
    "<p>Number of Points: ", number_of_points, "</p>"
  )
}
choropleth_map <- leaflet(data = contiguous_us_shapefiles) %>% 
  addTiles() %>% 
  addPolygons(
  stroke = TRUE,
  color = "#000000",
  smoothFactor = 0.2,
  fillOpacity = 0.6,
  fillColor = ~palette(Count.of.locations),
  weight = 1,
  popup = ~region_labeller(State_Name,Count.of.locations)
)
choropleth_map

Choropleth Legend

The leaflet library provides a very flexible legend, which can be provided with our custom colour scheme as follows:

choropleth_map %>% addLegend(
  position = 'bottomleft',
  colors = c("#cccccc",brewer.pal(5, "YlGnBu")),
  labels = c("0","1-5","5-10","10-20","20-50","50-350"),
  opacity = 0.6,
  title = "Number of Points in State")   ## title of the legend

htmlwidget Choropleth from Shapefiles and Points

Oxford University Interactive Data Network

Importing Shapefiles

Count Points Per Location

Basic Leaflet Map of Shapefiles

Custom Choropleth Color Palette

Choropleth Legend