Worksheet — Mapping Bigfoot Sightings

Turn latitude/longitude into maps with sf and geom_sf, then build a per-state choropleth

maps

spatial

Hands-on companion to the mapping lecture. Over two days students load ~4,600 Bigfoot sightings, make them spatial with sf, plot points on a US map, and join sighting counts to state polygons to build a choropleth.

Author

Bill Perry

Published

September 10, 2026

Mapping Bigfoot 🦶

Before you start — get the data

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,600 sightings with latitude, longitude, state, season, and classification.

Today’s Objectives (a 2-day worksheet)

Make latitude/longitude into a spatial sf object
Draw sighting points on a map of US states
Count sightings per state and join them to the map
Shade a choropleth and reproject it for a proper US look

How 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, name the region or state you expect to stand out. Then type the code yourself — st_as_sf(), geom_sf(), and coord_sf() all take arguments in a specific order (coords = c("longitude", "latitude"), crs = 4326), and typing them is what keeps that order from slipping the next time you write it from scratch.

🧩 Chunk 1 — Make it spatial (Day 1, after lecture Chunk 1)

Parts 1–3: load the sightings, clean coordinates, and build an sf object.

Part 1 · Load libraries and data

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

Part 2 · Clean the coordinates

▶ Run this:

# Drop missing coords + keep the lower 48 ----------------
bigfoot_df <- bigfoot_df %>%
  filter(!is.na(latitude), !is.na(longitude)) %>%
  filter(longitude > -125, longitude < -66,
         latitude  >   24, latitude  <  50)

nrow(bigfoot_df)

✏️ Your turn: How many rows remain? Why must we drop missing coordinates before the next step?

Rows remaining:
Why drop NA coords first:

Part 3 · Build the `sf` object

🔮 Predict first: You’re about to plot ~4,600 points with NO basemap. What shape will the cloud of points make?

▶ Run this:

# Turn lon/lat into a geometry column --------------------
bigfoot_sf <- bigfoot_df %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

ggplot(bigfoot_sf) +
  geom_sf(alpha = 0.2, size = 0.5) +
  theme_minimal()

✏️ 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.

Part 4 · Get US states as `sf`

▶ Run this:

# Built-in state outlines -> sf --------------------------
states_sf <- st_as_sf(
  maps::map("state", plot = FALSE, fill = TRUE)
) %>%
  st_set_crs(4326)

states_sf %>% head(3)

✏️ 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):

Part 5 · Layer the sightings on the map

🔮 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:

Part 6 · Make it yours

✏️ Your turn: Change theme_void() to theme_minimal(). What comes back? Which looks better for a map, and why?

# Write your modified plot here:

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.

Part 7 · Count sightings per state

🔮 Predict first: Which state has the most Bigfoot reports? Commit before you run count().

▶ Run this:

# Sightings per state, lowercased to match the map -------
sightings_by_state <- bigfoot_df %>%
  mutate(state = str_to_lower(state)) %>%
  count(state, name = "sightings") %>%
  arrange(desc(sightings))

head(sightings_by_state)

✏️ Your turn: Record the top three, and why we lowercase the state names now.

Top 3 states:
Why str_to_lower():

Part 8 · Join counts to the map

▶ Run this:

# Join sightings onto the state polygons -----------------
states_map <- states_sf %>%
  mutate(state = str_remove(ID, ":.*")) %>%   # drop ":north" etc.
  left_join(sightings_by_state, by = "state")

# Peek — every state should now have a sightings number
states_map %>% st_drop_geometry() %>% head()

✏️ 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:

Part 9 · Shade the choropleth

▶ Run this:

# The map IS the bar chart -------------------------------
ggplot(states_map) +
  geom_sf(aes(fill = sightings), color = "white") +
  scale_fill_viridis_c(option = "magma", direction = -1) +
  labs(title = "Bigfoot Sightings by State", fill = "Sightings") +
  theme_void()

✏️ 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.

Part 10 · Reproject for a proper US look

▶ Run this:

# Albers Equal-Area (EPSG:5070) --------------------------
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 (Albers)", fill = "Sightings") +
  theme_void()

✏️ Your turn: Compare this to the Part 9 map. What changed about the shape of the country?

Your answer:

Part 11 · Explore by season

🔮 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:

Part 12 · Save your map

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

Part 13 · Review and checkpoint

You should now be able to:

Explain points, polygons, and CRS
Build an sf object with st_as_sf(coords = ..., crs = 4326)
Draw a basemap + points with geom_sf() (basemap first!)
Count per group and left_join() the counts to polygons
Shade a choropleth with geom_sf(aes(fill = ...))
Reproject with coord_sf(crs = 5070)

✏️ 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:

Part 14 · Going further

Optional — do this if you finish early.

Class A only

▶ Try this: map only the clearest (“Class A”) sightings.

bigfoot_sf %>%
  filter(classification == "Class A") %>%
  ggplot() +
  geom_sf(data = states_sf, fill = "grey96", color = "grey80") +
  geom_sf(alpha = 0.3, size = 0.5, color = "darkgreen") +
  coord_sf(crs = 5070) +
  theme_void()

✏️ Your turn: Does restricting to Class A change the geographic pattern?

Your answer:

Per-capita thinking

✏️ 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:

Extension — out of class (~30–40 min)

