df_all_maps
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lubridate)
library(sf) ## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(tigris) ## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.Separate by season
df_all <- read_csv("~/Downloads/results.csv", show_col_types = FALSE) df_clean <- df_all |> filter(!is.na(obst))
df_clean <- df_clean |>
mutate(
date = as.Date(date),
month = month(date),
season = case_when(
month %in% c(12, 1, 2) ~ "Winter",
month %in% c(6, 7, 8) ~ "Summer",
month %in% c(3, 4, 5) ~ "Spring",
month %in% c(9, 10, 11) ~ "Fall"
)
)nc <- tigris::states(cb = TRUE) |>
filter(STUSPS == "NC")## Retrieving data for the year 2024## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%Tropical nights
Plotted across time
library(ggplot2)
library(viridis)## Loading required package: viridisLitesummer_stats <- df_clean |>
filter(season == "Summer") |>
group_by(year) |>
summarise(
total_tropical_nights = sum(mint > 68, na.rm = TRUE),
n_stations_with_data = n_distinct(uid[!is.na(mint)]),
.groups = "drop"
) |>
mutate(
tropical_nights_per_station = total_tropical_nights / n_stations_with_data
)
ggplot(summer_stats, aes(x = year, y = tropical_nights_per_station)) +
geom_col(fill = "blue", alpha = 0.7) +
geom_smooth(method = "loess", se = FALSE, color = "black", linewidth = 0.5) +
labs(
title = "Number of Summer Tropical Nights per Year",
subtitle = "Data averaged across NC stations",
x = "Year",
y = "Tropical Night Count"
) +
scale_y_continuous(limits = c(0, 92)) ## `geom_smooth()` using formula = 'y ~ x'
Linear model fit to plot above
stations_yearly <- df_clean |>
filter(season == "Summer" & !is.na(mint)) |>
group_by(uid, lat, lon, year) |>
summarise(
n_days = n(),
tropical_night_pct = 100 * sum(mint > 68, na.rm = TRUE) / n_days,
.groups = "drop"
)
stations_yearly <- stations_yearly |>
mutate(year_centered = year - min(year))
tropical_lm <- lm(tropical_night_pct ~ year_centered, data = stations_yearly)
summary(tropical_lm)##
## Call:
## lm(formula = tropical_night_pct ~ year_centered, data = stations_yearly)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.581 -28.644 1.049 23.238 64.759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.32988 0.93087 30.434 < 2e-16 ***
## year_centered 0.31413 0.04068 7.722 1.5e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.26 on 3429 degrees of freedom
## Multiple R-squared: 0.01709, Adjusted R-squared: 0.0168
## F-statistic: 59.62 on 1 and 3429 DF, p-value: 1.498e-14Plotted across state
stations_sf <- df_clean |>
filter(season == "Summer" & !is.na(mint)) |>
filter(year == 2024) |> group_by(uid, lat, lon) |>
summarise(n_days = n(),
tropical_night_pct = 100 * sum(mint > 68, na.rm = TRUE) / n_days,
.groups = "drop" ) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
ggplot() + geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = stations_sf, aes(color = tropical_night_pct), size = 3) +
scale_color_viridis(option = "plasma", name = "Tropical Nights (%)") +
labs( title = "Summer Tropical Nights Across North Carolina in 2024",
subtitle = "Percentage of days with Tmin > 65°F per station" ) +
theme_void()
Plotted across state by decade
# Select the years of interest
years_to_plot <- c(1985, 1995, 2005, 2015)
# Prepare station-level summary for each year
stations_sf <- df_clean |>
filter(season == "Summer" & !is.na(mint) & year %in% years_to_plot) |>
group_by(uid, lat, lon, year) |>
summarise(
n_days = n(),
tropical_night_pct = 100 * sum(mint > 68, na.rm = TRUE) / n_days,
.groups = "drop"
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
# Create panel plot
ggplot() +
geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = stations_sf, aes(color = tropical_night_pct), size = 1) +
scale_color_viridis(option = "plasma", name = "Tropical Nights (%)") +
facet_wrap(~year, ncol = 2) + # creates one panel per year
labs(
title = "Summer Tropical Nights Across North Carolina",
subtitle = "Percentage of days with Tmin > 68°F per station"
) +
theme_void()
Extreme temperatures
plotted across time, not further explored
summer_extreme <- df_clean |>
filter(season == "Summer" & !is.na(maxt)) |>
group_by(year) |>
summarise(
total_extreme_days = sum(maxt > 95, na.rm = TRUE),
n_stations_with_data = n_distinct(uid[!is.na(maxt)]),
.groups = "drop"
) |>
mutate(
extreme_days_per_station = total_extreme_days / n_stations_with_data
)
ggplot(summer_extreme, aes(x = year, y = extreme_days_per_station)) +
geom_col(fill = "red", alpha = 0.7) +
geom_smooth(method = "loess", se = FALSE, color = "black", linewidth = 0.5) +
labs(
title = "Number of Summer Extreme Heat Days per Year",
subtitle = "Data averaged across NC stations",
x = "Year",
y = "Extreme Heat Day Count"
) +
scale_y_continuous(limits = c(0, 92)) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
plotted across state
stations_extreme_sf <- df_clean |>
filter(season == "Summer" & !is.na(maxt)) |>
group_by(uid, lat, lon) |>
summarise(
n_days = n(),
extreme_heat_pct = 100 * sum(maxt > 95, na.rm = TRUE) / n_days,
.groups = "drop"
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
ggplot() +
geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = stations_extreme_sf, aes(color = extreme_heat_pct), size = 3) +
scale_color_viridis(option = "plasma", name = "Extreme Heat Days (%)") +
labs(
title = "Summer Extreme Heat Days Across North Carolina summed from 1985-2024",
subtitle = "Percentage of days with Tmax > 95°F per station"
) +
theme_void() 
Graphing tmin and tmax
tmin/tmax and filter of season changed to explore different temps/seasons
summer_data <- df_clean |>
filter(season == "Summer")
summer_tmin_year <- summer_data |>
group_by(year) |>
summarise(
avg_tmin = mean(mint, na.rm = TRUE),
.groups = "drop"
)
ggplot(summer_tmin_year, aes(x = year, y = avg_tmin)) +
geom_col(fill = "steelblue") +
geom_smooth(method = "loess", se = FALSE, color = "black", linewidth = 0.5) +
labs(
title = "Average Summer Min Temps Across Stations by Year",
x = "Year",
y = "Average Temp (°F)"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
for paper
library(dplyr)
library(ggplot2)
library(tidyr)
# Prepare data: Summer and Winter, Tmin and Tmax
season_data <- df_clean %>%
filter(season %in% c("Summer", "Winter")) %>%
group_by(season, year) %>%
summarise(
avg_tmin = mean(mint, na.rm = TRUE),
avg_tmax = mean(maxt, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(
cols = c(avg_tmin, avg_tmax),
names_to = "temp_type",
values_to = "temperature"
) %>%
mutate(
temp_type = recode(temp_type,
"avg_tmin" = "Min",
"avg_tmax" = "Max")
) %>%
# Compute baseline for 1985 per season and temp type
group_by(season, temp_type) %>%
mutate(
baseline = temperature[year == 1985],
temperature_rel = temperature - baseline
) %>%
ungroup()
# Plot panel
ggplot(season_data, aes(x = year, y = temperature_rel)) +
geom_col(fill = "red", width = 0.7, alpha = 0.7) +
geom_smooth(color = "black", method = "loess", se = FALSE, linewidth = 0.5, alpha = 0.6) +
facet_grid(temp_type ~ season, scales = "fixed") + # 2x2 panel: rows = Min/Max, cols = Summer/Winter
labs(
title = "Change in Seasonal Minimum and Maximum Temps",
x = "Year",
y = "Δ Temperature (°F from 1985)"
) +
theme_minimal() +
theme(
legend.position = "none",
strip.text = element_text(face = "bold"),
panel.grid = element_blank() # removes all grid marks
)## `geom_smooth()` using formula = 'y ~ x'
Looking at coeff from plot above,
more for exploration than anything, not great in terms of methodological approach
library(dplyr)
library(ggplot2)
library(tidyr)
library(broom)
# Prepare data: Summer and Winter, Tmin and Tmax
season_data <- df_clean %>%
filter(season %in% c("Summer", "Winter")) %>%
group_by(season, year) %>%
summarise(
avg_tmin = mean(mint, na.rm = TRUE),
avg_tmax = mean(maxt, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(
cols = c(avg_tmin, avg_tmax),
names_to = "temp_type",
values_to = "temperature"
) %>%
mutate(
temp_type = recode(temp_type,
"avg_tmin" = "Min",
"avg_tmax" = "Max")
) %>%
# Compute baseline for 1985 per season and temp type
group_by(season, temp_type) %>%
mutate(
baseline = temperature[year == 1985],
temperature_rel = temperature - baseline
) %>%
ungroup()
# 1️⃣ Fit linear regression for each panel (season × temp_type)
reg_results <- season_data %>%
group_by(season, temp_type) %>%
do(tidy(lm(temperature_rel ~ year, data = .))) %>%
select(season, temp_type, term, estimate, std.error, statistic, p.value) %>%
rename(coef = estimate, coef_se = std.error, t_value = statistic, p_value = p.value)
print(reg_results)## # A tibble: 8 × 7
## # Groups: season, temp_type [4]
## season temp_type term coef coef_se t_value p_value
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Summer Max (Intercept) 23.1 39.5 0.584 0.563
## 2 Summer Max year -0.0111 0.0197 -0.561 0.578
## 3 Summer Min (Intercept) -91.0 26.3 -3.46 0.00134
## 4 Summer Min year 0.0462 0.0131 3.52 0.00113
## 5 Winter Max (Intercept) -102. 71.9 -1.42 0.164
## 6 Winter Max year 0.0528 0.0359 1.47 0.149
## 7 Winter Min (Intercept) -138. 63.2 -2.18 0.0357
## 8 Winter Min year 0.0707 0.0315 2.24 0.0308tmin and max across state by season, can manipulate season and outcome in code
winter_stations <- df_clean |>
filter(season == "Winter" & !is.na(mint)) |>
group_by(uid, lat, lon) |>
summarise(
avg_winter_tmin = mean(mint, na.rm = TRUE),
n_days = n(),
.groups = "drop"
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
ggplot() +
geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = winter_stations, aes(color = avg_winter_tmin), size = 3) +
scale_color_viridis(option = "plasma", name = "Avg Winter Tmin (°F)") +
labs(
title = "Average Winter Minimum Temperatures Across NC Stations",
subtitle = "Data from 1985-2024"
) +
theme_void() 
Plotted across state by decade, can also change code to explore dif outcome or season
# Select the years of interest
years_to_plot <- c(1985, 1995, 2005, 2015)
# Prepare station-level summary for each year
winter_stations_sf <- df_clean |>
filter(season == "Summer" & !is.na(mint) & year %in% years_to_plot) |>
group_by(uid, lat, lon, year) |>
summarise(
avg_winter_tmin = mean(mint, na.rm = TRUE),
n_days = n(),
.groups = "drop"
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
# Create panel plot
ggplot() +
geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = winter_stations_sf, aes(color = avg_winter_tmin), size = 1) +
scale_color_viridis(option = "plasma", name = "Avg Winter Tmin (°F)") +
facet_wrap(~year, ncol = 2) + # one map per year
labs(
title = "Average Summer Minimum Temperatures"
) +
theme_void()
Linear model, can also manipulate
# Filter for winter data and remove missing Tmin
winter_tmin <- df_clean |>
filter(season == "Winter" & !is.na(mint))
# Aggregate to yearly average Tmin per station or overall
yearly_tmin <- winter_tmin |>
group_by(year) |>
summarise(avg_tmin = mean(mint, na.rm = TRUE)) |>
mutate(year_centered = year - min(year))
# Fit linear model
tmin_lm <- lm(avg_tmin ~ year_centered, data = yearly_tmin)
summary(tmin_lm)##
## Call:
## lm(formula = avg_tmin ~ year_centered, data = yearly_tmin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.363 -1.254 0.276 1.394 5.747
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.33073 0.71450 41.051 <2e-16 ***
## year_centered 0.07072 0.03153 2.243 0.0308 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.302 on 38 degrees of freedom
## Multiple R-squared: 0.1169, Adjusted R-squared: 0.09367
## F-statistic: 5.031 on 1 and 38 DF, p-value: 0.03082Across year
# Select the years of interest
years_to_plot <- c(1985, 1995, 2005, 2015)
# Prepare station-level summary for each year
yearly_stations_sf <- df_clean |>
filter(!is.na(mint) & year %in% years_to_plot) |> # no seasonal filter
group_by(uid, lat, lon, year) |>
summarise(
avg_year_tmin = mean(mint, na.rm = TRUE),
n_days = n(),
.groups = "drop"
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
# Create panel plot
ggplot() +
geom_sf(data = nc, fill = "white", color = "black") +
geom_sf(data = yearly_stations_sf, aes(color = avg_year_tmin), size = 1) +
scale_color_viridis(option = "plasma", name = "Tmin (°F)") +
facet_wrap(~year, ncol = 2) + # one map per year
labs(
title = "Average Yearly Minimum Temperatures Across NC Stations"
) +
theme_void() 
ai generated interactive map, commented out
```{# winter_year_stations_sf <- st_as_sf(} # winter_year_stations, # coords = c(“lon”, “lat”), # crs = 4326 # )
plot_ly(
data = winter_year_stations,
type = ‘scattergeo’,
lat = ~lat,
lon = ~lon,
color = ~avg_winter_tmax,
colors = viridis::plasma(100),
text = ~paste(
“Station:”, uid,
“
Avg Tmin:”, round(avg_winter_tmax, 1), “°F”
),
marker = list(size = 6),
frame = ~year
) %>%
layout(
geo = list(
scope = ‘usa’,
projection = list(type = ‘albers usa’),
lataxis = list(
range = c(min(winter_year_stations\(lat) - 0.5, # max(winter_year_stations\)lat) + 0.5)
),
lonaxis = list(
range = c(min(winter_year_stations\(lon) - 0.5, # max(winter_year_stations\)lon) + 0.5)
),
showland = TRUE,
landcolor = ‘white’,
subunitcolor = ‘black’
),
title = “Average Winter Maximum Temperatures Across NC Stations”
)
```