---
title: "Week 5 Tutorial"
subtitle: "Geospatial Data and Maps"
format:
  html:
    toc: true
    toc-location: left
    toc-depth: 3
    highlight-style: github
    code-copy: true
    embed-resources: true
    df-print: paged
    include-in-header:
      text: |
        <style>
          .cell-output pre { max-height: 250px; overflow-y: auto; }
        </style>
editor: visual
---

First things first: Make a new project in RStudio, *and be sure to place it somewhere where you can easily find it* (e.g., Documents, Desktop). Then download this file (the `.qmd` file) from the syllabus website and place it in the folder associated with that project. Then open this file in RStudio and follow along.

I encourage you to work with/talk to others throughout this tutorial, especially during times where I have you code on your own and use AI. When you code, think about the SPARCK framework:

**S: Set up** (What is your end goal?)

**P: Pseudocode** (What are the steps?)

**A: AI** (What are your AI prompts?)

**R: Run** (Run it!)

**C: Check** (Did your code work?)

**K: Knowledge** (What does your code mean?)

## Agenda

- What makes data "spatial"?

- Shape files and `tigris`

- Adding data to shape files

- Counts vs. rates

- Polishing your maps

- Practice

Learning how to work with spatial data could be an entire class in and of itself. It's a vast field, variously referred to as GIS (geographic information systems), geospatial mapping, etc. We don't have time in this class to learn all of it. Instead, you're going to learn how to make some cool visualizations of location-based data, along with the **core concepts** behind them, such as the difference between points and polygons, what a shape file is, or what a choropleth is. AI is very good at writing mapping code, but it can only help you if you know what to ask for, and you can only check its work if you know what a correct map should look like.

