Activity 12 — Real Climate Data in R

How fast is your city warming? (regression on time)

tidyverse
descriptive-stats

Hands-on companion to the Weather lecture: download NOAA station data with GSODR, clean the names, summarize to yearly means, and fit lm(temp ~ year) to get a warming rate — for Duluth, for each season, and for a city you pick.

Author

Bill Perry

Published

October 5, 2026

In-class Activity 12: Real Climate Data

Today’s Objectives

  1. Download real NOAA weather station data from inside R with GSODR
  2. Clean the column names with janitor::clean_names()
  3. Summarize ~28,000 daily records into one mean per year
  4. Fit lm(temp ~ year) and read the warming rate off summary()
  5. Split the year into seasons with case_when() and compare warming rates
  6. Find the station for a city you choose and run it again
Tip📂 Files for today
  • Skeleton script: 12_weather_skeleton.R → save it as scripts/12_weather.R.
  • No data file today — you download the data from NOAA inside R.
  • New package: install GSODR once, from the Console: install.packages("GSODR")
ImportantWhat you hand in
  • The code is already written in the skeleton — run it chunk by chunk as you work through Parts 1–10 below.
  • What gets graded is your answers, typed into the Q boxes in the script, plus the Extension at the bottom (your own city).
  • Hand in scripts/12_weather.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.

Answer boxes look like this — type after ANSWER: and start every line with #:

# # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Q5. What is the slope, with units?
# ANSWER: 0.027 degrees C per year
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Keep the plots plain today — no labels, colors, or themes. The goal is to read them.

🧩 Chunk 2 — Duluth’s warming rate (after lecture Chunk 2)

Parts 1–5: download, trim, plot, summarize, and fit the trend.

Part 1 · Load libraries and download the data

▶ Run this (already at the top of your skeleton):

# Load libraries ---------------------------------------
library(tidyverse)   # dplyr + ggplot2
library(janitor)     # clean_names()
library(GSODR)       # NOAA weather station data

# Download daily Duluth weather, 1948-2025 -------------
duluth_df <- get_GSOD(years = 1948:2025,
                      station = "727450-14913") %>%
  clean_names()

This takes a minute — it is pulling 78 years of daily records from NOAA.

▶ Run this:

# How big is it, and what is in it? --------------------
dim(duluth_df)

names(duluth_df)[13:18]

⚠️ Watch out! NOAA ships its columns as TEMP, YEARMODA, MXSPD — all capitals. clean_names() turns them into temp, yearmoda, mxspd at the moment you load the data, so you never have to think about capital letters again. Leave it off and select(temp) will error.

✏️ Q1: How many rows and columns did you get? What is one row — a day, a month, or a year?

Part 2 · Trim to what you need

▶ Run this:

# Keep only the columns we actually use ----------------
dlh_df <- duluth_df %>%
  select(name, yearmoda, year, month, temp)

head(dlh_df, 3)

47 columns is more than anyone needs to look at. Five answer our question.

✏️ Q2: What is in yearmoda that is not already in year and month? Why do we keep all three?

Part 3 · Plot the raw daily data

▶ Run this:

# Every single daily temperature -----------------------
dlh_df %>%
  ggplot(aes(x = yearmoda, y = temp)) +
  geom_point(alpha = 0.2, size = 0.5)

✏️ Q3a: How tall is the band of points, roughly, in °C? What causes that spread — climate change, or something else?

✏️ Q3b: Can you see a warming trend in this plot? Why not? (Think about the size of the trend compared with the size of the band.)

Part 4 · Summarize to one mean per year

▶ Run this:

# Mean temperature per year ----------------------------
dlh_year_df <- dlh_df %>%
  group_by(year) %>%
  summarize(temp = mean(temp, na.rm = TRUE))

head(dlh_year_df, 3)

Same group_by() + summarize() you have used since Summary Statistics — the grouping variable is just time now.

▶ Run this:

# Yearly means with the regression line ----------------
dlh_year_df %>%
  ggplot(aes(x = year, y = temp)) +
  geom_point() +
  geom_smooth(method = "lm")

✏️ Q4a: How many rows does dlh_year_df have? (Run nrow(dlh_year_df).) How many did you start with?

✏️ Q4b: Now the trend is visible. What did summarizing remove, and what did it keep?

💡 Key idea: Summarizing is not throwing data away — it is choosing which question you are answering. A heat wave question needs daily data. A climate question needs yearly.

Part 5 · Fit the trend and read the rate

▶ Run this:

# Temperature predicted by year ------------------------
dlh_year_model <- lm(temp ~ year, data = dlh_year_df)

summary(dlh_year_model)

Everything you need is in three places:

What you want Where it is
slope — °C per year Coefficients table, year row, Estimate column
p-value for the slope Coefficients table, year row, Pr(>|t|) column
R² Multiple R-squared near the bottom

✏️ Q5a: Write down the slope with units. Is it positive or negative — is Duluth warming or cooling?

✏️ Q5b: Multiply the slope by 10. That is the warming rate in °C per decade. Write it down. Over the 78 years of record, roughly how many degrees is that in total?

✏️ Q5c: Write down the p-value. Is the trend real, or could it be year-to-year luck?

✏️ Q5d: Write down R². Year explains what fraction of the variation in annual temperature? What is the rest of the variation?

⚠️ Watch out! The intercept says about −49 °C. That is the predicted temperature at “year 0,” which is 2,000 years outside your data. The intercept is bookkeeping here — the slope is the science.

🧩 Chunk 3 — Seasons, and your own city (after lecture Chunk 3)

Parts 6–10: label the seasons, fit a model for each, and find your own station.

Part 6 · Label each day with a season

▶ Run this:

# Northern hemisphere seasons --------------------------
dlh_season_df <- dlh_df %>%
  mutate(season = case_when(
    month %in% c(6, 7, 8) ~ "summer",
    month %in% c(12, 1, 2) ~ "winter",
    TRUE ~ "shoulder"
  ))

# Did it work? -----------------------------------------
dlh_season_df %>% count(season)

case_when() reads top to bottom — the first matching line wins. TRUE ~ "shoulder" is the catch-all for everything left over.

✏️ Q6a: How many days landed in each season? Does roughly a quarter / a quarter / a half make sense?

✏️ Q6b: Duluth is at latitude +46.84 — north of the equator. If you ran this exact code on a station in Melbourne, Australia (latitude −37.7), what would be wrong with the labels?

⚠️ Watch out! Seasons are not months — they are hemispheres. June–August is summer north of the equator and winter south of it. A station’s latitude tells you which: positive is northern, negative is southern.

Part 7 · Summarize by year and season

▶ Run this:

# Drop the shoulder, then one mean per year per season -
season_year_df <- dlh_season_df %>%
  filter(season != "shoulder") %>%
  group_by(year, season) %>%
  summarize(temp = mean(temp, na.rm = TRUE),
            .groups = "drop")

head(season_year_df, 4)
  • != means “is not” — this keeps only the two extremes
  • group_by(year, season) uses two grouping columns, so you get one mean per year per season

✏️ Q7: How many rows does season_year_df have? Where does that number come from?

Part 8 · Plot both seasons

▶ Run this:

# One line per season ----------------------------------
season_year_df %>%
  ggplot(aes(x = year, y = temp, color = season)) +
  geom_point() +
  geom_smooth(method = "lm")

Putting color = season inside aes() does two things at once: it colors the points and makes geom_smooth() fit a separate line for each season.

✏️ Q8a: Roughly what temperature is each band centered on?

✏️ Q8b: Which line looks steeper? Why is steepness the thing to compare, and not which band sits higher?

Part 9 · A model for each season

▶ Run this:

# Summer only ------------------------------------------
summer_df <- season_year_df %>%
  filter(season == "summer")

