Skip to main content

Lecture — Mapping: Bigfoot, Moose, and Live Maps

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

maps
sf
spatial
gbif
leaflet
ggplot

An intro to maps in R. Spatial data and coordinate systems, the sf package, and geom_sf() — warming up on ~5,000 Bigfoot sighting reports, then downloading real moose occurrence records from GBIF with rgbif, mapping them across the US including Alaska, and finishing with a live leaflet map.

Author

Bill Perry

Published

November 10, 2026

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:

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)

Scatter of Bigfoot sighting coordinates with no basemap, which already resembles the outline of the United States.

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

Map of the lower 48 US states with Bigfoot sighting points overlaid in red, denser in the Pacific Northwest.

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

Map of the United States including Alaska with moose occurrence points in dark green, clustered in Alaska, the northern Rockies, the upper Midwest, and New England.

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 map with moose occurrence points in dark green across northern North America, Scandinavia, and Russia, densest in Norway and Sweden.

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

Minnesota county map with moose occurrence points clustered in the northeast Arrowhead region along the North Shore of Lake Superior, with a few scattered in the northwest.

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

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.