EDA - Corn VPD and SM
2026-09-18
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggiraph)
library(tigris)
library(sf)
library(readr)
library(stringr)
library(gt)
library(tidyverse)
library(gstat)
library(stars)
library(sp)
library(tseries)options(tigris_use_cache = TRUE)Iowa Data
Source: USDA Corn Yields
Centroids
counties_map <- counties(state = "IA", cb = TRUE, resolution = "5m", year = 2024) %>%
st_transform(2163) %>%
st_centroid() %>%
mutate(
lon = st_coordinates(.)[,1],
lat = st_coordinates(.)[,2]
) %>%
st_drop_geometry() %>%
select(NAME, COUNTYFP, STATEFP, lon, lat)corn1 <- read.csv('usda/07ADC3C5-AAE5-3172-AA91-95CE5361B2E9.csv')
corn2 <- read.csv('usda/02661705-12EA-3FD5-ADB4-F9CEAF2AC953.csv')
corn3 <- read.csv('usda/A266FEE7-3BFC-303A-8816-2B36248CBF0E.csv')
corn <- rbind(corn1, corn2, corn3)corn <- corn %>%
dplyr::filter(Year >= 1984 & State == "IOWA") %>%
dplyr::select(c(Year, Ag.District, County, Value)) %>%
dplyr::mutate(County = stringr::str_to_title(County)) %>%
dplyr::filter(County != "Other Counties") %>%
dplyr::filter(County != "Other (Combined) Counties") %>%
dplyr::group_by(Ag.District) %>%
dplyr::mutate(
Yield_detrended = residuals(lm(Value ~ Year)) + mean(Value, na.rm = TRUE)
) %>%
dplyr::ungroup()Stationarity Tests
tseries::adf.test(corn$Value)##
## Augmented Dickey-Fuller Test
##
## data: corn$Value
## Dickey-Fuller = -7.5813, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationarytseries::kpss.test(corn$Value)##
## KPSS Test for Level Stationarity
##
## data: corn$Value
## KPSS Level = 24.752, Truncation lag parameter = 10, p-value = 0.01In the above analysis, it is clear that the KPSS and ADF tests disagree. This implies that corn yield is “trend stationary” (i.e. it has a deterministic trend rather than a stochastic trend).
Detrending Target
corn_plot <- corn %>%
tidyr::pivot_longer(cols = c(Value, Yield_detrended),
names_to = "Type", values_to = "Yield"
) %>%
dplyr::mutate(Type = ifelse(Type == "Value", "Original", "Detrended"))
p <- ggplot(corn_plot, aes(x = Year, y = Yield, color = Type, fill = Type)) +
geom_point(alpha = 0.55, size = 1.8, stroke = 0) +
geom_smooth(method = "lm", linewidth = 1.1, alpha = 0.18, se = TRUE) +
facet_wrap(~ Ag.District, ncol = 3) +
scale_color_manual(values = c("Original" = "#0072B2", "Detrended" = "#D55E00")) +
scale_fill_manual(values = c("Original" = "#0072B2", "Detrended" = "#D55E00")) +
labs(title = "Corn Yields: Original vs Detrended", x = "Year", y = "Yield", color = "Series", fill = "Series") +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.text.x = element_text(angle = 0, hjust = 1), # Rotates years so they don't overlap
axis.ticks.x = element_line(color = "grey80"),
legend.position = "none"
)
p
ggsave("detrending.png", dpi=150)ggplot(corn, aes(x = factor(Year), y = Yield_detrended)) +
geom_boxplot(fill = "#D55E00", color = "#0072B2", alpha = 0.5, outlier.alpha = 0.3) +
labs(title = "Corn Yield By Year", y = "Yield (Bushel per Acre)", x = "Year") +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.text.x = element_text(angle = 90, hjust = 1), # Rotates years so they don't overlap
axis.ticks.x = element_line(color = "grey80"),
legend.position = "none"
)
corn <- corn %>%
dplyr::mutate(yield_quartile = ntile(Yield_detrended, 4))map_data = corn %>%
dplyr::group_by(County) %>%
dplyr::summarise(Average_Yield = mean(Yield_detrended))# Get Iowa county polygons
iowa_counties <- counties(state = "IA", cb = TRUE, resolution = "5m", year = 2024) %>%
st_transform(2163) %>%
mutate(NAME = str_to_title(NAME))
iowa_map <- iowa_counties %>%
left_join(map_data, by = c("NAME" = "County"))
ggplot(iowa_map) +
geom_sf(aes(fill = Average_Yield), color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "YlOrRd", name = "Avg Yield\n(Bu/Acre)", na.value = "grey95") +
labs(
title = "Average Detrended Corn Yield by County",
subtitle = "Iowa, 1984-present",
caption = "Source: USDA"
) +
theme_void(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 8)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
legend.position = "right"
)
Key Notes:
- For the data below, leap years add an additional day. For all plots below, February 29 is dropped from the data sets.
Soil Moisture
Source: Copernius ERA5
aggregate_to_counties <- function(data, state = "IA") {
day_cols <- names(data)[4:ncol(data)]
iowa_counties <- counties(state = state, cb = TRUE, year = 2021)
iowa_counties_sf <- st_as_sf(iowa_counties)
data_sf <- st_as_sf(data,
coords = c("longitude", "latitude"),
crs = 4269)
data_with_county <- st_join(data_sf, iowa_counties_sf["NAME"])
data_with_county <- st_drop_geometry(data_with_county)
data_with_county <- data_with_county %>%
rename(county = NAME)
aggregated_data <- data_with_county %>%
group_by(county, year) %>%
summarise(across(all_of(day_cols), ~mean(.x, na.rm = TRUE)), .groups = 'drop') %>%
arrange(county, year)
return(aggregated_data)
}swvl1 <- read_csv("copernicus/swvl1_ave_timeseries.csv") %>%
dplyr::filter(year >= 1984) %>%
dplyr::filter(year <= 2024)swvl2 <- read_csv("copernicus/swvl2_ave_timeseries.csv") %>%
dplyr::filter(year >= 1984) %>%
dplyr::filter(year <= 2024)swvl3 <- read_csv("copernicus/swvl3_ave_timeseries.csv") %>%
dplyr::filter(year >= 1984) %>%
dplyr::filter(year <= 2024)pivot_swvl <- function(df, Yield_name) {
df %>%
pivot_longer(
cols = -c(latitude, longitude, year),
names_to = "date",
values_to = Yield_name
) %>%
mutate(date = as.Date(paste(year, date), format = "%Y %d-%b"))
}
long1 <- pivot_swvl(swvl1, "swvl1")
long2 <- pivot_swvl(swvl2, "swvl2")
long3 <- pivot_swvl(swvl3, "swvl3")
smvl_long <- long1 %>%
inner_join(long2, by = c("latitude", "longitude", "year", "date")) %>%
inner_join(long3, by = c("latitude", "longitude", "year", "date")) %>%
mutate(vwc_1m = ((swvl1 * 70) + (swvl2 * 210) + (swvl3 * 720)) / 1000) %>%
select(latitude, longitude, year, date, vwc_1m)
smvl <- smvl_long %>%
mutate(day_label = format(date, "%d-%b")) %>% # "01-Jan", "29-Feb", etc.
select(-date) %>%
pivot_wider(names_from = "day_label", values_from = "vwc_1m")smvl_counties <- aggregate_to_counties(smvl) %>%
dplyr::select(-c("29-Feb", "NA"))corn_quartiles <- corn %>%
dplyr::select(County, Year, yield_quartile) %>%
dplyr::rename(county = County, year = Year)
smvl_counties %>%
left_join(corn_quartiles, by = c("county", "year")) %>%
filter(!is.na(yield_quartile)) %>%
pivot_longer(cols = -c(county, year, yield_quartile),
names_to = "day_label", values_to = "vwc_1m") %>%
mutate(
day_label = as.Date(paste("2001", day_label), format = "%Y %d-%b"),
yield_quartile = factor(yield_quartile, labels = c("Quartile 1 (Low)", "Quartile 2", "Quartile 3", "Quartile 4 (High)"))
) %>%
ggplot(aes(x = day_label, y = vwc_1m, color = yield_quartile, fill = yield_quartile)) +
stat_summary(fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(fun.data = mean_se, geom = "ribbon", alpha = 0.15, color = NA) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
scale_color_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
labs(
title = "Soil Moisture by Day of Year and Yield Quartile",
subtitle = "Iowa counties - mean ± SE across county-years",
x = NULL,
y = "Soil Moisture (m³/m³)"
) +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.ticks.x = element_line(color = "grey80"),
legend.position = "bottom",
legend.title = element_blank()
)
Vapor Pressure Deficit
t2m <- read_csv("copernicus/t2m_ave_timeseries.csv") %>%
dplyr::filter(year >= 1984) %>%
dplyr::filter(year <= 2024)d2m <- read_csv("copernicus/d2m_ave_timeseries.csv") %>%
dplyr::filter(year >= 1984) %>%
dplyr::filter(year <= 2024)svp <- t2m %>%
mutate(
across(-c(latitude, longitude, year), ~ . - 273.15),
across(-c(latitude, longitude, year), ~ 6.112 * exp((17.67 * .) / (. + 243.5)), .names = "svp_{.col}")
) %>%
select(latitude, longitude, year, starts_with("svp_"))avp <- d2m %>%
mutate(
across(-c(latitude, longitude, year), ~ . - 273.15),
across(-c(latitude, longitude, year), ~ 6.112 * exp((17.67 * .) / (. + 243.5)), .names = "avp_{.col}")
) %>%
select(latitude, longitude, year, starts_with("avp_"))vpd <- svp %>%
left_join(avp, by = c("latitude", "longitude", "year")) %>%
mutate(across(
starts_with("svp_"),
~ . - get(str_replace(cur_column(), "svp_", "avp_")),
.names = "{str_replace(.col, 'svp_', '')}"
)) %>%
select(latitude, longitude, year, matches("^\\d{2}-[A-Z][a-z]{2}$"))vpd_counties <- aggregate_to_counties(vpd) %>%
dplyr::select(-c("29-Feb"))vpd_counties %>%
left_join(corn_quartiles, by = c("county", "year")) %>%
filter(!is.na(yield_quartile)) %>%
pivot_longer(cols = -c(county, year, yield_quartile),
names_to = "day_label", values_to = "vwc_1m") %>%
mutate(
day_label = as.Date(paste("2001", day_label), format = "%Y %d-%b"),
yield_quartile = factor(yield_quartile, labels = c("Quartile 1 (Low)", "Quartile 2", "Quartile 3", "Quartile 4 (High)"))
) %>%
ggplot(aes(x = day_label, y = vwc_1m, color = yield_quartile, fill = yield_quartile)) +
stat_summary(fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(fun.data = mean_se, geom = "ribbon", alpha = 0.15, color = NA) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
scale_color_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
labs(
title = "Vapor Pressure Deficit by Day of Year and Yield Quartile",
subtitle = "Iowa counties - mean ± SE across county-years",
x = NULL,
y = "Vapor Pressure Deficit"
) +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.ticks.x = element_line(color = "grey80"),
legend.position = "bottom",
legend.title = element_blank()
)
Great Plains Data
Source: UCI
Centroids
corn_belt_states <- c("IA", "IL", "IN", "KS", "KY", "MI",
"MN", "MO", "NC", "NE", "OH", "PA",
"WI", "SD", "ND", "TN", "VA", "WV",
"OK", "AR")
counties_sf <- counties(state = corn_belt_states, cb = TRUE, resolution = "5m", year = 2024) %>%
st_transform(4269) %>%
select(NAME, NAMELSAD, STATEFP, COUNTYFP)
states_sf <- states(cb = TRUE, resolution = "5m", year = 2024) %>%
st_transform(4269) %>%
filter(STUSPS %in% corn_belt_states)
counties_map_vb <- counties_sf %>%
st_point_on_surface() %>%
mutate(lon = st_coordinates(.)[,1],
lat = st_coordinates(.)[,2]) %>%
st_drop_geometry() %>%
select(NAME, NAMELSAD, STATEFP, COUNTYFP, lon, lat)
add_county_names <- function(pts_df) {
pts_sf <- pts_df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4269, remove = FALSE)
within_join <- st_join(pts_sf, counties_sf, join = st_within)
nearest_join <- st_join(pts_sf, counties_sf, join = st_nearest_feature)
within_join %>%
mutate(
NAME = coalesce(NAME, nearest_join$NAME),
NAMELSAD = coalesce(NAMELSAD, nearest_join$NAMELSAD),
STATEFP = coalesce(STATEFP, nearest_join$STATEFP),
COUNTYFP = coalesce(COUNTYFP, nearest_join$COUNTYFP)
) %>%
st_drop_geometry()
}load('/Users/michael/Library/CloudStorage/OneDrive-UniversityofMissouri/Hatley_research/Veronica_folder/Scalar-on-Function Regression (1 Year)/vpd_weekly.Rdata')
load('/Users/michael/Library/CloudStorage/OneDrive-UniversityofMissouri/Hatley_research/Veronica_folder/Scalar-on-Function Regression (1 Year)/sm_weekly.Rdata')
load('/Users/michael/Library/CloudStorage/OneDrive-UniversityofMissouri/Hatley_research/Veronica_folder/Scalar-on-Function Regression (1 Year)/yieldtouse.Rdata')
load('/Users/michael/Library/CloudStorage/OneDrive-UniversityofMissouri/Hatley_research/Veronica_folder/Scalar-on-Function Regression (1 Year)/locstouse.Rdata')yield <- as.data.frame(use_yield)
markers <- as.data.frame(lonlatstouse)reshape_weekly <- function(arr, value_name = "value") {
n_weeks <- dim(arr)[1]
n_counties <- dim(arr)[2]
n_years <- dim(arr)[3]
purrr::map_dfr(1:n_years, function(y) {
mat <- t(arr[, , y])
df <- as.data.frame(mat)
names(df) <- paste0("W", sprintf("%02d", 1:n_weeks))
df$year_index <- y
df$county_index <- 1:n_counties
df
})
}
sm_vb <- reshape_weekly(sm_weekly)
vpd_vb <- reshape_weekly(vpd_weekly)year_map <- c("1" = 2015, "2" = 2016, "3" = 2017, "4" = 2018)
yield <- as.data.frame(use_yield) %>%
mutate(county_index = row_number()) %>%
rename(Y1 = V1, Y2 = V2, Y3 = V3, Y4 = V4) %>%
pivot_longer(-county_index, names_to = "year_index", values_to = "yield") %>%
mutate(
year_index = as.integer(str_remove(year_index, "Y")),
year = year_map[as.character(year_index)],
yield_quartile = ntile(yield, 4)
)
markers <- as.data.frame(lonlatstouse) %>%
mutate(county_index = row_number())
sm_vb <- sm_vb %>%
left_join(markers, by = "county_index") %>%
left_join(yield, by = c("county_index", "year_index")) %>%
select(-year_index)
vpd_vb <- vpd_vb %>%
left_join(markers, by = "county_index") %>%
left_join(yield, by = c("county_index", "year_index")) %>%
select(-year_index)Yield Plots
map_data <- yield %>%
left_join(markers, by = "county_index") %>%
select(-c(county_index, year_index)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4269) %>%
st_join(counties_sf, join = st_nearest_feature) %>%
st_drop_geometry()
map_summary <- map_data %>%
group_by(NAME, STATEFP, COUNTYFP) %>%
summarise(avg_yield = mean(yield, na.rm = TRUE), .groups = "drop")
corn_belt_map <- counties_sf %>%
left_join(map_summary, by = c("NAME", "STATEFP", "COUNTYFP"))ggplot(corn_belt_map) +
geom_sf(aes(fill = avg_yield), color = "white", linewidth = 0.1) +
geom_sf(data = states_sf, fill = NA, color = "grey40", linewidth = 0.4) +
scale_fill_viridis_c(option = "YlOrRd", name = "Avg Yield\n(Bu/Acre)", na.value = "grey95") +
labs(title = "Average Corn Yield by County", subtitle = "Great Plains") +
theme_void(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 8)),
plot.margin = margin(16, 16, 16, 16),
legend.position = "right"
)
fips_lookup <- fips_codes %>%
distinct(state_code, state)
yield_plot <- yield %>%
dplyr::left_join(markers, by = "county_index") %>%
add_county_names() %>%
dplyr::left_join(fips_lookup, by = c("STATEFP" = "state_code")) %>%
dplyr::select(-c(county_index, year_index, yield_quartile, NAMELSAD, COUNTYFP, STATEFP)) %>%
dplyr::group_by(NAME, state) %>%
dplyr::mutate(
yield_detrended = residuals(lm(yield ~ year)) + mean(yield, na.rm = TRUE)
) %>%
dplyr::ungroup() %>%
relocate(year, yield, .after = state) %>%
tidyr::pivot_longer(cols = c(yield, yield_detrended),
names_to = "type", values_to = "yield") %>%
dplyr::mutate(type = ifelse(type == "yield", "Original", "Detrended"))
ggplot(yield_plot, aes(x = year, y = yield, color = type, fill = type)) +
geom_point(alpha = 0.55, size = 1.8, stroke = 0) +
geom_smooth(method = "lm", linewidth = 1.1, alpha = 0.18, se = TRUE) +
facet_wrap(~ state, ncol = 3) +
scale_color_manual(values = c("Original" = "#0072B2", "Detrended" = "#D55E00")) +
scale_fill_manual(values = c("Original" = "#0072B2", "Detrended" = "#D55E00")) +
labs(title = "Corn Yields: Original vs Detrended", x = "Year", y = "Yield", color = "Series", fill = "Series") +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey60", size = 9, hjust = 0),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.text.x = element_text(angle = 0, hjust = 1), # Rotates years so they don't overlap
axis.ticks.x = element_line(color = "grey80"),
legend.position = "none"
)
Vapor Pressure Deficit and Soil Moisture Analysis
sm_vb <- add_county_names(sm_vb)
vpd_vb <- add_county_names(vpd_vb)sm_vb <- sm_vb %>%
relocate((ncol(sm_vb)-6):ncol(sm_vb)) %>%
select(-c(COUNTYFP, NAMELSAD, county_index, yield, lon, lat)) %>%
relocate(year, .after = STATEFP)
vpd_vb <- vpd_vb %>%
relocate((ncol(vpd_vb)-6):ncol(vpd_vb)) %>%
select(-c(COUNTYFP, NAMELSAD, county_index, yield, lon, lat)) %>%
relocate(year, .after = STATEFP)vpd_vb %>%
pivot_longer(cols = starts_with("W"),
names_to = "week", values_to = "vpd") %>%
mutate(
week_num = as.integer(str_remove(week, "W")),
yield_quartile = factor(yield_quartile,
labels = c("Quartile 1 (Low)", "Quartile 2",
"Quartile 3", "Quartile 4 (High)"))
) %>%
ggplot(aes(x = week_num, y = vpd, color = yield_quartile, fill = yield_quartile)) +
stat_summary(fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(fun.data = mean_se, geom = "ribbon", alpha = 0.15, color = NA) +
scale_x_continuous(breaks = seq(1, 52, by = 4)) +
scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
scale_color_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
labs(
title = "Vapor Pressure Deficit by Week and Yield Quartile",
subtitle = "Great Plains - mean ± SE across county-years",
x = "Week",
y = "Vapor Pressure Deficit"
) +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.ticks.x = element_line(color = "grey80"),
legend.position = "bottom",
legend.title = element_blank()
)
sm_vb %>%
pivot_longer(cols = starts_with("W"),
names_to = "week", values_to = "vpd") %>%
mutate(
week_num = as.integer(str_remove(week, "W")),
yield_quartile = factor(yield_quartile,
labels = c("Quartile 1 (Low)", "Quartile 2",
"Quartile 3", "Quartile 4 (High)"))
) %>%
ggplot(aes(x = week_num, y = vpd, color = yield_quartile, fill = yield_quartile)) +
stat_summary(fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(fun.data = mean_se, geom = "ribbon", alpha = 0.15, color = NA) +
scale_x_continuous(breaks = seq(1, 52, by = 4)) +
scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
scale_color_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
labs(
title = "Soil Moisture by Week and Yield Quartile",
subtitle = "Great Plains - mean ± SE across county-years",
x = "Week",
y = "Soil Moisture"
) +
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 12)),
plot.margin = margin(16, 16, 16, 16),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(color = "grey30", size = 10),
axis.text = element_text(color = "grey30"),
axis.ticks.x = element_line(color = "grey80"),
legend.position = "bottom",
legend.title = element_blank()
)
Missing Values
missing_vpd <- vpd_vb |>
dplyr::mutate(across(where(is.character), ~ dplyr::na_if(.x, ""))) |>
dplyr::summarise(across(everything(), ~ round(mean(is.na(.x)) * 100, 2))) |>
tidyr::pivot_longer(cols = everything(), names_to = "Variable", values_to = "Missing_pct_vpd")
missing_sm <- sm_vb |>
dplyr::mutate(across(where(is.character), ~ dplyr::na_if(.x, ""))) |>
dplyr::summarise(across(everything(), ~ round(mean(is.na(.x)) * 100, 2))) |>
tidyr::pivot_longer(cols = everything(), names_to = "Variable", values_to = "Missing_pct_sm")
missing <- dplyr::left_join(missing_vpd, missing_sm, by = "Variable")
missing |> gt()| Variable | Missing_pct_vpd | Missing_pct_sm |
|---|---|---|
| yield_quartile | 0 | 0.00 |
| NAME | 0 | 0.00 |
| STATEFP | 0 | 0.00 |
| year | 0 | 0.00 |
| W01 | 0 | 57.76 |
| W02 | 0 | 49.41 |
| W03 | 0 | 52.70 |
| W04 | 0 | 43.97 |
| W05 | 0 | 41.50 |
| W06 | 0 | 53.02 |
| W07 | 0 | 48.86 |
| W08 | 0 | 33.97 |
| W09 | 0 | 31.62 |
| W10 | 0 | 30.37 |
| W11 | 0 | 27.43 |
| W12 | 0 | 26.57 |
| W13 | 0 | 16.42 |
| W14 | 0 | 2.90 |
| W15 | 0 | 2.08 |
| W16 | 0 | 0.90 |
| W17 | 0 | 0.00 |
| W18 | 0 | 0.00 |
| W19 | 0 | 0.00 |
| W20 | 0 | 0.00 |
| W21 | 0 | 0.00 |
| W22 | 0 | 0.00 |
| W23 | 0 | 0.00 |
| W24 | 0 | 0.00 |
| W25 | 0 | 0.00 |
| W26 | 0 | 0.16 |
| W27 | 0 | 0.16 |
| W28 | 0 | 0.16 |
| W29 | 0 | 0.16 |
| W30 | 0 | 0.16 |
| W31 | 0 | 0.16 |
| W32 | 0 | 0.16 |
| W33 | 0 | 0.16 |
| W34 | 0 | 0.16 |
| W35 | 0 | 0.16 |
| W36 | 0 | 0.16 |
| W37 | 0 | 0.00 |
| W38 | 0 | 0.00 |
| W39 | 0 | 0.00 |
| W40 | 0 | 0.00 |
| W41 | 0 | 0.00 |
| W42 | 0 | 0.00 |
| W43 | 0 | 0.00 |
| W44 | 0 | 0.00 |
| W45 | 0 | 0.31 |
| W46 | 0 | 0.74 |
| W47 | 0 | 4.43 |
| W48 | 0 | 6.70 |
| W49 | 0 | 9.99 |
| W50 | 0 | 15.36 |
| W51 | 0 | 14.89 |
| W52 | 0 | 30.49 |
sm_vb |>
summarise(across(starts_with("W"), ~ mean(is.na(.x)) * 100)) |>
pivot_longer(everything(), names_to = "week", values_to = "pct_missing") |>
mutate(week_num = as.integer(str_remove(week, "W"))) |>
ggplot(aes(x = week_num, y = pct_missing)) +
geom_col(fill = "#56B4E9") +
scale_x_continuous(breaks = seq(1, 52, by = 4)) +
labs(title = "Soil Moisture - % Missing by Week",
x = "Week", y = "% Missing") +
theme_minimal()
sm_vb |>
pivot_longer(cols = starts_with("W"), names_to = "week", values_to = "sm") |>
mutate(week_num = as.integer(str_remove(week, "W"))) |>
filter(week_num >= 18, week_num <= 44, is.na(sm)) |>
count(NAME, STATEFP, year, sort = TRUE) |>
gt()| NAME | STATEFP | year | n |
|---|---|---|---|
| Franklin | 47 | 2015 | 11 |
| Franklin | 47 | 2016 | 11 |
| Franklin | 47 | 2017 | 11 |
| Franklin | 47 | 2018 | 11 |
sm_vb |>
pivot_longer(cols = starts_with("W"), names_to = "week", values_to = "sm") |>
mutate(week_num = as.integer(str_remove(week, "W"))) |>
filter(week_num >= 14, week_num <= 45) |>
group_by(NAME, STATEFP) |>
summarise(pct_missing = mean(is.na(sm)) * 100, .groups = "drop") |>
filter(pct_missing > 0) |>
right_join(counties_sf, by = c("NAME", "STATEFP")) |>
st_as_sf() |>
ggplot() +
geom_sf(aes(fill = pct_missing), color = "white", linewidth = 0.1) +
geom_sf(data = states_sf, fill = NA, color = "grey40", linewidth = 0.4) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "% Missing\n(Wk 14-45)", na.value = "grey95") +
labs(title = "Soil Moisture - Mid-Season Missingness",
subtitle = "Weeks 14-45 across all county-years") +
theme_void(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 8)),
plot.margin = margin(16, 16, 16, 16),
legend.position = "right"
)
sm_vb |>
pivot_longer(cols = starts_with("W"), names_to = "week", values_to = "sm") |>
mutate(week_num = as.integer(str_remove(week, "W"))) |>
filter(week_num >= 14, week_num <= 45) |>
group_by(NAME, STATEFP, year) |>
summarise(pct_missing = mean(is.na(sm)) * 100, .groups = "drop") |>
filter(pct_missing > 0) |>
right_join(counties_sf, by = c("NAME", "STATEFP")) |>
filter(!is.na(year)) |>
st_as_sf() |>
ggplot() +
geom_sf(aes(fill = pct_missing), color = "white", linewidth = 0.1) +
geom_sf(data = states_sf, fill = NA, color = "grey40", linewidth=0.4) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "% Missing", na.value = "grey95") +
facet_wrap(~ year, ncol = 2) +
labs(title = "Soil Moisture - Mid-Season Missingness by Year",
subtitle = "Weeks 14-45") +
theme_void(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey40", size = 11, margin = margin(b = 8)),
plot.margin = margin(16, 16, 16, 16),
legend.position = "right",
strip.text = element_text(face = "bold", size = 11)
)