7 Leaflet

Main functions and concepts covered in this BP chapter:

census_api_key()
leaflet()
addPolygons()
GIS dataset
non GIS dataset
addLegend()

Packages used in this chapter:

## Load all packages used in this chapter
library(tidyverse) #includes dplyr, ggplot2, and other common packages

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(leaflet)
library(rmapshaper)
library(tidycensus)
library(tigris)

## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

options(tigris_use_cache = TRUE)
census_api_key("1666812db5ee2c29eac23238c5b2656fa3b42187")

## To install your API key for use in future sessions, run this function with `install = TRUE`.

## load census api key
## turn off scientific notation
options(scipen = 999)

7.1 GIS data download and preparation

Now we’re going to download the GIS/geographic data we need to create maps. See the comments in the code below for explanations of what each part is doing. I kept it all together in one code chunk with comments to make it easier to re-use elsewhere (e.g., as part of your CP).

Note that tidycensus/the Census API makes you download at least one variable of data. So, we are downloading the total population of each county. If we wanted to use Census data and download more than one variable, that would be part of our non-GIS dataset created first (i.e., the section above, just in this example, our only non-GIS data is the election data). This is what you’ll need to do for the CP. For the CP, you’ll download a bunch of variables using tidycensus and create a non-GIS dataset you use for the regressions and scatter plots. That dataset should NOT include the GIS data (i.e., use tidycensus options geometry = False, keep_geo_vars = FALSE). Then you’ll download separate a GIS dataset with just one variable, and merge the two at the end. Otherwise you will end up with really large file sizes and everything will be very slow (and might not work at all, especially when you try to commit/push to GitHub).

Make sure to also read the comments in the code.

Datasets used in this chapter:

## Load datasets used in this chapter

## Download GIS data for maps
##   geometry = TRUE --> include GIS shapefile data to create maps
##   B01001_001: total population (we have to download at least one variable, so we're using this one)
##   NOTE: If you want to download multiple variables for things like regressions, 
##     that's a non-GIS dataset and should use options: geometry = False, keep_geo_vars = FALSE
##     then merge that non-GIS dataset with a GIS dataset just like this one
##   We'll use year = 2023. The Census released the 2024 5-year ACS available in  December 2024, but it is not yet available via tidycensus as of when I created this file (January 19, 2025)
##   See the note at the bottom of this page about Connecticut and the difference between 2021 and earlier vs 2022 and later

countyGIS <- get_acs(geography = "county",
              variables = "B01001_001",
              year = 2023,
              geometry = TRUE,
              keep_geo_vars = TRUE)

## Getting data from the 2019-2023 5-year ACS

# State data (for displaying state borders on map)
stateGIS <- get_acs(geography = "state",
              variables = "B01001_001",
              year = 2023,
              geometry = TRUE,
              keep_geo_vars = FALSE)

## Getting data from the 2019-2023 5-year ACS

## In case there is ever a problem connecting to the Census data, when I ran this I ran the following save command:
# save(countyGIS, stateGIS, file = "countyGIS_stateGIS_2023.RData")
## In the event the API stops working, I'll be able to share the countyGIS_stateGIS_2023.RData object, which can be loaded with the following: 
# load(file = "countyGIS_stateGIS_2023.RData")

## Simplify GIS data to make file sizes smaller. ms_simplify essentially removes some details along coastlines and very-not-straight borders. It will make it look not quite as nice on some coasts, but overall it's not too noticeable and the reduced file size is important. Using object.size(countyGIS) before and after this step shows that it decreases countyGIS from around 20 MB to around 4 MB. This isn't a big deal in some contexts, but is important for GitHub and the HTML output loading well. 

# Note that if you use 0.01, the file size will be smaller, but it will remove a fair number of smaller counties (a lot of small city-counties in VA and a few others). Using 0.085 is the smallest number that doesn't remove counties in the continental US. 

stateGIS <- ms_simplify(stateGIS, keep = 0.085)
countyGIS <- ms_simplify(countyGIS, keep = 0.085)


countyGIS <- countyGIS %>% 
                dplyr::select(FIPS = GEOID, 
                       stFIPS = STATEFP, 
                       coFIPS = COUNTYFP, 
                       coNAME = NAME.x, 
                       pop = estimate, 
                       geometry)



## For mapping, let's drop the following: 
##   Puerto Rico (ST FIPS 72) (no election data)
##   Alaska (ST FIPS 02) (voting data isn't reported by county...we could also map the legislative districts, but we're not going to since we'd rather have smaller maps without those extra details)
##   Hawaii (ST FIPS 15) (so our map can zoom in on continental 48 states)
countyGIS <- countyGIS %>% filter(stFIPS != "72" & stFIPS != "02" & stFIPS != "15")
stateGIS <- stateGIS %>% filter(GEOID != "72" & GEOID != "02" & GEOID != "15")


## join 2-character state abbreviation and create name = "county, St" for labeling maps (e.g., Outagamie, WI) 
fipsToSTcode <- fips_codes %>% dplyr::select(stFIPS = state_code, stNAME = state) %>% unique()

countyGIS <- inner_join(countyGIS,fipsToSTcode,by="stFIPS")

countyGIS <- countyGIS %>% mutate(name = paste0(coNAME,", ", stNAME))



## NOTE: If you don't use keep_geo_vars = TRUE, you don't get separate STATEFP and COUNTYFP, but you can use mutate() and create stFIPS = substr(GEOID,1,2) and coFIPS = substr(GEOID,3,5)

7.2 Popup labels

We also need to make the text (with HTML for formatting) to use in popup labels for each county. First we create the text that will make up the labels, along with HTML formatting (e.g., <b> to make the county name bold, font color to make it red when Trump wins and blue when Haris wins). Then we pipe that to the HTML function from the htmltools package that turns our text into HTML code that leaflet can use to make popups in the map.

Note that you do this with the GIS dataset. This is crucial! If you use different datasets to make the map and create the popup labels, the popup labels might not match with the correct counties.

popupLabels <- paste0(
  "<b>", countyGIS$name, " (", countyGIS$FIPS, ")</b>",
  "<br><font color='black'>Population: ",
  format(countyGIS$pop, big.mark = ",", trim = TRUE),
  "</font>"
) %>%
  lapply(htmltools::HTML)

# big.mark = "," adds a comma to the number format

Note that in this example, we’re creating the popup labels one time and we’ll use them in each map (since all the maps show the same data, just using different color schemes). For the CP you will create several maps that show different variables. For that, you need to create new popup labels for each map.

Also note that it is important to use the same dataframe for each part. A common mistake people make is to do something like countyGIS$name for the name and countyGIS$totalVotes for total votes, but nonGIS$pctWinner for the other value. Mixing dataframes like this is very dangerous. Sometimes it works correctly, but other times you end up with data from the wrong county. If they have a different number of rows, usually it won’t work correctly, and you’ll end up displaying nonsense (i.e., the values for the wrong counties…picture a map that labels our county, instead of Outagamie County, WI, it labels it Travis County, TX…in other words, your map becomes complete nonsense). By using the same dataframe (i.e., countyGIS$pctWinner instead of nonGIS$pctWinner), you guarantee that the same county will be used for everything.

I suggest you try to make changes to the popup labels below and see what they look like in the map. Note that the parts in <> are HTML formatting. For example <b>text</b> makes the word “text” bold (like this). You aren’t expected to already know any HTML. You are, however, expected to be able to look online and figure it out. (Remember, becoming confident figuring out new skills on your own is one of the main goals of the course.) For the maps you create for the CP and PP, you are expected to properly update the formatting of the popup labels. Play around with the one above, playing close attention to every single symbol. Every “<”, “>”, “(”, “)”, “%”, etc, is important and must be correct.

Note that if when you create a map you get an error message similar to “cannot coerce type ‘closure’ to vector of type ‘character’”, it might be a problem with your popup labels.

7.3 Population Map

On this page we’ll build our first interactive web-based map of geospatial data using Leaflet. You have most likely interacted with maps made with Leaflet when reading articles online or on the websites of businesses, non-profits, universities, and government agencies. After working through this “book” you should be able to make leaflet maps yourself (and you’ll need to for the CP).

Leaflet was initially a javascript library used to include interactive maps on websites. The leaflet R package allows you to create Leaflet maps from R, combining the data-analytic capabilities of R with the ability to create maps communicating your analysis that can be put directly on websites.

## create color palette used by map
pal <- colorBin("Oranges", countyGIS$pop, bins = 10, reverse=TRUE,na.color = "black")

leaflet(countyGIS, options = leafletOptions(crsClass = "L.CRS.EPSG3857"), width="100%") %>%
  addPolygons(weight = 0.4, color = "lightgray", opacity = 0.7,
    fillColor = ~pal(pop), fillOpacity = 1, smoothFactor = 0.5,
    label = popupLabels,
    labelOptions = labelOptions(direction = "auto")) %>%
    addPolygons(data = stateGIS,fill = FALSE,color="black",weight = 1) %>%
    addLegend(pal = pal,values = ~countyGIS$pop, opacity = 0.7, title = "Population",position = "bottomright")

Using the leaflet() call we can create interactive maps of our data. Below see an outline of what each portion of the above code does.

The width argument handles how wide to print the graphic. 100% prints the standard size of a graphic in an Rmd file. crsClass = "L.CRS.EPSG3857" tells leaflet to use the standard map projection

The addPolygons function tells leaflet to draw in the shapes, in this case the county lines. weight controls the thickness of the lines, the larger the number the thicker the line. color determines the color of the outline. opacity controls the transparency with 1 being solid and 0, completely transparent. fillColor fills the counties appropriately by color pallette ~pal(). fillOpacity is similar to opacity with 1 being solid and 0 being transparent. smoothFactor does some additional smoothing to the shapes to help it render faster. label informs leaflet of what to include in the floating label. Finally, labelOptions = labelOptions(direction = "auto") tells leaflet to place the label wherever it fits best, keeping it from going off screen.

7.4 Note about which states are included

In this file, we are going to create county-level maps for the 48 continental US states. For the CP and PP, you are going to also create county-level maps of the 48 continental US states. However, for your regressions in the CP and PP, you should also include Alaska and Hawaii. We are only dropping Alaska and Hawaii for the maps (and Alaska doesn’t report voting on a county-level basis, so if we were doing regressions, we couldn’t include it…but for the CP and PP you can and should).

7.5 Note about Connecticut

Connecticut used to have 8 counties. For 2021 and earlier data, there will be 8 rows for Connecticut. Starting in 2022, they switched to 9 planning regions instead of the 8 counties. So, for 2022 and later, there are 9 rows for Connecticut. See this page for more details: https://www.federalregister.gov/documents/2022/06/06/2022-12063/change-to-county-equivalents-in-the-state-of-connecticut

Unfortunately, this makes working with data from different time periods difficult (or impossible). The 8 counties have different borders (and different FIPS codes) than the 9 planning regions (see the map in that link). If you wanted to make the same map but for 2020 data, you’d have to switch the census data to 2020 (or 2021) so that it would have 8 counties for Connecticut. To work with the 2024 election data, you have to use 2022 or 2023 Census data so that Connecticut has 9 counties. If you wanted to do something like compare the 2024 and 2020 election data (say, % Trump 2024 - % Trump 2020 to look at the change), you can’t do it for Connecticut with this data. The only way to do it would be to find another data source that has smaller areas within Connecticut and then aggregate them into the current 9 planning regions (or into the previous 8 counties). That is possible for some data, but not all data (and clearly works with counts, but not for things like percentages or medians).

This may come up in the PP and the RP, so it’s important to be aware of so if you later have problems with Connecticut, you know why.