A striking feature of many maps from early in the history of cartography is their linearity. Being primarily for travel (and given the technological limitations on how faithfully geographies could be understood and displayed), maps from centuries past often took the form of long scrolls or spans of parchment displaying the ordered features of routes (which river would be encountered after which village, and so on). The Tabula Peutingeriana, possibly the only extant map of the ancient Roman road system, is more than twenty feet long despite being a far slimmer foot or so in width. It reflects an approach to cartography that treats maps as itineraries, as ordered schematics. The linear design has a delightful functionality: The geographic realism we expect from a contemporary map (who nowadays would represent the journey from Denver to Aspen as a horizontal line?) is subjugated to the practical concern of safely finding one’s way from Florence to Rome.

But maps do much more than give travel instructions these days. Data scientists and journalists put them to work representing quantitative differences across space and time, from the disparate effects of climate change on richer and poorer neighborhoods to the amount of time people were spending at home across the US early on in the COVID-19 pandemic. For these purposes, linearity just won’t do: We don’t need to represent the journey between two points; we need to represent data about many areas simultaneously.

Doing so in R is alluringly easy with Leaflet. Leaflet is an open-source JavaScript library for making interactive maps. Its use is simple: The user creates a map widget and then layers features onto that map widget until the display and interactivity are as desired. (This process will feel familiar to users of the R community’s darling package ggplot2.)

In this article, I’ll review two common ways of representing data on a map using Leaflet in R: (1) Defining regions on a map and distinguishing them based on their value on some measure using colors and shading (“choropleths,” these maps are called), and (2) marking individual points on a map (e.g., archaeological dig sites, baseball stadiums, voting locations, etc.).

If you plan to code along with the examples below, be sure to install and load the leaflet package. I also recommend installing and/or loading the magrittr package so that you can make use of the pipe operator (%>%) to chain functions, which proves to simplify the map-building process greatly. Install and/or load the following packages too; I call on them at various points in the code below: sf, geojsonio, htmltools, htmlwidgets, stringi, and RColorBrewer.



The goal: Build a choropleth to represent the average price of electricity (in cents per kilowatt-hour) across each of the 48 contiguous US states (+ DC) in 2018. The data are sourced from the US Energy Information Administration, and you can access them using the code below.1

A peek at the data:

dat <- read.csv('https://static.lib.virginia.edu/statlab/materials/data/state_electricity_data_2018.csv')
        NAME centskWh
1    Alabama     9.63
2     Alaska    19.36
3    Arizona    10.85
4   Arkansas     7.78
5 California    16.58
6   Colorado    10.02

Whenever we’re working with Leaflet, we use the function leaflet() to initialize a map widget.

m <- leaflet()

Running m at this point will simply load a gray panel void of any geographic features. This is the hollow shell of the eventual choropleth.

Before we define our areas of interest (states) and layer on data (electricity costs), the map widget needs a “basemap” to serve as its foundation. Basemaps consist of map tiles—individual map sections that join together to form a composite picture. The simplest way to add a basemap to the map widget is with addTiles(). Running this without adjustment will result in an overhead view of the whole world. Since the data at hand are for US states, the map should be centered on America: The setView() function allows the user to specify a longitude (east-west) and latitude (north-south) to serve as the initial center of the map. (The geographic center of the contiguous US? Not far from the Kansas-Nebraska border.)

By default, addTiles() generates a basemap using OpenStreetMap tiles. They’re suitable, but the options of which tiles to use are extensive. Check out a lengthy list here. I’m partial to the washed simplicity of the “CartoDB.PositronNoLabels” tiles, so I’ll use those in the maps that follow. You can use addProviderTiles() instead of addTiles() to pick from among a prespecified set of third-party options, or you can call a tile set from a URL by providing a link in addTiles().

m <- leaflet() %>%
  addProviderTiles(providers$CartoDB.PositronNoLabels)  %>% 
  setView(lng = -96.25, lat = 39.50, zoom = 4)

With that, we have a basemap. But we need to represent each state as a distinct area that we can associate with an electricity cost. We can feed Leaflet the necessary spatial data using shapefiles. The Census Bureau is an excellent source of shapefiles for the US, and the shapefiles it offers cover many levels of geographic granularity (states, counties, school districts, congressional districts, and more).

For these maps, we’ll use the Census Bureau’s 2019 1:5,000,000 US states shapefile, which you can directly download here. The file you’ll need to read in ends with .shp, and you can read it in with read_sf() from the sf package. (Note: The folder you download from the Census Bureau will contain other files—.dbf, .shx, and more—in addition to the .shp file. You’ll need to keep those files in the same directory from which you read in the .shp file. You can learn more about the structure of shapefiles here.)

states <- read_sf('cb_2019_us_state_5m/cb_2019_us_state_5m.shp')

The shapefile has been assigned to the object states, and it contains a geometry column with the spatial data. With it, we can now define each state on the basemap using the addPolygons() function:

m <- leaflet() %>%
  addProviderTiles(providers$CartoDB.PositronNoLabels)  %>% 
  setView(lng = -96.25, lat = 39.50, zoom = 4) %>% 
  addPolygons(data = states,
              weight = 1)