Activity 17 — Mapping: Bigfoot, Moose, and Live Maps

Points on a map with sf, real occurrence records from GBIF, and a map you can zoom

maps

spatial

gbif

leaflet

Hands-on companion to the mapping lecture. Students warm up by mapping ~4,000 Bigfoot sightings, then pull real moose occurrence records from GBIF with rgbif, map them across the US including Alaska, and finish with an interactive leaflet map they can pan and zoom. Homework is the same recipe on a species they pick.

Author

Bill Perry

Published

November 10, 2026

In-class Activity 17: Bigfoot, Moose, and Live Maps 🦶🫎

Today’s Objectives

Turn latitude/longitude into a spatial sf object and plot it
Pull real occurrence records from GBIF with rgbif
Map a species across the whole US — Alaska included
Change the scale — the whole world, or one state
Build a live map you can pan and zoom
Do it again on a species you pick, for homework

⚠️ Install these BEFORE class

Two new packages. They pull in a pile of dependencies, so do this at home, not at 9:02 on class day:

install.packages("rgbif")
install.packages("leaflet")

📂 Files for today

In-class skeleton: 17_mapping_skeleton.R → save it as scripts/17_mapping.R. All of today’s code is already in it.
Homework skeleton: 17_mapping_homework.R → save it as scripts/17_mapping_homework.R. This one is what you turn in.
Data: bfro_reports_geocoded.csv in your data/ folder. The moose data downloads itself — you need an internet connection.

What you hand in

Today is a run-it-and-look-at-it day. Every line of code is in this activity and in the in-class skeleton. Run it, read it, change things and see what breaks. You do not turn in scripts/17_mapping.R.
The Q boxes are things to think about. Not graded, not written down.
You turn in scripts/17_mapping_homework.R — a species you pick, mapped with the same recipe. It is a separate file so there is no confusion about what is being graded.

🧩 Chunk 1 — Bigfoot, for fun (after lecture Chunk 1)

Parts 1–3: lat/long becomes a map.

Part 1 · Load up and read the sightings

▶ Run this:

# Load libraries ---------------------------------------
library(tidyverse)
library(sf)       # spatial data frames
library(maps)     # ready-made country and state outlines
library(rgbif)    # download species records from GBIF
library(leaflet)  # live, zoomable maps

▶ Run this:

# Bigfoot reports, trimmed to the lower 48 -------------
bigfoot_df <- read_csv("data/bfro_reports_geocoded.csv") %>%
  filter(!is.na(latitude), !is.na(longitude)) %>%
  filter(longitude > -125, longitude < -66,
         latitude  >   24, latitude  <  50)

nrow(bigfoot_df)

⚠️ Watch out! Those four numbers are a bounding box — a crude rectangle around the lower 48. It throws Alaska and Hawaii away on purpose. Remember that: in Part 7 we will want Alaska back, and this is the line that would have eaten it.

Part 2 · Make it spatial

🔮 Predict first: you are about to plot ~4,000 points with no basemap at all. 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)

✏️ Q1: In coords = c("longitude", "latitude"), which comes first — x or y? What happens if you swap them?

✏️ Q2: The points drew a recognizable USA with no map underneath. What does that tell you about what a “map” actually is?

💡 Key idea: an sf object is a data frame with one extra column — geometry. Every dplyr verb you already know still works on it.

Part 3 · Put it on a real map

▶ Run this:

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

# Basemap first, points on top -------------------------
ggplot() +
  geom_sf(data = states_sf, fill = "grey96", color = "grey75") +
  geom_sf(data = bigfoot_sf, alpha = 0.2, size = 0.5, color = "firebrick") +
  theme_void()

✏️ Q3: Why does the order of the two geom_sf() layers matter — what would you see if you drew them the other way around?

✏️ Q4: Is this map about Bigfoot, or about people — where they live, hike, and file reports? What would tell those apart?

That last question is the whole reason we start here. Bigfoot is a joke dataset with a serious lesson in it: an occurrence map shows you where someone looked, not only where the thing is. Hold onto that — the moose data has exactly the same problem.

🧩 Chunk 2 — Real occurrence data: GBIF (after lecture Chunk 2)

Parts 4–5: download actual species records.

Part 4 · Meet GBIF

GBIF — the Global Biodiversity Information Facility — pools occurrence records from museums, herbaria, survey programs, and apps like iNaturalist. Billions of records, free, and the rgbif package queries it straight from R.

