How to make web-ready US county-level maps

Here, I will describe 3 workflows to get you up and running in making fast, web-ready maps.

Asmae Toumi https://twitter.com/asmae_toumi
10-13-2020

Table of Contents


Code is up on my GitHub at https://github.com/asmae-toumi/spatialR and here I’ll rapidly go over the two workflows I like. If you see something wrong or need more explaining, hit me up on Twitter.

Workflow 1: Mapbox

Mapbox has been picking up tons of steam as of late. You might’ve seen that it’s what the New York Times uses for their COVID-19 maps: https://www.nytimes.com/interactive/2020/us/coronavirus-us-cases.html. I was impressed with how rapidly they loaded in my browser and how fast the zooming was. I decided to figure out how to make a map with Mapbox because the Leaflet ones I was making for US county-level data were just too slow.

Here’s how I did it with Mapbox:

Step 1: Set up

Kyle Walker made an awesome package that allows us to interface with Mapbox’s API called {mapboxapi}. Install it in R:


remotes::install_github("walkerke/mapboxapi")

Go to Mapbox and create an account. Generate a token and save it with mb_access_token() like so:


library(mapboxapi)
mb_access_token("pk.eyas...", install = TRUE)

Step 2: Get geometries

Fetch the geometries you need. Here, I get US counties from {tigris}. Then, I join them to my data and finally simplify the geometries to make them more lightweight using ms_simplify() from the {rmapshaper} package.


library(tidyverse)
library(tigris)
library(sf)
library(rmapshaper)

data <- 
  tigris::counties(class = "sf") %>% 
  geo_join(data, by_sp = "GEOID", by_df = "GEOID") %>% 
  ms_simplify()

Your data is of class sf now!

Step 3: Optimize geometries

Our data is big, not that big but big enough that if you wanted to upload it directly to Mapbox Studio it probably wouldn’t let you. This is when Tippecanoe comes in. Tippecanoe converts your data into much smaller files. It does so by creating vector tilesets from large GeoJSON feature collections. The output is an “mbtiles” file that can be uploaded to Mapbox. I think this tool is the secret behind why Mapbox maps are lightning fast.

You need to install Tippecanoe via your terminal using homebrew. If you don’t have homebrew, go here and download it: https://brew.sh.

In your terminal:


brew install tippecanoe

Once you have it installed on your computer, use the tippecanoe() command from {mapboxapi} to create your mbtiles file like so:


tippecanoe(
  input = my_tiles,
  output = "my_tiles.mbtiles",
  layer_name = "my_tiles") 

Step 4: Upload to API

Once your mbtiles file is created, you need to upload it to Mapbox using upload_tiles():


upload_tiles(input = "my_tiles.mbtiles",
             username = "YOUR_USERNAME", 
             tileset_id = "my_tiles",
             multipart = TRUE)

Go to Mapbox Studio. You should see your mbtiles uploaded. Start styling!

If you want to start from something, you can tinker with mine by adding my map to your styles using this link: https://api.mapbox.com/styles/v1/atoumi/ckddsbc903xey1io2y1nxiowx.html?fresh=true&title=copy&access_token=pk.eyJ1IjoiYXRvdW1pIiwiYSI6ImNrZGFqM3lhZzBpYjkyeXJ5djVxdXozZ2gifQ.N6IaEvOeboUvhKxjgSDPnQ

Step 5: Bring it back to R

When you’re done styling, you can either publish directly or bring it back to R using the {mapdeck} package.


mapdeck(token = Sys.getenv("MAPBOX_PUBLIC_TOKEN"),
        style = "mapbox://styles/atoumi/ckddsbc903xey1io2y1nxiowx",
        zoom = 6,
        location = c(-98.7382803, 31.7678448)) 
        # if it doesn't show up in your viewer, view in browser

Now that it’s in R, you can further edit it and even place it in a {Shiny} app using mapdeckOutput() in the UI and renderMapdeck() in the server!

The downside with this workflow is that if you want to add a legend or hover/popup elements, you need to get acquainted with Mapbox GL JS, which is Mapbox’s JavaScript library that uses WebGL to render interactive maps from vector tiles and Mapbox styles. You can see how I made my map interactive with some JS by checking out the index.html in my GitHub. Be sure to include your own access token before rendering it. You can also shove the entire index.html file inside {Shiny} like this:


# UI ----------------------------------------------------------------------
ui <- navbarPage("Fun with Mapbox", id="nav",

                 tabPanel(
                   "Interactive map",
                   includeHTML("index.html")
                 )
)


# Server ------------------------------------------------------------------

server <- function(input, output) { }

shinyApp(ui = ui, server = server)

