Lecture — Mapping: Bigfoot, Moose, and Live Maps

Spatial data and sf, real species records from GBIF, and an interactive leaflet map

Bill Perry

2026-11-10

Where we left off — Joins

• Joins — combined two fisher tables by a shared key, and used anti_join() to find what didn’t match
• Keys had to be cleaned before they would match — case, type, format
• That same skill comes back at the end of today, pointed at a map

Note

✅ Transition

Today has two halves. First, Bigfoot — a joke dataset with a serious lesson in it. Then the real thing: we download actual occurrence records for a species from GBIF, map them across the country, and finish with a map you can pan and zoom. 🦶🫎

Goals for today

• Understand spatial data — points, polygons, and coordinate systems
• Meet sf — a data frame with a geometry column
• Turn latitude/longitude into a map with geom_sf()
• Pull real species records from GBIF with rgbif
• Map a species across the US — Alaska included
• Change the scale — the whole world, or one state
• Build a live, zoomable map with leaflet

Tools today:

• tidyverse, sf, maps
• rgbif, leaflet — install

Data:

• bfro_reports_geocoded.csv — in your data/ folder
• Moose records — downloaded live from GBIF

Reference:

• 🌐 rgbif
• 🌐 leaflet for R

Tip

Four chunks: Bigfoot → GBIF → the moose map → going live.

How to Use These Slides

This lecture runs in four short chunks. After each chunk you switch to the activity and run the code yourself.

Today is different: every line is already written, in the activity and in the skeleton. Your job is to run it, read it.

Important

⚠️ Install first

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

Do this at start of class — they pull in a lot of dependencies.

Note

The homework is a separate file

17_mapping_homework.R. That one you turn in — a species you pick.

🧩 Chunk 1 of 4 · Bigfoot, for Fun

We will cover: what spatial data is, and how latitude/longitude becomes a map.

What Is Spatial Data?

Three ideas and you have it:

• Points — one coordinate pair. A sighting, a tree, a weather station.
• Polygons — a closed ring of coordinates. A state, a lake, a park.
• CRS — the coordinate reference system: what the numbers mean.

4326 is plain latitude/longitude in degrees — what a phone GPS gives you, and what almost every dataset arrives in.

Note

📖 Vocabulary

Two datasets will not line up unless they share a CRS. Most mismatched-map bugs are two layers in two different systems.

The sf Package — a Data Frame with Geometry

# Load packages + the sighting data ------------------
library(tidyverse)
library(sf)       # spatial data frames
library(maps)     # country + state outlines
library(rgbif)    # species records from GBIF
library(leaflet)  # live maps

# The file is in your data/ folder ------------------
bigfoot_df <- read_csv("data/bfro_reports_geocoded.csv")

dim(bigfoot_df)

[1] 5021   29

• An sf object is a tibble with one extra column, geometry
• Every dplyr verb you know still works on it
• ~5,000 real BFRO reports, each geocoded to a lat/long

Clean the Coordinates First

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

[1] 4031

• st_as_sf() errors on missing coordinates — filter them first
• The four numbers are a bounding box around the lower 48

Warning

⚠️ Remember this line. It throws Alaska away on purpose. In Chunk 3 we want Alaska back.

Make It Spatial — st_as_sf()

Note

🔮 Predict first: ~4,000 rows of lat/long, plotted with no basemap. What shape will the cloud make?

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

# Just the points — the USA appears on its own -------
ggplot(bigfoot_sf) +
  geom_sf(alpha = 0.2, size = 0.5)

• coords = c("longitude", "latitude") — x first, always
• crs = 4326 — tell sf these are degrees
• No map needed: the data drew the country

Put It on a Real Map

# Basemap first, then points on top ------------------
states_sf <- st_as_sf(
  maps::map("state", plot = FALSE, fill = TRUE)
) %>%
  st_set_crs(4326)

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

• Order matters — states first, points on top
• theme_void() drops axes and grid

Note

🤔 The real question

Is this a map of Bigfoot, or a map of where people live, hike, and file reports? Hold that thought — the moose data has the same problem.

🛑 Pause — Do Activity Parts 1–3 Now

Read the sightings, make them spatial, and put them on the map. Then ask yourself what the map is really showing.

🧩 Chunk 2 of 4 · Real Occurrence Data — GBIF

We will cover: where real species records come from, and how to pull them into R.

Tip

🖐 After this chunk: Activity Parts 4–5 (download moose, then check what you got).

Where Real Occurrence Data Lives

GBIF — the Global Biodiversity Information Facility — pools occurrence records from:

• museum and herbarium collections
• government survey programs
• citizen science apps like iNaturalist

Billions of records. Free. And rgbif queries it straight from R — no website, no download button.

Warning

⚠️ Scientific names only

GBIF does not know “moose.” It knows Alces alces. Common names are ambiguous across languages and regions; the genus species is not.

Downloading Moose 🫎

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

dim(moose_raw)

[1] 2000  119

• scientificName — genus + species, spelled exactly
• country = "US" — two-letter code
• hasCoordinate = TRUE — skip records we can’t map
• limit — how many to bring back
• $data — occ_data() returns a list; the records are in its data slot

119 columns arrive. We want five.

What Did limit Actually Give You?

# How many US moose records EXIST? -------------------
occ_count(scientificName = "Alces alces",
          country = "US")

[1] 28020

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

# A tibble: 1 × 2
   year     n
  <int> <int>
1  2026  2000

About 28,000 records exist. We asked for 2,000 — and they all came from one year.

Important

limit is a cutoff, not a sample.

It hands back the first n records GBIF returns, newest first. Map those and call it “the range” and you are making a claim your data cannot support.

🛑 Pause — Do Activity Parts 4–5 Now