Today’s species is the moose, Alces alces. You need the scientific name — GBIF does not take “moose.”

▶ Run this:

# Up to 2000 US moose records that have coordinates ----
moose_raw <- occ_data(
  scientificName = "Alces alces",
  country        = "US",
  hasCoordinate  = TRUE,
  limit          = 2000
)$data

dim(moose_raw)

scientificName — genus and species, spelled exactly
country = "US" — a two-letter country code
hasCoordinate = TRUE — skip records with no lat/long; they are useless to us
limit — how many records to bring back
$data on the end — occ_data() hands back a list, and the records are in its data slot

⚠️ Watch out! 119 columns arrive. GBIF records carry everything anyone ever recorded — collector, institution, elevation, licence. We want five of them.

Part 5 · What did `limit` actually give you?

▶ Run this:

# How many US moose records does GBIF actually have? ---
occ_count(scientificName = "Alces alces", country = "US")

# And what years did OUR 2000 records come from? -------
moose_raw %>%
  count(year, sort = TRUE)

✏️ Q5: GBIF holds about 28,000 US moose records. You asked for 2,000 — and they all came from one year. So what does limit actually do: take a random sample, or take the first 2,000 GBIF hands back?

✏️ Q6: If you mapped these 2,000 and called it “the range of moose in the US,” what would be wrong with that claim?

💡 Key idea: limit is not a sample, it is a cutoff. Every occurrence map you will ever make is a map of sampling effort as much as of the species. Say so out loud when you present one.

🧩 Chunk 3 — Map the moose, Alaska included (after lecture Chunk 3)

Parts 6–8: clean it, map the whole country, then change the scale.

Part 6 · Trim and make it spatial

▶ Run this:

# Five columns is plenty -------------------------------
moose_df <- moose_raw %>%
  select(species, stateProvince, decimalLatitude,
         decimalLongitude, year) %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude))

moose_df %>%
  count(stateProvince, sort = TRUE) %>%
  head()

▶ Run this:

# Same st_as_sf() as Bigfoot, different column names ---
moose_sf <- moose_df %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"),
           crs = 4326)

moose_sf

GBIF calls its coordinates decimalLongitude and decimalLatitude. The recipe did not change — only the names you hand to coords. Longitude still goes first.

✏️ Q7: Which state has the most moose records? Does that surprise you?

Part 7 · A basemap with Alaska in it

Our states_sf from Part 3 is the lower 48 only — it has no Alaska. The single biggest cluster of moose records in the country is in Alaska, so we need a different basemap.

▶ Run this:

# Country outlines, Alaska included --------------------
usa_sf <- st_as_sf(
  maps::map("world", regions = "USA", plot = FALSE, fill = TRUE)
) %>%
  st_set_crs(4326)

# Moose across the whole country -----------------------
ggplot() +
  geom_sf(data = usa_sf, fill = "grey96", color = "grey75") +
  geom_sf(data = moose_sf, alpha = 0.3, size = 0.7, color = "darkgreen") +
  coord_sf(xlim = c(-170, -60), ylim = c(25, 72)) +
  theme_void()

The coord_sf() limits are doing real work here. Without them the map stretches from Maine to the far end of the Aleutian Islands — which cross the international date line — and the whole country squashes into a sliver.

✏️ Q8: Where are the moose? Name the regions. Does that match what you know about the species — and about where it is cold?

✏️ Q9: There are no moose in the southern US on this map. Is that a real range limit, or a sampling artifact like the Bigfoot map? How is this case different?

Part 8 · Change the scale — the whole world, or one state

Every map so far has been “the US.” Changing the scale means turning two knobs, and they are separate:

Knob	What it controls	Where you turn it
The download	which records GBIF sends you	the arguments to `occ_data()`
The basemap	which outlines get drawn underneath	the `maps::map()` call (and `coord_sf()` limits)

Turn one without the other and you get world-wide points on a US map, or US points floating on an empty planet.

8a · The whole world

Knob 1: just leave country out. No country means everywhere.

▶ Run this:

# Moose anywhere on Earth ------------------------------
moose_world_raw <- occ_data(
  scientificName = "Alces alces",
  hasCoordinate  = TRUE,
  limit          = 2000
)$data

moose_world_raw %>%
  count(country, sort = TRUE) %>%
  head()

Knob 2: maps::map("world") with no regions = is every country.

▶ Run this:

# Same recipe: drop missing coords, make it spatial ----
moose_world_sf <- moose_world_raw %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"),
           crs = 4326)

# Every country on Earth -------------------------------
world_sf <- st_as_sf(
  maps::map("world", plot = FALSE, fill = TRUE)
) %>%
  st_set_crs(4326)

ggplot() +
  geom_sf(data = world_sf, fill = "grey96", color = "grey75") +
  geom_sf(data = moose_world_sf, alpha = 0.3, size = 0.7, color = "darkgreen") +
  theme_void()

🌍 Want an actual globe? Add one more line to the end of that plot: + coord_sf(crs = "+proj=ortho +lat_0=60 +lon_0=-30"). That is an orthographic projection — the Earth as seen from space. lat_0 and lon_0 are the spot you are looking straight down on; change them and the globe turns.

✏️ Q10: Most of these 2,000 records come from Norway and Sweden. Are there more moose in Scandinavia than in Canada and Russia — or more people logging them? (Remember what limit does.)

8b · One state

Knob 1: add stateProvince. Spell it the way GBIF stores it — full name, capitalized.

▶ Run this:

# Moose in Minnesota only ------------------------------
moose_mn_raw <- occ_data(
  scientificName = "Alces alces",
  country        = "US",
  stateProvince  = "Minnesota",
  hasCoordinate  = TRUE,
  limit          = 2000
)$data

nrow(moose_mn_raw)

moose_mn_raw %>%
  count(year) %>%
  arrange(year)

Knob 2: maps::map("county", regions = "minnesota") draws the state with its counties. This time the name is lower case — the maps package spells things its own way.

▶ Run this:

# Same recipe again ------------------------------------
moose_mn_sf <- moose_mn_raw %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"),
           crs = 4326)

# One state, with its counties -------------------------
mn_sf <- st_as_sf(
  maps::map("county", regions = "minnesota", plot = FALSE, fill = TRUE)
) %>%
  st_set_crs(4326)

ggplot() +
  geom_sf(data = mn_sf, fill = "grey96", color = "grey75") +
  geom_sf(data = moose_mn_sf, alpha = 0.4, size = 1, color = "darkgreen") +
  theme_void()

⚠️ Watch out! Two spellings of the same state: "Minnesota" for GBIF, "minnesota" for maps. Want a different state? Change both. Want state outlines only, no counties? Use maps::map("state", regions = "minnesota", ...).

✏️ Q11: The US download was all one year. The Minnesota download goes back to the 1980s — with the same limit = 2000. Why? (Hint: compare nrow() to the limit.)

✏️ Q12: Find Duluth on the map. What is the big cluster up the North Shore telling you — about moose, and about who is out there looking?

💡 Key idea: scale is two decisions — what you download (country, stateProvince, or neither) and what you draw underneath ("world", "state", "county"). A small area often gets you the whole record; a big one gets you only the newest slice.

🧩 Chunk 4 — Go live

(after lecture Chunk 4)

Parts 9–10: a map you can actually zoom.

Part 9 · Your first live map

Everything so far has been a picture. leaflet makes a real web map — tiles, pan, zoom, scroll — right inside RStudio’s Viewer pane.

▶ Run this:

# A map you can zoom -----------------------------------
leaflet(moose_sf) %>%
  addTiles() %>%
  addCircleMarkers(radius = 4, stroke = FALSE,
                   fillOpacity = 0.6, color = "darkgreen")

Three lines:

leaflet(moose_sf) — start a map from your sf object
addTiles() — put the actual world map underneath, streamed from the web
addCircleMarkers() — drop a dot on every record

Zoom in. Scroll onto Anchorage, onto the Colorado Rockies, onto Maine. This is the same data as Part 7 — but now you can ask it questions by looking.

✏️ Q13: Zoom all the way in on a cluster. Do the points land in forest and wetland, or do they line up along roads and towns? What does that tell you about how these records were collected?

Part 10 · Labels on the dots

▶ Run this:

# Click any dot to see what it is ----------------------
leaflet(moose_sf) %>%
  addTiles() %>%
  addCircleMarkers(
    radius      = 4,
    stroke      = FALSE,
    fillOpacity = 0.6,
    color       = "darkgreen",
    popup       = ~paste0(species, "<br>", stateProvince, "<br>", year)
  )

The ~ makes popup a formula — “look these columns up in the data for each row.” <br> is a line break, because the popup is a little web page.

✏️ Q14: Click a few dots. Is every record the species you asked for?