Another downside is that Mapbox is going to start charging money soon for uploading tiles to their API. This is why I like workflow 2, because nothing is better than FREE.

Workflow 2: r2d3map

Using d3_map():

The amazing dreamRs team is behind the {r2d3maps} package which allows you to make D3 maps in a flash. Check out their documentation to get a taste of what you can make: https://github.com/dreamRs/r2d3maps. You can also try out my fully worked out examples on Github. I have a regular R script, and a quick demo on {Shiny}. It’s absurdly easy to get your map up on the web using {Shiny}:


# UI ----------------------------------------------------------------------

ui <- navbarPage("Fun with r2d3map", id="nav", 
                 
                 tabPanel(
                   "Interactive map",
                   withSpinner(d3Output(
                     outputId = "mymap", 
                     width = "900px", 
                     height = "500px"))),
                 
                 tabPanel("Explore the data",
                          DT::dataTableOutput("table"))
)


# Server ------------------------------------------------------------------

server <- function(input, output) {
  
  # map panel 
  output$mymap <- renderD3({
    d3_map(shape = cty_sf_joined, projection = "Albers") %>%
      add_labs(caption = "Viz: Asmae Toumi | Data: CDC") %>% 
      add_continuous_breaks(var = "rpl_theme1", palette = "Reds") %>%
      add_legend(title = "Socioeconomic Vulnerability (percentile)") %>%
      add_tooltip("<b>{location}</b>: {rpl_theme1}")
  })
  
  # data panel
  output$table <- DT::renderDataTable({
     DT::datatable(svi, rownames = F,  filter = 'top',
                  extensions = c('Buttons', 'FixedHeader', 'Scroller'),
                  options = list(pageLength = 15, lengthChange = F,
                                 fixedHeader = TRUE,
                                 dom = 'lfBrtip',
                                 list('copy', 'print', list(
                                   extend = 'collection',
                                   buttons = c('csv', 'excel', 'pdf'),
                                   text = 'Download'
                                 ))
                  ))
  })
  
}

shinyApp(ui = ui, server = server)

Make your own

You don’t have to use d3_map(), you can make your own using pure D3 code.

First, you need to convert your data from an sf object to a topojson object using the handy r2d3map() command:


r2d3map(
  data = my_data,
  script = "my_map.js"
)

Then, call use_r2d3map() to create a minimal template. It will create the 3 scripts you need to get your map working: an R script, a JS script and a CSS script.


use_r2d3map("my_map.js")

A fun exercise would be to try and recreate this striking chloropleth using use_r2d3map(): http://bl.ocks.org/syntagmatic/623a3221d3e694f85967d83082fd4a77.

Let me know if you do!

Workflow 3: Leaflet

I was wrong! {Leaflet} doesn’t have to be slow. The {tigris} package offers county shapefiles, and the secret is in specifying that you want lower resolution shapefiles. Let’s make a county-level map of U.S. income data from the {tidycensus} package. First, we obtain the county-level income data from the 2014-2018 5-year American Community Survey (ACS) which we then join to the cty dataset using tigris::geo_join :


library(tigris)
library(tidyverse)

cty <- counties(cb = TRUE, resolution = "20m") %>% rename(fips = GEOID)

library(tidycensus)

us_county_income <- 
  get_acs(geography = "county", variables = "B19013_001") %>% 
  select(fips = GEOID, name = NAME, income = "estimate") %>% 
  drop_na()

cty <- cty %>% 
  geo_join(us_county_income, by_df = "fips", by_sp = "fips")

Let’s make the map now:


library(leaflet)
library(RColorBrewer)
library(htmltools)

pal <- colorNumeric("YlOrRd", domain = cty$income)

labels <- 
  sprintf(
    "<strong>%s</strong><br/> Income: %s $",
    cty$name, cty$income) %>% 
  lapply(htmltools::HTML)

leaflet(cty) %>%
  setView(-96, 37.8, 4) %>%
  addProviderTiles("CartoDB.PositronNoLabels") %>%
  addPolygons(
        fillColor = ~pal(income), 
        weight = 1,
        opacity = 1,
        color = "white",
        fillOpacity = 0.7,
        highlight = highlightOptions(
          weight = 2,
          color = "#666",
          fillOpacity = 0.7,
          bringToFront = TRUE),
        label = labels,
        labelOptions = labelOptions(
          style = list("font-weight" = "normal"),
          textsize = "15px",
          direction = "auto")) %>%
  addLegend(pal = pal,
            values = cty$income,
            position = "bottomright",
            title = "Income (5-year ACS)",
            labFormat = labelFormat(suffix = "$"),
            opacity = 0.8,
            na.label = "No data")