Week 5 Tutorial

Geospatial Data and Maps

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.

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.

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.)

library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.4.3

Warning: package 'tidyr' was built under R version 4.4.3

Warning: package 'purrr' was built under R version 4.4.3

Warning: package 'dplyr' was built under R version 4.4.3

-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
v dplyr     1.2.0     v readr     2.1.5
v forcats   1.0.0     v stringr   1.5.1
v ggplot2   4.0.2     v tibble    3.2.1
v lubridate 1.9.4     v tidyr     1.3.2
v purrr     1.2.1     
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

Warning: package 'sf' was built under R version 4.4.3

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(tigris)

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

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.)

ny_counties <- counties(state = 'NY', cb = TRUE)

Retrieving data for the year 2024

Take a look at the object you just created.

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 aesthetics. 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.

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?)

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.

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. It can download three main kinds of shape files:

Census geographies: the areas the Census uses to report its data (states, counties, census tracts, etc.)
Political and administrative boundaries: congressional districts, state legislative districts, school districts, cities and towns, etc.
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.

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

Retrieving data for the year 2024

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.

ny_cd <- congressional_districts(state = 'NY', cb = TRUE)

Retrieving data for the year 2024

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.)

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.)

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.

manhattan_tracts_land <- erase_water(manhattan_tracts)

Fetching area water data for your dataset's location...

Erasing water area...
If this is slow, try a larger area threshold value.

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.

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!)

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().)

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.)

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.

library(janitor)


Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

schools <- read_csv('2019_schools.csv') |>
  clean_names()

Rows: 2190 Columns: 41

-- Column specification --------------------------------------------------------
Delimiter: ","
chr  (32): system_code, location_code, location_name, Managed_by_name, locat...
dbl   (8): fiscal_year, BEDS, X_COORDINATE, Y_COORDINATE, LONGITUDE, LATITUD...
dttm  (1): open_date

i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.

Warning in FUN(X[[i]], ...): unable to translate '<U+00C4>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00D6>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00DC>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00E4>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00F6>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00FC>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00DF>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00C6>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00E6>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00D8>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00F8>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00C5>' to native encoding

Warning in FUN(X[[i]], ...): unable to translate '<U+00E5>' to native encoding

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.)

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.

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.

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.

nyc_trees <- st_read('nyc_trees')

Reading layer `geo_export_abdecf61-f1a2-4fa2-b666-4a223b6a0134' from data source `/Users/benjaminsilver/Library/CloudStorage/GoogleDrive-bms2202@columbia.edu/My Drive/Postdoc/Teaching/RM_Fall2026/Coding/Tutorials/Week5_GeospatialData/nyc_trees' 
  using driver `ESRI Shapefile'
Simple feature collection with 20513 features and 8 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -74.22339 ymin: 40.52271 xmax: -73.71178 ymax: 40.91247
Geodetic CRS:  WGS 84

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.)

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.)

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.)

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.)

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.

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.

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.

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.

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).

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), 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!

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/