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

November 10, 2026

In-class Activity 12: Real Climate Data

Today’s Objectives

Download real NOAA weather station data from inside R with GSODR
Clean the column names with janitor::clean_names()
Summarize ~25,000 daily records into one mean per year
Fit lm(temp ~ year) and read the warming rate off summary()
Split the year into seasons with case_when() and compare warming rates
Find the station for a city you choose and run it again

📂 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 packages: install these two once, from the Console: install.packages("GSODR") and install.packages("nasapower")
- loading GSODR prints a notice that the GSOD data retired in Aug 2025. That is expected — not an error. The historical archive still downloads fine.

What 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 Homework 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 station data - Duluth, in class
library(nasapower)   # NASA data by lat/long - another source, example only

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

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

You asked for 78 years but you will get 70. NOAA has no files at all for 1965–1972 — a hole in the GSOD archive itself, not something about Duluth; most long-running US stations are missing exactly those eight years. get_GSOD() skips what isn’t there and carries on.

That is worth seeing rather than hiding. The gap shows up as a hole in your plots, real archives have holes, and the regression doesn’t care — it just has 70 points instead of 78.

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

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.

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)

✏️ Q1: Can you see a warming trend in this cloud of daily points? Why not?

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

✏️ Q2: The trend is suddenly 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

✏️ Q3: Is Duluth warming or cooling? Give the rate in °C per decade (the slope × 10).

✏️ Q4: R² is small but p is tiny. What is each one telling you?

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

✏️ Q5: Run this same 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

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.

✏️ Q6: Which line is 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.

✏️ Q7: Which season is warming faster, and by roughly how much?

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

✏️ Q9: Write the results sentence for the annual trend — slope with units, per-decade rate, p, and R².

Part 10 · Pick a place of your own

Duluth is done. For homework you run the same analysis somewhere else — with the same GSODR package you just used, so you get a long station record to compare against Duluth’s.

What GSOD actually has — any country

GSODR ships a table of every station in the world and the years each one covers. This is how you find out what exists before you ask for it.

▶ Run this, with Australia as the example:

# Every GSOD station and the years it covers ------------
load(system.file("extdata", "isd_history.rda", package = "GSODR"))

stations_df <- isd_history %>%
  clean_names() %>%
  filter(country_name == "AUSTRALIA") %>%
  mutate(begin_year = as.integer(str_sub(begin, 1, 4)),
         end_year   = as.integer(str_sub(end, 1, 4))) %>%
  filter(begin_year <= 1970,
         end_year >= 2025) %>%
  select(stnid, name, begin_year, end_year) %>%
  arrange(begin_year)

stations_df

Put your country in filter(country_name == ), in CAPITALS. Australia has 151 stations starting in 1970 or earlier and still reporting; Kenya has 19, India 138, Brazil 49.

Scan the name column for a city you want, and note two things: its stnid and its begin_year.

Download it

stnid goes into get_GSOD(), and begin_year tells you what to put at the start of the year range:

▶ Run this, with Melbourne filled in as an example:

# Your city's whole station record ----------------------
city_df <- get_GSOD(years = 1966:2025,
                    station = "959360-99999") %>%
  clean_names()

dim(city_df)

Replace both the years and the stnid with your own before you go any further.

💡 Key idea: get_GSOD() hands back the same columns as Duluth — year, month, temp, already lower-case from clean_names(). Nothing to rename. Every line of the homework is code you have already run, with city_df in place of dlh_df.

⚠️ Watch out! Stations outside the US are gappier than Duluth. If get_GSOD() returns far fewer years than you asked for, or errors, go back to stations_df and pick another station in the same country — a capital city or a major airport is usually the safest bet.

Another source worth knowing — NASA POWER

You are not using this for the homework, but it is worth seeing that a station is not the only way to get climate data.

nasapower needs no station at all — just a longitude and a latitude — and gives you 1981 to now for anywhere on Earth with no gaps.

▶ Run this once, with London filled in, just to see it work:

# NASA POWER: lonlat = c(LONGITUDE, LATITUDE) -----------
power_df <- get_power(community = "ag",
                      lonlat = c(-0.13, 51.51),
                      pars = "T2M",
                      dates = c("1981-01-01", "2024-12-31"),
                      temporal_api = "daily") %>%
  clean_names() %>%
  rename(month = mm, temp = t2m)

dim(power_df)

Watch the order: longitude first, then latitude. The columns come back as MM and T2M, so rename() is needed here to make them match the names you used all class — that is the kind of small translation every new data source needs.

The two sources side by side

	GSOD station	NASA POWER
find it by	country + `stnid`	longitude, latitude
record	1930s–2025, with gaps	1981–now, no gaps
what it is	a real weather station	a modelled grid cell (~50 km)
columns	`year`, `month`, `temp` already	need `rename()`