Download the moose records, then find out how many exist and what years you actually got.

🧩 Chunk 3 of 4 · Map the Moose — Alaska Included

We will cover: the same st_as_sf() recipe on new column names, a basemap that doesn’t stop at the Canadian border, and how to zoom out to the world or in to one state.

Tip

🖐 After this chunk: Activity Parts 6–8 (trim it, map the whole country, then change the scale).

Trim It and Make It Spatial

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

moose_sf <- moose_df %>%
  st_as_sf(coords = c("decimalLongitude",
                      "decimalLatitude"),
           crs = 4326)

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

# A tibble: 4 × 2
  stateProvince     n
  <chr>         <int>
1 Alaska          641
2 Colorado        469
3 Wyoming         261
4 Idaho           107

GBIF names its coordinates decimalLongitude / decimalLatitude.

The recipe did not change — only the names you hand to coords. Longitude still goes first.

Alaska has more moose records than anywhere else. Which is a problem…

The Lower-48 Basemap Problem

states_sf from Chunk 1 is the lower 48 only. It has no Alaska at all — so the biggest cluster of moose in the country would simply not be drawn.

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

Warning

⚠️ Two traps at once

1. The Aleutians cross the date line, so the raw extent runs −178 to +180 and squashes the country to a sliver. coord_sf() limits fix it.
2. crs = 5070 is CONUS Albers — it mangles Alaska. Leave this map unprojected.

Moose Across the 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()

Read the biology off it:

• Alaska, the northern Rockies, the upper Midwest, New England
• Nothing in the south — moose are a cold-adapted animal with a real southern limit

Unlike Bigfoot, this pattern is mostly the species, not only the observers. Mostly.

Change the Scale — Two Knobs

# The whole world: leave country OUT -----------------
moose_world_raw <- occ_data(
  scientificName = "Alces alces",
  hasCoordinate  = TRUE,
  limit          = 2000
)$data

# One state: ADD stateProvince -----------------------
moose_mn_raw <- occ_data(
  scientificName = "Alces alces",
  country        = "US",
  stateProvince  = "Minnesota",
  hasCoordinate  = TRUE,
  limit          = 2000
)$data

nrow(moose_world_raw)

[1] 2000

nrow(moose_mn_raw)

[1] 481

Scale is two separate decisions:

1. The download — what occ_data() asks for: country, stateProvince, or neither
2. The basemap — what maps::map() draws: "world", "state", "county"

Warning

⚠️ Turn both knobs

World points on a US basemap, or Minnesota points on the whole planet — the map is technically right and useless.

Knob 2 — The Whole World

moose_world_sf <- moose_world_raw %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"),
           crs = 4326)

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.5, color = "darkgreen") +
  theme_void()

• "world" with no regions = → every country
• Over half the 2,000 are Norway + Sweden — more moose, or more people logging them?

🌍 Want a globe? Add one line:

+ coord_sf(crs = "+proj=ortho
             +lat_0=60 +lon_0=-30")

lat_0 / lon_0 = the spot you look straight down on.

Knob 2 — One State, With Counties

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

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 = 0.8, color = "darkgreen") +
  theme_void()

• Two spellings: "Minnesota" for GBIF, "minnesota" for maps — change both to switch states
• The North Shore cluster — moose habitat, and Highway 61

Important

A small area gets you the whole record.

The US pull was all one year. Minnesota has only a few hundred records — under the limit — so you get them all, back to the 1980s.

🛑 Pause — Do Activity Parts 6–8 Now

Trim the records, map them on a basemap with Alaska in it, then change the scale — the whole world, and one state.

🧩 Chunk 4 of 4 · Go Live 🗺️

We will cover: turning a static picture into a map you can pan, zoom, and click.

Tip

🖐 After this chunk: Activity Parts 9–10 (make it live, then add popups).

leaflet — a Real Web Map

# 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 from your sf object
• addTiles() — the actual world map, streamed from the web
• addCircleMarkers() — a dot per record

Important

▶ Run this one yourself

It opens in the Viewer pane. It can’t be printed onto a slide — the whole point is that it moves.

Zoom in: Anchorage. The Colorado Rockies. Maine.

Click a Dot — Popups

leaflet(moose_sf) %>%
  addTiles() %>%
  addCircleMarkers(
    radius      = 4,
    stroke      = FALSE,
    fillOpacity = 0.6,
    color       = "darkgreen",
    popup       = ~paste0(species, "<br>",
                          stateProvince, "<br>", year)
  )

• ~ makes popup a formula — look these columns up for each row
• <br> is a line break, because a popup is a tiny web page

Note

🤔 Zoom all the way in

Do the points sit in habitat — or do they trace roads, towns, and trailheads? That is sampling effort, visible with your own eyes.

🛑 Pause — Do Activity Parts 9–10 Now

Make the map live, add popups, and zoom into a cluster until you can see what the records are really tracking.

What We Learned Today

Concepts:

• Spatial data = points/polygons + a CRS
• sf objects are tibbles with a geometry column
• geom_sf() draws layers; order matters (basemap first)
• GBIF holds real occurrence records for nearly any species
• limit is a cutoff, not a sample — every occurrence map is partly a map of sampling effort
• A lower-48 basemap silently drops Alaska
• Scale = two knobs: what you download and what you draw underneath

R skills:

• st_as_sf(), st_set_crs(), maps::map(), coord_sf()
• occ_data(), occ_count() from rgbif
• leaflet(), addTiles(), addCircleMarkers()

References:

• 🌐 rgbif
• 🌐 leaflet for R
• 🌐 sf package

Important

Homework — a separate file

17_mapping_homework.R. Pick a species you care about, run today’s recipe on it, and make it live. The moose code is sitting right there to copy from.

Up next: your final project pitches.