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

# Write your code here

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

# Write your code here

Take a look at the object you just created.

# Write your code here

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.

# Write your code here

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

# Write your code here

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.

# Write your code here

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:

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.

# Write your code here

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:

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.

# Write your code here

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

# Write your code here

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

# Write your code here

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.

# Write your code here

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.

# Write your code here

What column in ny_counties matches the fips column in countypov?

Answer:

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

# Write your code here

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

# Write your code here

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

Answer:

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

# Write your code here

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.

# Write your code here

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

# Write your code here

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

# Write your code here

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.

# Write your code here

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.

# Write your code here

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:

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

# Write your code here

Which borough looks like it had the most tree planting?

Answer:

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

# Write your code here

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

# Write your code here

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

# Write your code here

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

# Write your code here

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

Answer:

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:

P:

A:

R:

# Write your code here

C:

# Write your code here

K:

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

Answer:

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:

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.

# Write your code here

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

# Write your code here

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!

S:

P:

A:

R:

# Write your code here

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/