Turn this in with your worksheet. In class you mapped Bigfoot reports. Now you map a real species you pick, pulled live from GBIF with rgbif.

Important

E2 and E3 must be handwritten on paper, photographed, and embedded (![caption](my_photo.jpg)). Typed answers get at most half credit, even if correct — E2 is your prediction before you download.

E1 · Map your own species (4 pts)

Pick any species with a scientific name (a plant, bird, mammal, insect — your choice). Then, with the same recipe as class:

library(rgbif)
library(sf)
library(tidyverse)

# 1. up to ~1500 US records that have coordinates
occ <- occ_data(scientificName = "Genus species",
                country = "US", hasCoordinate = TRUE, limit = 1500)$data

# 2. keep real coordinates, build the sf object (x = longitude first!)
sp_sf <- occ %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

# 3. points on the same states basemap you used in class
ggplot() +
  geom_sf(data = states_sf, fill = "grey96", color = "grey80") +
  geom_sf(data = sp_sf, alpha = 0.3, size = 0.6, color = "darkgreen") +
  coord_sf(crs = 5070) +
  theme_void()

Save the map to figures/ at dpi = 300. Record your species and how many records you got.

E2 · Predict, then check — ✍️ by hand (3 pts)

Before running E1: sketch, on a rough US outline, where you expect your species to occur and why (climate, habitat, range you know). Then run E1 and write 2–3 sentences comparing your prediction to the real map — and note whether anything looks like a sampling artifact rather than real range (a dense blob around one city, a hard line at a state border).

E3 · Explain it — ✍️ by hand (3 pts)

GBIF records are where people reported the species, not where it truly lives. Give one concrete way that could bias your map.
Washington and California had many Bigfoot reports partly because they are big and populous. What extra data would make a state-to-state comparison of your species fair?
The choropleth recipe (polygons → counts → left_join() → geom_sf(fill=)) didn’t change from Bigfoot to your species. In your own words, why is that recipe dataset-independent?

Getting unstuck

st_as_sf error about missing coordinates → filter out NA lat/long first (Part 2).
Points and states don’t line up → both need the same CRS; set the states to st_set_crs(4326).
Whole choropleth is grey → your join keys don’t match (case!). Lowercase both sides.
Empty states in the choropleth → a state name spelled differently in the two tables; check with anti_join().
could not find function "st_as_sf" → install.packages("sf"), then library(sf).
Cheat sheet — https://r-spatial.github.io/sf/

💡 Key idea: the choropleth recipe doesn’t change with the dataset — get polygons, get counts, left_join() them on a matching key, then shade with geom_sf(). You just ran it on Bigfoot sightings; the Extension runs the same four steps on a species you pick from GBIF.

End of the Mapping worksheet.

Other Formats

Worksheet — Mapping Bigfoot Sightings

Mapping Bigfoot 🦶

Before you start — get the data

Today’s Objectives (a 2-day worksheet)

How to use this worksheet

🧩 Chunk 1 — Make it spatial (Day 1, after lecture Chunk 1)

Part 1 · Load libraries and data

Part 2 · Clean the coordinates

Part 3 · Build the `sf` object

🧩 Chunk 2 — Points on a real map (Day 1, after lecture Chunk 2)

Part 4 · Get US states as `sf`

Part 5 · Layer the sightings on the map

Part 6 · Make it yours

🧩 Chunk 3 — Choropleth by state (Day 2, after lecture Chunk 3)

Part 7 · Count sightings per state

Part 8 · Join counts to the map

Part 9 · Shade the choropleth

🧩 Chunk 4 — Polish & explore (Day 2, after lecture Chunk 4)

Part 10 · Reproject for a proper US look

Part 11 · Explore by season

Part 12 · Save your map

Part 13 · Review and checkpoint

Part 14 · Going further

Class A only

Per-capita thinking

Extension — out of class (~30–40 min)

E1 · Map your own species (4 pts)

E2 · Predict, then check — ✍️ by hand (3 pts)

E3 · Explain it — ✍️ by hand (3 pts)

Getting unstuck

Other Formats

Mapping Bigfoot 🦶

Before you start — get the data

Today’s Objectives (a 2-day worksheet)

How to use this worksheet

🧩 Chunk 1 — Make it spatial (Day 1, after lecture Chunk 1)

Part 1 · Load libraries and data

Part 2 · Clean the coordinates

Part 3 · Build the sf object

🧩 Chunk 2 — Points on a real map (Day 1, after lecture Chunk 2)

Part 4 · Get US states as sf

Part 5 · Layer the sightings on the map

Part 6 · Make it yours

🧩 Chunk 3 — Choropleth by state (Day 2, after lecture Chunk 3)

Part 7 · Count sightings per state

Part 8 · Join counts to the map

Part 9 · Shade the choropleth

🧩 Chunk 4 — Polish & explore (Day 2, after lecture Chunk 4)

Part 10 · Reproject for a proper US look

Part 11 · Explore by season

Part 12 · Save your map

Part 13 · Review and checkpoint

Part 14 · Going further

Class A only

Per-capita thinking

Extension — out of class (~30–40 min)

E1 · Map your own species (4 pts)

E2 · Predict, then check — ✍️ by hand (3 pts)

E3 · Explain it — ✍️ by hand (3 pts)

Getting unstuck

Part 3 · Build the `sf` object

Part 4 · Get US states as `sf`