Use GSOD for the homework. The longer record is the point — more years means a trend you can actually detect, and it lets you compare your city against the Duluth data you already have in dlh_year_df.

✏️ Your turn: find a station and download it.

✏️ Q10: Which city and stnid did you pick? Northern or southern hemisphere — do you need to swap the season labels?

Save your script. The rest is homework.

Homework — out of class

Add this to the bottom of the same script. Answers as # comments.

H1 · Your city’s warming rate

Parts 4 and 5 again, on city_df:

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)

✏️ H1: Slope in °C per year and per decade, p, and R². Is your city warming?

H2 · Your city’s seasons

Parts 6–9 again. The labels below are northern hemisphere — if your city is south of the equator, swap summer and winter.

city_season_df <- city_df %>%
  mutate(season = case_when(
    month %in% c(6, 7, 8) ~ "summer",
    month %in% c(12, 1, 2) ~ "winter",
    TRUE ~ "shoulder"
  )) %>%
  filter(season != "shoulder") %>%
  group_by(year, season) %>%
  summarize(temp = mean(temp, na.rm = TRUE),
            .groups = "drop")

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

Then fit a model for each season, exactly as in Part 9.

✏️ H2: Both seasonal rates in °C per decade. Does winter warm faster in your city too?

H3 · Is your slope different from Duluth’s?

Your city and Duluth almost certainly gave you different numbers. That does not mean the two rates are really different — the difference might be within the scatter. This tests it.

Both records came from the same source, GSODR, so you already have everything you need: dlh_year_df from Part 4 and city_year_df from H1.

Step 1 · Match the years

The two stations do not start in the same year. Compare them over the span they share, or a difference in place gets mixed up with a difference in years.

# The later of the two start years ----------------------
start_year <- max(min(dlh_year_df$year), min(city_year_df$year))

start_year

Step 2 · Label each one and stack them

# A column saying which place each row came from --------
duluth_labeled_df <- dlh_year_df %>%
  filter(year >= start_year) %>%
  mutate(station = "duluth")

city_labeled_df <- city_year_df %>%
  filter(year >= start_year) %>%
  mutate(station = "my_city")

# Stack them into one data frame ------------------------
both_df <- bind_rows(duluth_labeled_df, city_labeled_df)

# Did it work? -----------------------------------------
both_df %>% count(station)

bind_rows() puts one data frame underneath the other. It matches up columns by name, which is why both frames having year, temp, and station matters — and why using GSODR for both places saved you a rename().

The station column is the whole trick: stacking alone would throw away which rows came from where, so you add the label before you stack.

Step 3 · Plot both places

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

Step 4 · Test the two slopes

# A separate slope for each station --------------------
slope_model <- lm(temp ~ year * station, data = both_df)

summary(slope_model)

temp ~ year * station fits a separate slope for each station. In the Coefficients table:

Row	What it is
`year`	Duluth’s slope
`year:stationmy_city`	how much your city’s slope differs

So your city’s slope is year + year:stationmy_city. The p-value on that second row is the answer: it is a t-test on the difference between two slopes, reported in the same t value and Pr(>|t|) columns you have been reading all class.

p < 0.05 → the two places warm at measurably different rates
p > 0.05 → you cannot tell the two rates apart, even if the numbers look different

⚠️ If that row is exactly 0 with p = 1, you are comparing Duluth with Duluth — go back to Part 10 and put in your own stnid.

✏️ H3a: Both slopes in °C per decade, and the difference between them.

✏️ H3b: One sentence reporting the comparison — both rates, the difference, and the p-value on the year:stationmy_city row. Can you tell the two rates apart?

H4 · Explain your numbers

✏️ H4a: How can a low R² sit next to a p-value below 0.001? What is each one telling you?

✏️ H4b: Why does a longer record make a trend easier to detect?

✏️ H4c: A significant trend at one place does not prove why it is warming. Name one thing that could inflate a single location’s trend with nothing to do with global climate — and how you would check for it.

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

Other Formats

Activity 12 — Real Climate Data in R

In-class Activity 12: Real Climate Data

Today’s Objectives

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

Part 1 · Load libraries and download the data

Part 2 · Trim to what you need

Part 3 · Plot the raw daily data

Part 4 · Summarize to one mean per year

Part 5 · Fit the trend and read the rate

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

Part 6 · Label each day with a season

Part 7 · Summarize by year and season

Part 8 · Plot both seasons

Part 9 · A model for each season

Part 10 · Pick a place of your own

What GSOD actually has — any country

Download it

Another source worth knowing — NASA POWER

Homework — out of class

H1 · Your city’s warming rate

H2 · Your city’s seasons

H3 · Is your slope different from Duluth’s?

Step 1 · Match the years

Step 2 · Label each one and stack them

Step 3 · Plot both places

Step 4 · Test the two slopes

H4 · Explain your numbers