Turn latitude/longitude into maps with sf and geom_sf, then build a per-state choropleth
2026-07-05
Download bfro_reports_geocoded.csv from Kaggle (https://www.kaggle.com/datasets/thedevastator/unlocking-mysteries-of-bigfoot-through-sightings) and save it into your project’s data/ folder. It has ~4,586 sightings with latitude, longitude, state, season, and classification.
sf objectHow to use this worksheet
- Work through the parts in order. Type the code into a new R script and run it line by line.
- Blocks marked ▶ Run this should be executed as written.
- Blocks marked ✏️ Your turn ask you to write, modify, or interpret.
- Day 1 = Parts 1–6. Day 2 = Parts 7–12.
🔮 Predict before you run — and type, don’t paste
Before each map draws, predict where the points or the dark states will be. Then type the code yourself. Predicting the pattern (“Pacific Northwest?”) is what makes a map mean something when it appears — and typing builds the sf + geom_sf habits.
🧩 Chunk 1 — Make it spatial (Day 1, after lecture Chunk 1)
Parts 1–3: load the sightings, clean coordinates, and build an
sfobject.
▶ Run this:
library(tidyverse)
library(sf) # spatial data frames
library(maps) # ready-made state outlines
# Swap this path for your download -----------------------
bigfoot_df <- read_csv("data/bfro_reports_geocoded.csv")
bigfoot_df %>%
select(title, state, latitude, longitude, classification, season) %>%
head()✏️ Your turn: How many rows and columns does the raw file have? Run dim(bigfoot_df).
Rows: Columns:
What does one row represent?▶ Run this:
✏️ Your turn: How many rows remain? Why must we drop missing coordinates before the next step?
Rows remaining:
Why drop NA coords first:sf object🔮 Predict first: You’re about to plot 4,000+ points with NO basemap. What shape will the cloud of points make?
▶ Run this:
✏️ Your turn: In coords = c("longitude", "latitude"), which comes first — x or y? What did the bare points end up looking like?
First coordinate is (x / y):
The points look like:🧩 Chunk 2 — Points on a real map (Day 1, after lecture Chunk 2)
Parts 4–6: add a US basemap and layer the sightings on top.
sf▶ Run this:
✏️ Your turn: What is in the ID column, and what case are the state names in? (You’ll need this for the join later.)
ID column holds:
Case (UPPER / lower / Title):🔮 Predict first: Which region will have the densest cluster of points? Write your guess.
▶ Run this:
# Basemap first, points on top ---------------------------
ggplot() +
geom_sf(data = states_sf, fill = "grey96", color = "grey75") +
geom_sf(data = bigfoot_sf, aes(color = classification),
alpha = 0.3, size = 0.6) +
labs(title = "Bigfoot Sightings Across the Lower 48",
color = "Report class") +
theme_void()✏️ Your turn: Was your prediction right? Why does the order of the two geom_sf() layers matter?
Densest region:
Why layer order matters:✏️ Your turn: Change theme_void() to theme_minimal(). What comes back? Which looks better for a map, and why?
What theme_minimal() adds:
Better for a map and why:🛑 End of Day 1. You can put points on a map. Day 2: the choropleth.
🧩 Chunk 3 — Choropleth by state (Day 2, after lecture Chunk 3)
Parts 7–9: count per state, join to the map, shade it.
🔮 Predict first: Which state has the most Bigfoot reports? Commit before you run
count().
▶ Run this:
✏️ Your turn: Record the top three, and why we lowercase the state names now.
Top 3 states:
Why str_to_lower():▶ Run this:
✏️ Your turn: What would happen to the map if you did NOT lowercase the state names before joining? (Try it — join count(bigfoot_df, state) directly.)
Your prediction / what happened:▶ Run this:
✏️ Your turn: Does the darkest state match your Part 7 prediction? Name one state that surprised you.
Darkest state:
A surprise:🧩 Chunk 4 — Polish & explore (Day 2, after lecture Chunk 4)
Parts 10–12: reproject, explore by season, and save.
▶ Run this:
✏️ Your turn: Compare this to the Part 9 map. What changed about the shape of the country?
Your answer:🔮 Predict first: Which season has the most reported sightings? Why might that be (think about who is outside)?
▶ Run this:
# One small map per season -------------------------------
bigfoot_sf %>%
filter(season %in% c("Spring","Summer","Fall","Winter")) %>%
ggplot() +
geom_sf(data = states_sf, fill = "grey96", color = "grey80") +
geom_sf(alpha = 0.2, size = 0.4, color = "firebrick") +
facet_wrap(~ season) +
coord_sf(crs = 5070) +
theme_void()✏️ Your turn: Which season dominates? Is that about Bigfoot, or about people?
Busiest season:
What it really tells you:▶ Run this:
choropleth_plot <- ggplot(states_map) +
geom_sf(aes(fill = sightings), color = "white") +
scale_fill_viridis_c(option = "magma", direction = -1) +
coord_sf(crs = 5070) +
labs(title = "Bigfoot Sightings by State", fill = "Sightings") +
theme_void()
ggsave("figures/bigfoot_choropleth.png",
plot = choropleth_plot,
width = 7, height = 5, units = "in", dpi = 300)You should now be able to:
✏️ Your turn — before you move on: Run your whole script top to bottom. Does it run cleanly?
Ran cleanly? Y / N
If not, what error appeared:Optional — do this if you finish early.
▶ Try this: map only the clearest (“Class A”) sightings.
✏️ Your turn: Does restricting to Class A change the geographic pattern?
Your answer:✏️ Your turn: Washington and California have many sightings — but they’re also big, populous states. What extra data would you need to make this a fair comparison between states?
Your answer:st_as_sf error about missing coordinates → filter out NA lat/long first (Part 2).st_set_crs(4326).anti_join().could not find function "st_as_sf" → install.packages("sf"), then library(sf).💡 Key idea: every map is the same recipe — get geometry, get data, join them, and shade or plot. You just did it with Bigfoot; the homework does it with real species from GBIF.
End of Worksheet 13. Homework: map your own species from GBIF with rgbif.