💡 Key idea: the recipe never changed. Get coordinates → st_as_sf() → draw them. Bigfoot, moose, or anything you download next semester — same four steps, different file.

Part 11 · If there is time — sightings per state

This is the left_join() from last class, pointed at a map. Counts are one table, polygons are another, and state is the key.

▶ Run this:

# Count Bigfoot sightings per state --------------------
sightings_by_state <- bigfoot_df %>%
  mutate(state = str_to_lower(state)) %>%
  count(state, name = "sightings")

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

▶ 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) +
  coord_sf(crs = 5070) +
  theme_void()

direction = -1 flips the scale so more sightings = darker; the default reads backwards to most people. coord_sf(crs = 5070) is an equal-area projection for the lower 48 — which is why we did not use it on the moose map, where Alaska would have been mangled.

✏️ Q15: Exactly one polygon comes back NA and draws grey. Is that a broken join or a true zero — and how would you tell?

✏️ Q16: Washington and California are big, populous, and forested. What extra data would make this a fair state-to-state comparison?

Homework — map a species of your own

Open 17_mapping_homework.R. That file — not today’s script — is what you turn in.

You pick one species you actually care about (or several, if you want to compare them) and run today’s recipe on it:

Look up its scientific name — GBIF will not take a common name
occ_data() it from GBIF
Check what you actually got — how many records, and from when
st_as_sf() and map it on the US basemap
Change the scale — map it on the whole world, or in one state
Make it live with leaflet
Say whether the map matches the range you expected, and where sampling effort is showing through

Everything you need is in the homework skeleton, with the moose code sitting right there as the worked example to copy from. All the graded questions are in one block at the top of the file — answer them there as # comments, so they are all in one place.

Review and checkpoint

At this point you should be able to:

Build an sf object with st_as_sf(coords = ..., crs = 4326)
Draw a basemap + points with geom_sf() — basemap first
Pull records for any species from GBIF with occ_data()
Say what limit does, and why an occurrence map is a sampling map
Use a basemap that includes Alaska, and crop it with coord_sf()
Change the scale: leave out country for the world, add stateProvince for one state, and match the basemap
Make a live, zoomable map with leaflet()

💡 Key idea: you started with a joke dataset and finished by downloading real biodiversity records for any species on Earth and putting them on a map you can zoom. That is the same twenty lines of code either way.

End of Activity 17. Next: your final project pitches.

Other Formats

Activity 17 — Mapping: Bigfoot, Moose, and Live Maps

In-class Activity 17: Bigfoot, Moose, and Live Maps 🦶🫎

Today’s Objectives

🧩 Chunk 1 — Bigfoot, for fun (after lecture Chunk 1)

Part 1 · Load up and read the sightings

Part 2 · Make it spatial

Part 3 · Put it on a real map

🧩 Chunk 2 — Real occurrence data: GBIF (after lecture Chunk 2)

Part 4 · Meet GBIF

Part 5 · What did `limit` actually give you?

🧩 Chunk 3 — Map the moose, Alaska included (after lecture Chunk 3)

Part 6 · Trim and make it spatial

Part 7 · A basemap with Alaska in it

Part 8 · Change the scale — the whole world, or one state

8a · The whole world

8b · One state

🧩 Chunk 4 — Go live

Part 9 · Your first live map

Part 10 · Labels on the dots

Part 11 · If there is time — sightings per state

Homework — map a species of your own

Review and checkpoint

Other Formats

In-class Activity 17: Bigfoot, Moose, and Live Maps 🦶🫎

Today’s Objectives

🧩 Chunk 1 — Bigfoot, for fun (after lecture Chunk 1)

Part 1 · Load up and read the sightings

Part 2 · Make it spatial

Part 3 · Put it on a real map

🧩 Chunk 2 — Real occurrence data: GBIF (after lecture Chunk 2)

Part 4 · Meet GBIF

Part 5 · What did limit actually give you?

🧩 Chunk 3 — Map the moose, Alaska included (after lecture Chunk 3)

Part 6 · Trim and make it spatial

Part 7 · A basemap with Alaska in it

Part 8 · Change the scale — the whole world, or one state

8a · The whole world

8b · One state

🧩 Chunk 4 — Go live

Part 9 · Your first live map

Part 10 · Labels on the dots

Part 11 · If there is time — sightings per state

Homework — map a species of your own

Review and checkpoint

Part 5 · What did `limit` actually give you?