summer_model <- lm(temp ~ year, data = summer_df)
summary(summer_model)
# Winter only ------------------------------------------
winter_df <- season_year_df %>%
  filter(season == "winter")

winter_model <- lm(temp ~ year, data = winter_df)
summary(winter_model)

Filter first, then the same lm() as Part 5. Two small models are much easier to read than one complicated one.

✏️ Q9a: Fill in this table by reading the year row of each summary():

            slope (C/yr)   x10 (C/decade)   p-value
annual
summer
winter

✏️ Q9b: Which season is warming faster, and by roughly how much? Does the annual rate from Part 5 sit between them?

✏️ Q9c: Why does an annual mean hide this? Name one biological thing that depends on winter temperature but not on the annual average.

✏️ Q9d: Write a results sentence for the annual trend. Give the slope with units, the per-decade rate, F with both df, p, and R².

Part 10 · Find the station for your city

▶ Run this:

# Find stations near a place -- Duluth is 46.84, -92.19
nearest_stations(LAT = 46.84, LON = -92.19,
                 distance = 50) %>%
  clean_names() %>%
  select(stnid, name, begin, end) %>%
  head(3)

nearest_stations() takes a latitude, a longitude, and a radius in km.

✏️ Your turn: pick a city — anywhere in the world. Look up its latitude and longitude (any search engine will tell you). Remember:

  • south of the equator → negative latitude
  • west of Greenwich → negative longitude

Then run nearest_stations() with your coordinates.

✏️ Q10a: What city did you pick? What is its latitude and longitude?

✏️ Q10b: Give the stnid, the name, and the begin and end years of the station you chose.

✏️ Q10c: Is your city in the northern or southern hemisphere? Will you need to swap the season labels?

⚠️ Watch out! A trend needs decades. Pick a station whose begin year is 1970 or earlier, and only ask get_GSOD() for years the station actually has — requesting years before begin makes the download fail.

Save your script, then run it top to bottom one last time before you hand it in.

Extension — out of class (~30 min)

Turn this in with your script, as # comments at the bottom of the file. Everyone analyses Duluth — the city is yours.

E1 · Your city’s warming rate

Only one line changes from the in-class code — the station in get_GSOD(). Everything after it is identical:

# The ONLY line that changes ---------------------------
city_df <- get_GSOD(years = 1973:2025,
                    station = "YOUR-STATION-ID") %>%
  clean_names()

# Everything below is the same as Parts 2, 4 and 5 -----
city_year_df <- city_df %>%
  group_by(year) %>%
  summarize(temp = mean(temp, na.rm = TRUE))

city_year_df %>%
  ggplot(aes(x = year, y = temp)) +
  geom_point() +
  geom_smooth(method = "lm")

city_model <- lm(temp ~ year, data = city_year_df)
summary(city_model)

Report, as comments: your city and station ID, the years of record, the slope in °C per decade, the p-value, R², and whether your city is warming faster or slower than Duluth.

E2 · Your city’s seasons

Repeat Parts 6–9 on your city. If your city is in the southern hemisphere, swap the labels in case_when():

# SOUTHERN hemisphere seasons --------------------------
mutate(season = case_when(
  month %in% c(12, 1, 2) ~ "summer",
  month %in% c(6, 7, 8) ~ "winter",
  TRUE ~ "shoulder"
))

Report both seasonal rates in °C per decade. Is the winter-warms-faster pattern you found in Duluth also true in your city?

E3 · Explain your own numbers

  1. Using your R², what fraction of the year-to-year variation in your city’s temperature does year alone explain? What is the rest?
  2. How can a low R² sit next to a p-value below 0.001? What is each one telling you?
  3. One station’s significant trend does not prove why a city is warming. Name one thing that could inflate a single station’s trend with nothing to do with global climate — and say how you might check for it.

End of Activity 12. Next: the final project — the same pipeline, on data you choose.