#' this section is a bit overly complicated as I wanted to run the correlations for all combinations of engineered features for ASI, VHI, NDSI, SWE.
# In the end I only show a few of the resulting plots, but I think the rest may be nice to have at some point - there are tons!!!
df_snow_ready <- df_snow_proc_adj |>
mutate(
mo_label = month(mo, label = T, abbr = T),
snow_type = type
) |>
select(-yr) |>
rename(
yr = yr_adj
)
merge_snow_veg <- function(
df_snow,
snow_param,
df_veg,
veg_param,
start_year = 2000) {
inclusion_range <- df_veg |>
pull(yr) |>
range()
if (is.null(start_year)) {
start_year <- inclusion_range[1]
}
snow_filt <- df_snow |>
filter(
snow_type == snow_param
) |>
filter(yr >= start_year, yr <= inclusion_range[2]) |>
pivot_wider(
id_cols = c(
"adm1_name",
"yr"
), # pivot so 1 row per year
names_from = mo_label,
values_from = c("value", "anom", "anom_z", "value_change"), names_prefix = "snow_"
)
veg_filt <- df_veg |>
filter(
type == veg_param
)
left_join(snow_filt, veg_filt, by = c(
"adm1_name",
"yr"
))
}
corr_all_snow <- function(df_snow = df_snow_ready,
snow_param = "NDSI",
snow_feature = "value",
df_veg = df_fao_yr,
veg_param = "asi",
fao_feature = "mam_mean_value",
start_year = 2000) {
mo_labs <- month(1:12, label = T, abbr = T)
mo_rgx <- glue_collapse(mo_labs, sep = "|")
p_title <- glue("Snow ({snow_param}: {snow_feature}) vs Vegetation ({toupper(veg_param)}: {fao_feature})")
df_merged <- merge_snow_veg(
df_snow = df_snow,
snow_param = snow_param,
df_veg = df_veg,
veg_param = veg_param,
start_year = start_year
) |>
filter(yr < 2025)
df_corrs <- df_merged |>
group_by(
type,
adm1_name,
) |>
summarise(
across(
.cols = any_of(ends_with(c("Nov", "Dec", "Jan", "Feb", "Mar", "Apr"))),
.fns = ~ cor(., !!sym(fao_feature), use = "pairwise.complete.obs")
),
.groups = "drop"
) |>
mutate(
fao_feature = {{ fao_feature }}
)
col_rgx <- glue_collapse(glue("^{snow_feature}_snow_{mo_labs}"), sep = "|")
df_corrs_long <- df_corrs |>
select(
all_of(c(
"type",
"adm1_name"
)), matches(col_rgx)
) |>
pivot_longer(
cols = matches(col_rgx)
) |>
mutate(
mo = factor(str_extract(name, mo_rgx), levels = c("Nov", "Dec", "Jan", "Feb", "Mar", "Apr"))
)
df_corrs_long |>
ggplot(
aes(
x = mo,
y = adm1_name,
fill = value
)
) +
geom_tile() +
scale_fill_gradient2(
low = hdx_hex("tomato-hdx"), # Color for negative values
mid = "white", # Color for zero
high = hdx_hex("mint-hdx"), # Color for positive values
midpoint = 0 # Set midpoint at zero
) +
geom_tile(
data = df_corrs_long |>
filter(
adm1_name == "Faryab",
name %in% c("Nov", "Dec", "Jan", "Feb", "Mar", "Apr")
),
fill = NA, color = "black", lwd = 1.5
) +
geom_text(
aes(label = round(value, 2))
) +
labs(
title = p_title,
# subtitle = "Afghanistan by Province",
y = "Province"
) +
theme(
axis.title.x = element_blank(),
legend.title = element_blank(),
plot.title = element_text(size = 12),
plot.subtitle = element_text(size = 12),
legend.text = element_text(angle = 90)
)
}
all_feature_combos <- expand_grid(
snow_features = c(
"value",
"anom",
"anom_z",
"value_change"
),
fao_features = c(
"mam_mean_value",
"mamj_mean_value",
"may_last_value",
"june_last_value",
"mam_mean_anom",
"mamj_mean_anom",
"mam_mean_z",
"mamj_mean_z"
)
) |>
mutate(
label = glue("Snow: {snow_features} vs Vegetation: {fao_features}")
)
lps_ndsi_asi <- map(
set_names(all_feature_combos$label, all_feature_combos$label),
\(feature_combo){
df_feature <- all_feature_combos |>
filter(label == feature_combo)
corr_all_snow(
df_snow = df_snow_ready,
snow_param = "NDSI",
snow_feature = df_feature$snow_features,
df_veg = df_fao_yr,
veg_param = "asi",
fao_feature = df_feature$fao_features,
start_year = 2000
)
}
)
lps_swe_asi <- map(
set_names(all_feature_combos$label, all_feature_combos$label),
\(feature_combo){
df_feature <- all_feature_combos |>
filter(label == feature_combo)
corr_all_snow(
df_snow = df_snow_ready,
snow_param = "SWE",
snow_feature = df_feature$snow_features,
df_veg = df_fao_yr,
veg_param = "asi",
fao_feature = df_feature$fao_features,
start_year = 2000
)
}
)