If you want to learn more about working with spatial data in R, a great (free) place to start is *Geocomputation with R* (Lovelace et al., 2025), which you can find [here](https://r.geocompx.org/).

## What makes data "spatial"?

Spatial data is data where **location** is part of the information. Every row is tied to some place on Earth. Most spatial data you'll encounter describes places using three basic kinds of shapes, or **geometries** (Lovelace et al., 2025):

- **Points**: a single location. (E.g., the address of a school.)

- **Lines**: a series of connected points. (E.g., a street, a subway line, or a river.)

- **Polygons**: a series of connected points that closes back on itself to create an area. (E.g., a state, a county, or a police precinct.)

The locations that make up these shapes are usually described with **longitude** (how far east or west, the x axis) and **latitude** (how far north or south, the y axis).

## Shape files and `tigris`

Shape files are a special file type used in GIS. They contain information about the locations of borders and objects, and are the basis of drawing *geographic shapes*. You don't need to know exactly how they work, because an R package called `sf` (which stands for **simple features**) can read them for you (Pebesma, 2018). `sf` is the main package that people use for working with spatial data and building maps in R.

Now, where do we *get* these shape files? There are many sources online. But the most important source of shape files for the US is the Census Bureau, which publishes the official boundaries for just about every kind of geographic area in the country. The `tigris` package lets you download them directly into R (Walker, 2023).

If you haven't already, install the packages we'll need today in the console. (Remember: install packages in the console, not in your Quarto document.) We'll use `usmap` later on, just for some of the datasets that come with it.

```{r}
#| eval: false
install.packages(c('sf', 'tigris', 'usmap'))
```

Load `tidyverse`, `sf`, and `tigris`. Then, set the option `tigris_use_cache` to `TRUE`. This tells `tigris` to save the files it downloads, so it doesn't have to download the same files every time you run your code. (Hint: Use the `options()` function.)

```{r}
library(tidyverse)
library(sf)
library(tigris)
options(tigris_use_cache = TRUE)
```

Use the `counties()` function from `tigris` to download all of the counties in New York State. Set the argument `cb = TRUE`, which downloads a simplified version of the boundaries that is smaller and faster to work with. (We'll use `cb = TRUE` for most of today.)

```{r}
ny_counties <- counties(state = 'NY', cb = TRUE)
```

Take a look at the object you just created.

```{r}
ny_counties
```

The **big idea** of the `sf` package is that a spatial data frame is just **a regular data frame with a geometry column**. That means that all of the Data Moves we learned (filtering, calculating, joining, etc.) work exactly the same way on spatial data. The geometry column just comes along for the ride.

### Making a map with `geom_sf()`

Last week we learned that every `ggplot` graph needs `data`, a `geom`, and `aes`thetics. Maps are no different! `sf` comes with its own geom, `geom_sf()`, which knows how to draw whatever shapes are in the geometry column. You don't need to specify x and y aesthetics, because the geometry column already contains the locations.

Use `ggplot()` and `geom_sf()` to make a map of the counties in New York State.

```{r}
ggplot(data = ny_counties) +
  geom_sf()
```

Notice the axes are longitude and latitude. `ggplot` figured that out from the geometry column. Usually we don't need axes on a map, though. Add a theme that gets rid of them. (Hint: Last week we learned about `theme_` functions. Which one would get rid of *everything* except the map itself?)

```{r}
ggplot(data = ny_counties) +
  geom_sf() +
  theme_void()
```

### FIPS codes

Let's zoom in on New York City. To do this, we need a way to identify the five counties that make up NYC (one for each borough). We *could* use county names, but names are messy: The Bronx or Bronx? Is Manhattan called Manhattan or New York County? (It's New York County.)

Instead, spatial data usually uses **FIPS codes**, which are standardized ID numbers that the federal government assigns to every state and county (and smaller areas too). A county's FIPS code has 5 digits: 2 for the state (New York is 36), and 3 for the county. It's good practice to use FIPS codes whenever you're working with spatial data, and you can easily look them up on Google. In `tigris` data, the FIPS code is in the `GEOID` column.

The FIPS codes for the five NYC counties are `36005` (Bronx), `36047` (Kings/Brooklyn), `36061` (New York/Manhattan), `36081` (Queens), and `36085` (Richmond/Staten Island). Use `filter()` to create a new object that contains only the five NYC counties, and map it.

```{r}
nyc_counties <- ny_counties |>
  filter(GEOID %in% c('36005', '36047', '36061', '36081', '36085'))

ggplot(data = nyc_counties) +
  geom_sf() +
  theme_void()
```

Notice that we just used `filter()` like we would with any data frame!

### Exploring `tigris`

`tigris` can download a lot more than counties. You can see the full list of what it can download [here](https://github.com/walkerke/tigris). It can download three main kinds of shape files:

1.  **Census geographies**: the areas the Census uses to report its data (states, counties, census tracts, etc.)

2.  **Political and administrative boundaries**: congressional districts, state legislative districts, school districts, cities and towns, etc.

3.  **Physical features**: roads, railroads, rivers, lakes, coastlines, etc.

Every `tigris` function works in a similar way: the name of the function is the kind of area you want (e.g., `counties()`), and the arguments let you narrow it down (e.g., `state =`, `county =`, `year =`). Once you've learned a few, you can use any of them.

#### Census geographies

Census geographies are **nested**, like Russian dolls. Each level fits perfectly inside the level above it:

| Level        | Roughly how big?                         | `tigris` function |
|--------------|------------------------------------------|-------------------|
| State        | (You know this one)                      | `states()`        |
| County       | NYC has 5, one for each borough          | `counties()`      |
| Census tract | About 4,000 people (a neighborhood)      | `tracts()`        |
| Block group  | About 600 to 3,000 people (a few blocks) | `block_groups()`  |
| Block        | Often a single city block                | `blocks()`        |

Because the areas are nested, their FIPS codes are too: a county's code starts with its state's code, a tract's code starts with its county's code, and so on.

Use `tracts()` to download the census tracts for Manhattan (New York County). Then map them.

```{r}
manhattan_tracts <- tracts(state = 'NY', county = 'New York', cb = TRUE)

ggplot(data = manhattan_tracts) +
  geom_sf() +
  theme_void()
```

Look at the `GEOID` column for the tracts. How many digits does each code have? The first five digits should look familiar. What do they mean?

**Answer:** Tract codes have 11 digits. The first 5 digits are the FIPS code for New York County (`36` for New York State, and `061` for New York County). The last 6 digits identify the tract within the county.

Why does this matter? The Census (and lots of other sources) report data **at specific levels**, like "median income by census tract." To put that data on a map, you need to download shapes at the **same level**, and match them up using the FIPS code.

#### Political and administrative boundaries

Census geographies are useful for data, but they're not the boundaries that determine who **represents** you, where you go to school, or which government provides your services. `tigris` can download those too. A few examples:

- `congressional_districts()`: districts for the US House of Representatives

- `state_legislative_districts()`: districts for the state legislature (use `house = 'upper'` for the State Senate and `house = 'lower'` for the State Assembly)

- `school_districts()`: school districts

- `places()`: cities, towns, and villages

Use `congressional_districts()` to download New York's congressional districts, and map them.

```{r}
ny_cd <- congressional_districts(state = 'NY', cb = TRUE)

ggplot(data = ny_cd) +
  geom_sf() +
  theme_void()
```

Let's zoom in on the districts in NYC. Use `st_filter()` to keep only the districts that overlap with `nyc_counties`. Then make a map of those districts, with the outlines of the five boroughs on top. (Hint: `ggplot` works in **layers**, so you can add two `geom_sf()` layers, each with its own `data =`. Give the boroughs no fill, using `fill = NA`, and a different border color, so you can see the districts underneath.)

```{r}
nyc_cd <- ny_cd |>
  st_filter(nyc_counties)

ggplot() +
  geom_sf(data = nyc_cd) +
  geom_sf(data = nyc_counties, fill = NA, color = 'red') +
  theme_void()
```

Layers are drawn in order, so whatever comes last is drawn on top. We'll use this a lot today.

You'll notice the congressional districts don't fit neatly inside the boroughs, the way census tracts fit inside counties. Political boundaries don't have to follow census geographies. This matters if you want to combine data from both. For example, you can't just add up borough-level data to get data for a congressional district.

#### Physical features

`tigris` can also download physical features of the landscape. Notice that these are often **lines**, not polygons. A few examples:

- `roads()`: every road in a county (`primary_roads()` gets just the major highways)

- `area_water()`: lakes, rivers, bays, etc.

- `rails()`: railroads

Use `roads()` to download all of the roads in Manhattan. Then make a map with Manhattan's census tracts on the bottom and the roads on top. (Hint: There are a *lot* of roads, so you'll probably want to make the lines thin, using the `linewidth =` argument.)

```{r}
manhattan_roads <- roads(state = 'NY', county = 'New York')

ggplot() +
  geom_sf(data = manhattan_tracts, fill = 'grey90') +
  geom_sf(data = manhattan_roads, linewidth = .1) +
  theme_void()
```

Look closely at the edges of the island. Some tracts stick out into the water, past where the roads end. That's because census tracts often include the water around them!

`tigris` has a function to fix this: `erase_water()`, which cuts the water out of any shape file. Use it to make a new version of the Manhattan tracts without water, and remake the map.

```{r}
manhattan_tracts_land <- erase_water(manhattan_tracts)

ggplot() +
  geom_sf(data = manhattan_tracts_land, fill = 'grey90') +
  geom_sf(data = manhattan_roads, linewidth = .1) +
  theme_void()
```

Finally, `tigris` data comes with some useful columns besides the shapes. Look at `nyc_counties` again. `ALAND` is the land area and `AWATER` is the water area (both in square meters). They could be useful for calculating things like population density (people per square mile). We'll use `ALAND` later on.

## Adding data to shape files

### Choropleth maps

Maps of boundaries are nice, but what's better is plotting **data** on the map. One of the most common types of maps is a **choropleth map**, in which each area is filled in with a color that represents some value (e.g., poverty rate by county, or election results by state).

To make a choropleth map, you need two things: (1) the shapes of your areas, which we can get from `tigris`, and (2) the data you want to show for each area. Usually these come from different places, so we need to **join** them together, using a column they have in common.

The `usmap` package comes with a few useful datasets. One of them is called `countypov`, which has the percentage of people living in poverty in every US county in 2021. Load `usmap` and take a look at `countypov`.

```{r}
library(usmap)
countypov
```

What column in `ny_counties` matches the `fips` column in `countypov`?

**Answer:** The `GEOID` column. Both contain FIPS codes.

Use `left_join()` to join `countypov` to `ny_counties`. Since the matching columns have *different names* in each data frame, you'll need to tell it which columns to join by. (Make sure `ny_counties` comes first in your join, so that the result is still a spatial data frame!)

```{r}
ny_counties_pov <- ny_counties |>
  left_join(countypov, by = join_by(GEOID == fips))
```

Now make a choropleth map of New York State, where the fill color of each county represents its poverty rate. (Hint: Just like last week, if a visual feature like fill color **conveys information**, it goes inside `aes()`.)

```{r}
ggplot(data = ny_counties_pov) +
  geom_sf(aes(fill = pct_pov_2021)) +
  theme_void()
```

Which county has the highest poverty rate? Was it easy to spot on the map?

**Answer:** The Bronx (26.4%). But it's so small on the map that it's hard to see! This is a common problem with choropleth maps: big areas (like rural counties) take up most of the map, even if very few people live there, and small areas (like cities) are easy to miss, even if lots of people live there.

Now make the same kind of map, but just for the five NYC counties. (Hint: You can either filter `ny_counties_pov`, or join `countypov` to `nyc_counties`.)

```{r}
nyc_counties_pov <- nyc_counties |>
  left_join(countypov, by = join_by(GEOID == fips))

ggplot(data = nyc_counties_pov) +
  geom_sf(aes(fill = pct_pov_2021)) +
  theme_void()
```

### Points

Let's now plot specific **points** on our map. The file `2019_schools.csv` (download it from the course website and put it in the same folder as this tutorial) comes from NYC Open Data, and contains location data for all public schools in NYC in the 2019-2020 school year. Load it in, clean the column names, and take a look at it.

```{r}
library(janitor)
schools <- read_csv('2019_schools.csv') |>
  clean_names()
```

To put the schools on our map, we need to turn `schools` into a spatial data frame. We do this with `st_as_sf()`. (Most `sf` functions start with `st_`, which stands for "spatial type.") You'll need two arguments: `coords = c('longitude', 'latitude')`, which tells it which columns have the coordinates (longitude always comes first!), and `crs = 4326`, which tells R that the coordinates are regular longitude and latitude. Use AI to help you create a spatial version of `schools`. (Hint: `st_as_sf()` won't work if any rows are missing coordinates, so you'll need to filter those out first.)

```{r}
schools_sf <- schools |>
  filter(!is.na(longitude), !is.na(latitude)) |>
  st_as_sf(coords = c('longitude', 'latitude'), crs = 4326)
```

Now plot the schools on top of your poverty map, by adding another `geom_sf()` layer.

```{r}
ggplot() +
  geom_sf(data = nyc_counties_pov, aes(fill = pct_pov_2021)) +
  geom_sf(data = schools_sf) +
  theme_void()
```

Ooh la la! Let's change a few visual features to make it look nicer. Make the dots smaller, white, and partially transparent (using `alpha =`), so we can see where schools are more and less dense.

```{r}
ggplot() +
  geom_sf(data = nyc_counties_pov, aes(fill = pct_pov_2021)) +
  geom_sf(data = schools_sf, size = .5, color = 'white', alpha = .3) +
  theme_void()
```

### Lines: reading in shape files

So far, we've accessed our shape files using `tigris`. But shape files are just regular files that can be downloaded from anywhere (like NYC Open Data), and read in using the `st_read()` function from `sf`.

Let's try that with a shape file from NYC Open Data that has information on trees that were planted as part of an initiative called MillionTreesNYC. Download `nyc_trees.zip` from the course website, unzip it, and put the `nyc_trees` folder that's inside it in your project folder.

Look inside the folder. A shape file is actually a *bundle* of several files with the same name and different extensions (`.shp`, `.shx`, `.dbf`, `.prj`, etc.). Each file holds a different piece of the puzzle (the shapes, the data, etc.), and they **must be kept together** in the same folder. To read in the shape file, use `st_read()` with the folder name as the only argument.

```{r}
nyc_trees <- st_read('nyc_trees')
```

Look at the message that printed when you read it in, and take a look at the data. What type of geometry is this data? What do you think each row represents?

**Answer:** The geometry type is `LINESTRING`, so each row is a line. Each row is a stretch of street (a block) where trees were planted.

Plot the trees on top of your poverty map, in green. (You can drop the schools for now.)

```{r}
ggplot() +
  geom_sf(data = nyc_counties_pov, aes(fill = pct_pov_2021)) +
  geom_sf(data = nyc_trees, color = 'lightgreen') +
  theme_void()
```

Which borough looks like it had the most tree planting?

**Answer:** It's a bit hard to tell, but the Bronx and Queens look like they have a lot of tree planting. Manhattan looks like it has less.

The `plnt_seas` column tells us the season and year each tree was planted (e.g., "Fall11" means fall 2011). Use AI to help you create a new column with just the year of planting. (Hint: The `str_sub()` function can grab the last two characters of a piece of text.) Then make a map of the five boroughs where the color of each tree planting site shows the year it was planted. (Hint: You'll probably want to make the lines a bit thinner.)

```{r}
nyc_trees <- nyc_trees |>
  mutate(plant_year = paste0('20', str_sub(plnt_seas, -2)))

ggplot() +
  geom_sf(data = nyc_counties) +
  geom_sf(data = nyc_trees, aes(color = plant_year), linewidth = .3) +
  theme_void()
```

## Counts vs. rates

It's hard to compare boroughs by eyeballing thousands of lines. Let's count the number of tree planting sites in each borough, and make a choropleth map of the counts.

The `boroughcod` column tells us which borough each tree planting site is in. First, use `count()` to count the number of tree planting sites for each value of `boroughcod`. (Hint: Use `st_drop_geometry()` on the trees data before counting. This turns it back into a regular data frame. We don't need the lines anymore, and it makes counting much faster.)

```{r}
tree_counts <- nyc_trees |>
  st_drop_geometry() |>
  count(boroughcod, name = 'num_trees')

tree_counts
```

To join these counts to our shape file, we need a column with each borough's FIPS code. NYC's borough codes are `1` = Manhattan (`36061`), `2` = Bronx (`36005`), `3` = Brooklyn (`36047`), `4` = Queens (`36081`), and `5` = Staten Island (`36085`). Add a `GEOID` column to `tree_counts` with the right FIPS code for each borough. (Hint: `case_when()` is useful for this.)

```{r}
tree_counts <- tree_counts |>
  mutate(GEOID = case_when(boroughcod == 1 ~ '36061',
                           boroughcod == 2 ~ '36005',
                           boroughcod == 3 ~ '36047',
                           boroughcod == 4 ~ '36081',
                           boroughcod == 5 ~ '36085'))
```

Now join `tree_counts` to `nyc_counties_pov`, and make a choropleth map of the number of tree planting sites.

```{r}
nyc_trees_map <- nyc_counties_pov |>
  left_join(tree_counts, by = 'GEOID')

ggplot(data = nyc_trees_map) +
  geom_sf(aes(fill = num_trees)) +
  theme_void()
```

Which borough had the most tree planting? Which had the least?

**Answer:** The Bronx had the most (5,533), closely followed by Queens (5,257). Manhattan had the least (2,200).

But wait. Is this a fair comparison? The boroughs are very different sizes. Queens is almost five times as big as Manhattan, so of course it has room for more trees! This is one of the most important principles of choropleth maps: a map of **counts** (raw numbers) often ends up just being **a map of which areas are biggest** (or have the most people). Big places will have more of almost everything.

If we want to compare places fairly, we usually need to map a **rate** instead: a count divided by some measure of size. For tree planting, the best measure of size is probably land area. Let's calculate the number of tree planting sites **per square mile** in each borough, using SPARCK.

**S:** Make a choropleth map of the number of tree planting sites per square mile in each borough.

**P:** Convert the `ALAND` column (square meters) into square miles. Create a new column that divides the number of tree planting sites by the number of square miles. Map the new column.

**A:** I have an sf data frame in R called `nyc_trees_map`, with one row for each NYC borough. It has a column called `num_trees` with the number of tree planting sites, and a column called `ALAND` with the land area in square meters. Using tidyverse, create a new column with the land area in square miles, and another new column called `trees_per_sq_mi` with the number of tree planting sites per square mile. Then make a choropleth map of `trees_per_sq_mi` using ggplot and geom_sf, with theme_void.

**R:**

```{r}
nyc_trees_map <- nyc_trees_map |>
  mutate(sq_mi = ALAND / 2589988,
         trees_per_sq_mi = num_trees / sq_mi)

ggplot(data = nyc_trees_map) +
  geom_sf(aes(fill = trees_per_sq_mi)) +
  theme_void()
```

**C:** I looked at the `sq_mi` column to make sure the areas make sense. Manhattan is about 23 square miles and Queens is about 109, which matches what I found when I Googled the size of each borough. The rates are between about 50 and 130 tree planting sites per square mile.

```{r}
nyc_trees_map |>
  st_drop_geometry() |>
  select(NAME, num_trees, sq_mi, trees_per_sq_mi)
```

**K:** `mutate()` created two new columns. First, it converted land area from square meters into square miles (there are about 2,589,988 square meters in a square mile). Then it divided the number of tree planting sites in each borough by the borough's area. The Bronx had the most tree planting per square mile (about 131), followed by Manhattan (about 97). Queens had the least (about 48).

Compare the count map to the rate map. What changed? Which map tells a more accurate story?

**Answer:** In the count map, Queens looks like it got almost as much tree planting as the Bronx, and Manhattan looks like it got the least. But once we account for size, Queens drops to last place, and Manhattan jumps to second. The rate map is a fairer comparison, because it accounts for the fact that some boroughs have much more land than others.

Now compare the rate map to the poverty map. Do you notice any pattern? What questions would you want to ask about how the city decided where to plant trees?

**Answer:** The Bronx has both the highest poverty rate and the most tree planting per square mile. We might want to know whether the program deliberately targeted lower-income neighborhoods or neighborhoods that had fewer trees to begin with (which is often the same thing). With only five boroughs, though, we can't say much. We'd want to look at smaller areas, like census tracts, to see whether the pattern holds.

## Polishing your maps

All of the `ggplot` polishing tools we learned last week work on maps too. A few that are especially useful for maps:

- **Color scales.** `scale_fill_viridis_c()` and `scale_fill_distiller()` give you color scales that are easy to read and colorblind-friendly. A *sequential* scale (light to dark) is best for values that go from low to high, like rates. A *diverging* scale (two colors with a neutral middle) is best when there's a meaningful midpoint, like change from last year.

- **Borders.** The `color =` and `linewidth =` arguments in `geom_sf()` control the borders between areas. (`color = NA` removes them entirely, which can be useful when you have lots of small areas.)

- **Labels.** `labs()` lets you add a title, subtitle, and caption, and set the legend title (using the name of the aesthetic, e.g., `fill =`). It's good practice to include a caption with the source of your data.

- **Legend position.** `theme(legend.position = ...)` lets you move the legend (e.g., `'right'`, `'bottom'`, or `'none'` to get rid of it), which can help if it's covering up part of your map.

- **Text labels.** `geom_sf_text()` and `geom_sf_label()` add text on top of each shape (like the name of each area), using the `label =` aesthetic.

Make a polished version of your trees per square mile map. Use a green color scale, white borders between boroughs, a title, a clear legend title, a caption with your data source, and a bigger font.

```{r}
ggplot(data = nyc_trees_map) +
  geom_sf(aes(fill = trees_per_sq_mi), color = 'white', linewidth = .5) +
  scale_fill_distiller(palette = 'Greens', direction = 1) +
  labs(title = 'MillionTreesNYC tree planting by borough',
       fill = 'Planting sites\nper square mile',
       caption = 'Data: NYC Open Data; US Census Bureau') +
  theme_void(base_size = 14)
```

Like last week, you don't need to memorize all of these options. But you should know what's possible, so that you can **describe to AI exactly what you'd like your map to look like**.

### Saving your maps

Save your map with `ggsave()`, just like any other `ggplot` graph. Maps often need a bit more space, so it can help to specify the width and height (in inches).

```{r}
ggsave('trees_map.jpg', width = 7, height = 7)
```

## Practice

Now it's your turn. Make your own map below, using SPARCK. It must:

- Use at least one `tigris` function that we **didn't** use in this tutorial (you can find the full list [here](https://github.com/walkerke/tigris)), **or** a `tigris` function we did use, but for a different place or year.

- Show some kind of data, not just boundaries.

For example, you could map poverty rates by county in a state of your choice, map population density (people per square mile) in each county of a state, or show the major highways (`primary_roads()`) on top of a choropleth map. **Show us something interesting!**

**S:**

**P:**

**A:**

**R:**

```{r}

```

**C:**

**K:**

## References

Lovelace, R., Nowosad, J., & Muenchow, J. (2025). *Geocomputation with R* (2nd ed.). CRC Press. <https://r.geocompx.org/>

Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. *The R Journal, 10*(1), 439–446. <https://doi.org/10.32614/RJ-2018-009>

Walker, K. (2023). *Analyzing US Census data: Methods, maps, and models in R*. CRC Press. <https://walker-data.com/census-r/>