mapbox dot density

These days some of the best, most complicated interactive maps are being published using Mapbox.

Raster maps are so 2012.


It feels like a natural evolution. In fact, a lot of the same people who developed Leaflet went on to work on Mapbox.

What’s the difference?

Leaflet uses Raster Tiles while Mapbox uses Vector Tiles.

It’s like creating graphics in Adobe Photoshop versus Adobe Illustrator.

Raster versus Vector

  • Elements of a vector tile can be easily styled, like roads.

  • Vector tiles are rendered on the user end, not on the server end

  • It’s just faster. Vector files load 20-50 percent quicker.

Mapbox on R with Mapdeck

Just as someone made an R package that interfaces R with Leaflet, someone made a package in 2019 to interface with Mapbox: %>% %>% Mapdeck

It’s still really new. And it’ll take some time for people to create new packages that build on top of it in R. But it’s something that will definitely continue to grow in popularity.

Downside: It’s not unlimited. You’ll be charged if your map gets a lot of use.

Let’s reproduce this Washington Post map

Remember this story about diversity and segregation?

Let’s see if we can get close to reproducing it.

The repo containing the data and scripts for this section is on Github. To follow along, simply run the lines of code below in R.


# There is no need to run these install lines below if you already have done so
install.packages("usethis")
usethis::use_course("https://github.com/andrewbtran/NICAR-2020-mapping/archive/master.zip")

# This section is in scripts/05_mapbox_dotdensity.R
file.edit("scripts/05_mapbox_dotdensity.R")

Get an API key

  1. Sign up for a Mapbox account and generate an Access token

Next, we need to load the right packages.

# This function checks if you don't have the correct packages installed yet
# If not, it will install it for you
packages <- c("sf", "tidyverse", "tigris",
              "tidycensus", "mapdeck")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(packages, rownames(installed.packages())), repos = "https://cran.us.r-project.org")  
}
library(sf)
library(tidyverse)
library(tigris)
library(tidycensus)
library(mapdeck)

The next steps involve querying the Census API with tidycensus for blockgroup-level data on population by race.

It takes a while to download and requires your own Census API key, so go ahead and just load the pre-downloaded data with the code below.

acs <- readRDS("data/acs_geoid.RDS")

Here’s the R code to query and download that data on your own if you want to pull for another state (or ALL of them if you’ve got a lot of time set aside). Be sure 1) you have a solid internet connection and 2) you have a Census API key that you can substitute in.

census_key <- "YourCensusKeyGoesHere"

census_api_key(census_key)

## this function downloads the data by variables we've set
acs <- get_acs("block group", table = "B02001", cache_table = TRUE,
               geometry = TRUE, state = "22",
               year = 2016, output = "tidy")

## These lines cleans up the downloaded data frame
acs <- acs %>%
  mutate(
    id = str_extract(variable, "[0-9]{3}$") %>% as.integer
  ) %>%
  # variable 1 is the "total", which is just the sum of the others
  filter(id > 1 & id < 8) %>%
  mutate(race =case_when(
    id == 2 ~ "White",
    id == 3 ~ "Black",
    id == 4 ~ "Native",
    id == 5 | id == 6 ~ "Asian",
    id == 7 ~ "Other"
  )) %>% 
  group_by(GEOID, race) %>% 
  summarise(estimate = sum(estimate))

This is what the downloaded data from the Census looks like:

glimpse(acs)
Rows: 17,355
old-style crs object detected; please recreate object with a recent sf::st_crs()
Columns: 4
Groups: GEOID [3,471]
$ GEOID    <chr> "220019601001", "220019601001", "220019601001", "220019601001…
$ race     <chr> "Asian", "Black", "Native", "Other", "White", "Asian", "Black…
$ estimate <dbl> 0, 109, 0, 0, 1121, 0, 1088, 0, 0, 412, 0, 0, 0, 0, 587, 0, 3…
$ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-92.22295 3..., MULTIPOLYGON (((…

Notice the geometry column? That’s where the polygon data resides but in a tidy format.

Let’s quickly render this on a map to see what we’ve downloaded from the Census API and what we’re working with.

First, let’s look at just one variable in the map.

# let's look at one variable so it doesn't throw off the map
acs_black <- acs %>%
  filter(race=="Black")
old-style crs object detected; please recreate object with a recent sf::st_crs()
glimpse(acs_black)
Rows: 3,471
Columns: 4
Groups: GEOID [3,471]
$ GEOID    <chr> "220019601001", "220019601002", "220019601003", "220019601004…
$ race     <chr> "Black", "Black", "Black", "Black", "Black", "Black", "Black"…
$ estimate <dbl> 109, 1088, 0, 378, 140, 0, 117, 0, 31, 13, 0, 39, 104, 0, 0, …
$ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-92.22295 3..., MULTIPOLYGON (((…

Let’s load the mapdeck library and put in the Mapbox API key.

library(mapdeck)

mb_key <- "PutYourKeyHere"
mapdeck(token = mb_key, style = mapdeck_style("dark")) %>%
  add_polygon(
    data = acs_black,
    fill_colour = "estimate",
    fill_opacity = .4
  )
Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite