R and RStudio

• R: programming language and software environment for data analysis (and wherever else your imagination can take you!)

• RStudio: integrated development environment (IDE) for R developed by Posit

• Posit Cloud: run RStudio with today’s workshop pre-configured at https://posit.cloud/content/7549022

The Decennial US Census

• Complete count of the US population mandated by Article 1, Sections 2 and 9 in the US Constitution

• Directed by the US Census Bureau (US Department of Commerce); conducted every 10 years since 1790

• Used for proportional representation / congressional redistricting

• Limited set of questions asked about race, ethnicity, age, sex, and housing tenure

The American Community Survey

• Annual survey of 3.5 million US households

• Covers topics not available in decennial US Census data (e.g. income, education, language, housing characteristics)

• Available as 1-year estimates (for geographies of population 65,000 and greater) and 5-year estimates (for geographies down to the block group)

• Data delivered as estimates characterized by margins of error

How to get Census data

• data.census.gov is the main, revamped interactive data portal for browsing and downloading Census datasets

• The US Census Application Programming Interface (API) allows developers to access Census data resources programmatically

tidycensus

tidycensus: key features

Getting started with tidycensus

• To get started, install the packages you’ll need for today’s workshop

• If you are using the Posit Cloud environment, these packages are already installed for you

install.packages(c("tidycensus", "tidyverse", "mapview", "ggspatial", "leafsync"))

Optional: your Census API key

• tidycensus (and the Census API) can be used without an API key, but you will be limited to 500 queries per day

• Power users: visit https://api.census.gov/data/key_signup.html to request a key, then activate the key from the link in your email.

• Once activated, use the census_api_key() function to set your key as an environment variable

library(tidycensus)

census_api_key("YOUR KEY GOES HERE", install = TRUE)

Spatial Census data: the old way

Traditionally, getting “spatial” Census data required:

Getting started with tidycensus

• Your core functions in tidycensus are get_decennial() for decennial Census data, and get_acs() for ACS data

• Required arguments are geography and variables

texas_income <- get_acs(
  geography = "county",
  variables = "B19013_001",
  state = "TX",
  year = 2022
)

• ACS data are returned with five columns: GEOID, NAME, variable, estimate, and moe

texas_income

# A tibble: 254 × 5
   GEOID NAME                    variable   estimate   moe
   <chr> <chr>                   <chr>         <dbl> <dbl>
 1 48001 Anderson County, Texas  B19013_001    57445  4562
 2 48003 Andrews County, Texas   B19013_001    86458 16116
 3 48005 Angelina County, Texas  B19013_001    57055  2484
 4 48007 Aransas County, Texas   B19013_001    58168  6458
 5 48009 Archer County, Texas    B19013_001    69954  8482
 6 48011 Armstrong County, Texas B19013_001    70417 14574
 7 48013 Atascosa County, Texas  B19013_001    67442  4309
 8 48015 Austin County, Texas    B19013_001    73556  4757
 9 48017 Bailey County, Texas    B19013_001    69830 13120
10 48019 Bandera County, Texas   B19013_001    70965  5710
# ℹ 244 more rows

Spatial Census data with tidycensus

• Use the argument geometry = TRUE to get pre-joined geometry along with your data!

texas_income_sf <- get_acs(
  geography = "county",
  variables = "B19013_001",
  state = "TX",
  year = 2022,
  geometry = TRUE
)

• We can now make a simple map with plot():

plot(texas_income_sf['estimate'])

Looking under the hood: simple features in R

• Spatial data are returned with five data columns: GEOID, NAME, variable, estimate, and moe, along with a geometry column representing the shapes of locations

texas_income_sf

Simple feature collection with 254 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -106.6456 ymin: 25.83738 xmax: -93.50829 ymax: 36.5007
Geodetic CRS:  NAD83
First 10 features:
   GEOID                     NAME   variable estimate   moe
1  48273    Kleberg County, Texas B19013_001    52487  4987
2  48391    Refugio County, Texas B19013_001    54304  2650
3  48201     Harris County, Texas B19013_001    70789   482
4  48443    Terrell County, Texas B19013_001    52813 11107
5  48229   Hudspeth County, Texas B19013_001    35163  8159
6  48205    Hartley County, Texas B19013_001    78065 21076
7  48351     Newton County, Texas B19013_001    38871  6573
8  48373       Polk County, Texas B19013_001    57315  2976
9  48139      Ellis County, Texas B19013_001    93248  2485
10 48491 Williamson County, Texas B19013_001   102851  1462
                         geometry
1  MULTIPOLYGON (((-97.3178 27...
2  MULTIPOLYGON (((-97.54085 2...
3  MULTIPOLYGON (((-94.97839 2...
4  MULTIPOLYGON (((-102.5669 3...
5  MULTIPOLYGON (((-105.998 32...
6  MULTIPOLYGON (((-103.0422 3...
7  MULTIPOLYGON (((-93.91113 3...
8  MULTIPOLYGON (((-95.20018 3...
9  MULTIPOLYGON (((-97.08703 3...
10 MULTIPOLYGON (((-98.04989 3...

Interactive viewing with mapview()

library(mapview)

mapview(
  texas_income_sf, 
  zcol = "estimate"
)

US Census Geography

Geography in tidycensus

• Information on available geographies, and how to specify them, can be found in the tidycensus documentation

Searching for variables

• To search for variables, use the load_variables() function along with a year and dataset

• For the 2022 5-year ACS, use "acs5" for the Detailed Tables; "acs5/profile" for the Data Profile; "acs5/subject" for the Subject Tables; and "acs5/cprofile" for the Comparison Profile

• The View() function in RStudio allows for interactive browsing and filtering

vars <- load_variables(2022, "acs5")

View(vars)

Small-area spatial demographic data

• Smaller areas like Census tracts or block groups are available with geography = "tract" or geography = "block group"; one or more counties can optionally be specified to focus your query

nyc_income <- get_acs(
  geography = "tract",
  variables = "B19013_001",
  state = "NY",
  county = c("New York", "Kings", "Queens",
             "Bronx", "Richmond"),
  year = 2022,
  geometry = TRUE
)

“Tidy” or long-form data

• The default data structure returned by tidycensus is “tidy” or long-form data, with variables by geography stacked by row

• For spatial data, this means that geometries will also be stacked, which is helpful for group-wise analysis and visualization

“Tidy” or long-form data

san_diego_race <- get_acs(
  geography = "tract",
  variables = c(
    Hispanic = "DP05_0073P",
    White = "DP05_0079P",
    Black = "DP05_0080P",
    Asian = "DP05_0082P"
  ),
  state = "CA",
  county = "San Diego",
  year = 2022,
  geometry = TRUE
)

“Wide” data

• The argument output = "wide" spreads Census variables across the columns, returning one row per geographic unit and one column per variable

• This will be a more familiar data structure for traditional desktop GIS users

“Wide” data

san_diego_race_wide <- get_acs(
  geography = "tract",
  variables = c(
    Hispanic = "DP05_0073P",
    White = "DP05_0079P",
    Black = "DP05_0080P",
    Asian = "DP05_0082P"
  ),
  state = "CA",
  county = "San Diego",
  geometry = TRUE,
  output = "wide",
  year = 2022
)

Census geometry and the tigris R package

Problem: interior water areas

• Let’s re-visit the NYC income map

• Water areas throughout the city are included within tract polygons, making patterns difficult to observe for those who know the area

• erase_water() in the tigris package offers a solution; it automates the removal of water areas from your shapes

• Tips: tune the area_threshold to selectively remove water area; use cb = FALSE to avoid sliver polygons

• area_threshold = 0.5 means we are removing the largest 50% of water areas in the area

library(tigris)
library(sf)
sf_use_s2(FALSE)

nyc_erase <- erase_water(
  nyc_income_tiger,
  area_threshold = 0.5,
  year = 2022
)

Part 1 exercises

1. Use the load_variables() function to find a variable that interests you that we haven’t used yet.

2. Use get_acs() to fetch spatial ACS data on that variable for a geography and location of your choice, then use mapview() to display your data interactively.

Mapping in R

• R has a robust set of tools for cartographic visualization that make it a suitable alternative to desktop GIS software in many instances

• Popular packages for cartography include ggplot2, tmap, and mapsf

• Today, we’ll be focusing on ggplot2; see Chapter 6 of my book for similar examples using tmap

ggplot2 and geom_sf()

• ggplot2: R’s most popular visualization package (over 137 million downloads!)

• ggplot2 graphics are defined by an aesthetic mapping and one or more geoms

• geom_sf() is a special geom that interprets the geometry type of your spatial data and visualizes it accordingly

• As a result, we can make attractive maps using familiar ggplot2 syntax

Continuous choropleth

• By default, ggplot2 will apply a continuous color palette to create choropleth maps

• Choropleth maps: the shading of a polygon / shape is mapped to a data attribute

Continuous choropleth with styling

• We can apply some styling to customize our choropleth maps

• Used here: a viridis color palette, which is built-in to ggplot2

cont_choro <- ggplot(san_diego_asian, aes(fill = estimate)) + 
  geom_sf() + 
  theme_void() + 
  scale_fill_viridis_c(option = "rocket") + 
  labs(title = "Percent Asian by Census tract",
       subtitle = "San Diego County, CA",
       fill = "ACS estimate",
       caption = "2018-2022 ACS | tidycensus R package")

Classed choropleth

• We can also create a binned choropleth; ggplot2 will identify “pretty” breaks, or custom breaks can be supplied

classed_choro <- ggplot(san_diego_asian, aes(fill = estimate)) + 
  geom_sf() + 
  theme_void() + 
  scale_fill_viridis_b(option = "rocket", n.breaks = 6) + 
  labs(title = "Percent Asian by Census tract",
       subtitle = "San Diego County, CA",
       fill = "ACS estimate",
       caption = "2018-2022 ACS | tidycensus R package")

Faceted choropleth maps

• Spatial ACS data in tidy (long) format can be faceted by a grouping variable, allowing for comparative mapping

faceted_choro <- ggplot(san_diego_race, aes(fill = estimate)) + 
  geom_sf(color = NA) + 
  theme_void() + 
  scale_fill_viridis_c(option = "rocket") + 
  facet_wrap(~variable) + 
  labs(title = "Race / ethnicity by Census tract",
       subtitle = "San Diego County, California",
       fill = "ACS estimate (%)",
       caption = "2018-2022 ACS | tidycensus R package")

Mapping count data

san_diego_race_counts <- get_acs(
  geography = "tract",
  variables = c(
    Hispanic = "DP05_0073",
    White = "DP05_0079",
    Black = "DP05_0080",
    Asian = "DP05_0082"
  ),
  state = "CA",
  county = "San Diego",
  geometry = TRUE,
  year = 2022
)

Graduated symbol maps

• Graduated symbol maps show difference on a map with the size of symbols (often circles)

• They are better for count data than choropleth maps as the shapes are directly comparable (unlike differentially-sized polygons)

• We’ll need to convert our data to centroids to plot graduated symbols in ggplot2

Graduated symbol maps

• We’ll first plot a base layer of Census tracts, then a layer of graduated symbols on top

• Use scale_size_area() to plot proportional symbols

Dot-density mapping

• It can be difficult to show heterogeneity or mixing of different categories on maps

• Dot-density maps scatter dots proportionally to data size; dots can be colored to show mixing of categories

• Traditionally, dot-density maps are slow to make in R; tidycensus’s as_dot_density() function addresses this

Dot-density mapping

san_diego_race_dots <- as_dot_density(
  san_diego_race_counts,
  value = "estimate",
  values_per_dot = 200,
  group = "variable"
)

Dot-density mapping

• Like the graduated symbol map, we plot points over a base layer, but in this case with a much smaller size

• Use override.aes in guide_legend() to plot visible colors in the legend

dot_density_map <- ggplot() + 
  geom_sf(data = san_diego_hispanic, color = "lightgrey", fill = "white") + 
  geom_sf(data = san_diego_race_dots, aes(color = variable), size = 0.01) + 
  scale_color_brewer(palette = "Set1") + 
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  theme_void() + 
  labs(color = "Race / ethnicity",
       caption = "2018-2022 ACS | 1 dot = approximately 200 people")

Adding a basemap

• You may instead want a basemap as a reference layer; this is straightforward with the ggspatial package and annotation_map_tile()

library(ggspatial)

dot_density_with_basemap <- ggplot() + 
  annotation_map_tile(type = "cartolight", zoom = 9) + 
  geom_sf(data = san_diego_race_dots, aes(color = variable), size = 0.01) + 
  scale_color_brewer(palette = "Set1") + 
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  theme_void() + 
  labs(color = "Race / ethnicity",
       caption = "2018-2022 ACS | 1 dot = approximately 200 people")

Customizing interactive maps with mapview()

• mapview() accepts custom color palettes and labels, making it a suitable engine for interactive maps for presentations!

library(viridisLite)

colors <- rocket(n = 100)

mv1 <- mapview(san_diego_asian, zcol = "estimate", 
        layer.name = "Percent Asian<br>2018-2022 ACS",
        col.regions = colors)

Linked interactive maps with mapview()

• In mapview, layers can be stacked with the + operator or swiped between with the | operator

• leafsync takes this one step further by creating side-by-side synced maps

Installing the rdeck package

• To run this example, you’ll need the rdeck package - which is not on CRAN

install.packages("remotes")
library(remotes)
install_github("qfes/rdeck")

Getting migration flow data

• The get_flows() function obtains origin-destination flow data from the 5-year ACS

library(tidyverse)

fulton_inflow <- get_flows(
  geography = "county",
  state = "GA",
  county = "Fulton",
  geometry = TRUE,
  year = 2020
) %>%
  filter(variable == "MOVEDIN") %>%
  na.omit()

fulton_top_origins <- fulton_inflow %>%
  slice_max(estimate, n = 30) 

Mapping in 3D with rdeck!

library(rdeck)

Sys.setenv("MAPBOX_ACCESS_TOKEN" = "pk.eyJ1Ijoia3dhbGtlcnRjdSIsImEiOiJjbHRoYm12eDQwMzZ1MnNvN2JyMzZqYXBpIn0.Sdd16XvAh70IwBrqDD7MzQ")

fulton_top_origins$centroid1 <- st_transform(fulton_top_origins$centroid1, 4326)
fulton_top_origins$centroid2 <- st_transform(fulton_top_origins$centroid2, 4326)

flow_map <- rdeck(map_style = mapbox_light(), 
      initial_view_state = view_state(center = c(-98.422, 38.606), zoom = 3, 
                                      pitch = 45)) %>%
  add_arc_layer(get_source_position = centroid2,
          get_target_position = centroid1,
          data = as_tibble(fulton_top_origins),
          get_source_color = "#274f8f",
          get_target_color = "#274f8f",
          get_height = 1,
          get_width = scale_linear(estimate, range = 1:5),
          great_circle = TRUE
          )

Part 2 exercises

• Try reproducing one of the maps you created in this section with a different state / county combination. What patterns do you observe?

Challenge: making 100 maps at once

Let’s tackle a hypothetical example you might get at work. Your boss has assigned you the following task:

I’d like you to look at geographic patterns in working from home for the 100 largest counties by population in the US. Make me maps for each county.

Related blog post: https://walker-data.com/posts/iterative-mapping/

Warning: may be too memory-intensive for Posit Cloud users

Challenge: making 100 maps at once

At first, this may seem like a significant task! You’ll need to:

Step 1: find the 100 largest counties in the US by population

library(tidycensus)
library(tidyverse)

top100counties <- get_acs(
  geography = "county",
  variables = "B01003_001",
  year = 2022,
  survey = "acs1"
) %>%
  slice_max(estimate, n = 100)

Step 2: pull tract-level Census data by county

• We now need to organize tract data county-by-county using the information in the table. But how can we do this efficiently?

• Something to remember (from Downey, Think Python: )

Repeating identical or similar tasks without making errors is something that computers do well and people do poorly

• My preference for iteration in R: the purrr package

Step 2: pull tract-level data by county

• Here’s how I’d do it:

wfh_tract_list <- top100counties %>%
  split(~NAME) %>%
  map(function(county) {
    state_fips <- str_sub(county$GEOID, 1, 2)
    county_fips <- str_sub(county$GEOID, 3, 5)
    
    get_acs(
      geography = "tract",
      variables = "DP03_0024P",
      state = state_fips,
      county = county_fips,
      year = 2022,
      geometry = TRUE
    )
  })

What just happened?

Now, let’s make some maps:

• We can again use map() to… well… make some maps

library(mapview)

wfh_maps <- map(wfh_tract_list, function(county) {
  mapview(
    county, 
    zcol = "estimate",
    layer.name = "% working from home"
  ) 
})

Prompt: analyzing change over time

Nice work! Now, you have a new follow-up from your boss:

Interesting stuff. However, I’m really interested in knowing how work-from-home is changing over time in our target markets. Could you show me where work-from-home has increased the most in Salt Lake City?

Time-series and the ACS: some notes

• 5-year ACS data represent overlapping survey years, meaning that direct comparison between overlapping datasets is not typically recommended

• The Census Bureau recommends comparing non-overlapping years in the 5-year ACS

• For 2018-2022 data, available 5-year intervals for comparison are 2008-2012 and 2013-2017

The ACS Comparison Profile

• The Comparison Profile includes data on the current 5-year ACS and the 5-year ACS that ends 5 years prior to help with time-series comparisons

utah_wfh_compare <- get_acs(
  geography = "county",
  variables = c(
    wfh17 = "CP03_2017_024",
    wfh22 = "CP03_2022_024"
  ),
  state = "UT",
  year = 2022
)

Geography and making comparisons

• Data in the Comparison Profile tables is only available down to the county level (though counties can change from ACS to ACS!)

• Comparing neighborhood-level change across ACS datasets can introduce additional challenge as Census tract and block group boundaries may differ from previous years

New boundaries drawn with the 2020 Cenus

• Locations with significant population growth (e.g. suburban Collin County, Texas) will have Census tracts / block groups with large populations subdivided in the 2020 Census

• Example: total population by Census tract in Collin County in the 2013-2017 ACS (on the left) and the 2018-2022 ACS (on the right)

Data setup

• Let’s grab data on the population working from home by Census tract in Salt Lake City

• Use a projected coordinate reference system to speed things up!

library(sf)

wfh_17 <- get_acs(geography = "tract", variables = "B08006_017", year = 2017,
                  state = "UT", county = "Salt Lake", geometry = TRUE) %>%
  st_transform(6620)

wfh_22 <- get_acs(geography = "tract", variables = "B08006_017", year = 2022,
                  state = "UT", county = "Salt Lake", geometry = TRUE) %>%
  st_transform(6620)

Method 1: area-weighted interpolation

• Interpolating data between sets of boundaries involves the use of weights to re-distribute data from one geography to another

• A common method is area-weighted interpolation, which allocates information from one geography to another weighted by the area of overlap

• Typically more accurate when going backward, as many new tracts will “roll up” within parent tracts from a previous Census (though not always)

Area-weighted interpolation

library(mapview)
library(leafsync)

m22a <- mapview(wfh_22, zcol = "estimate", layer.name = "2020 geographies")
m17a <- mapview(wfh_22_to_17, zcol = "estimate", layer.name = "2015 geographies")

sync(m22a, m17a)

Method 2: population-weighted interpolation

• If you need to go forward, area-weighted interpolation may be very inaccurate as it can incorrectly allocate large values to low-density / empty areas

• Population-weighted interpolation will instead use an underlying dataset explaining the population distribution as weights

• Example: population-weighted interpolation using block population weights from the 2020 Census

In tidycensus: interpolate_pw()

library(tigris)
options(tigris_use_cache = TRUE)

salt_lake_blocks <- blocks(
  "UT", 
  "Salt Lake", 
  year = 2020
)

wfh_17_to_22 <- interpolate_pw(
  from = wfh_17,
  to = wfh_22,
  to_id = "GEOID",
  weights = salt_lake_blocks,
  weight_column = "POP20",
  crs = 6620,
  extensive = TRUE
)

Evaluating the result

m17b <- mapview(wfh_17, zcol = "estimate", layer.name = "2017 geographies")
m22b <- mapview(wfh_17_to_22, zcol = "estimate", layer.name = "2022 geographies")

sync(m17b, m22b)

Evaluating change over time

• With consistent geographies, we can now join 2017 and 2022 data and evaluate change over time

• We covered this workflow in more detail in the February 22nd workshop

Evaluating change over time

wfh_shift <- wfh_17_to_22 %>%
  select(GEOID, estimate17 = estimate) %>%
  left_join(
    select(st_drop_geometry(wfh_22), 
           GEOID, estimate22 = estimate), by = "GEOID"
  ) %>%
  mutate(
    shift = estimate22 - estimate17,
    pct_shift = 100 * (shift / estimate17)
  )