library(sf)
library(terra)
library(tidyverse)
library(ggplot2)
library(tidyterra)
library(ggspatial)
library(ggnewscale)
library(viridis)
library(patchwork)
library(cowplot)
library(RColorBrewer)
library(landscapemetrics)
library(furrr)
library(progressr)
library(future)
library(pastclim)
library(caret)
library(MuMIn)
library(corrplot)
library(car)
library(ggeffects)
library(AICcmodavg)
library(rstatix)
library(dunn.test)
library(broom)
library(spdep)
library(dplyr)
library(tidyr)
library(patchwork) This document describes the workflow used to model german wide acceptance of pollinator friendly measures, and how it is used to predict acceptance changes with climate change. The model is built on landscape metrics calculated from a german wide map of polliantor friendly meassures, CHELSA climate data and answers from a survey conducted with farmers.
Landcover, agricultural and biotope data were combined to create a german wide map looking at pollinator friendly meaussures. The data was reclassified to 8 classes relevant for pollinators. The classes are:
For a detailed explanation look here: https://rpubs.com/jonasagri4pol/1440941.
Using german postal codes as a boundary, several landscape metrics were calculated from the pollinator friendly map. All of them were specified to ‘landscape level’, to look at landscape structures between different patches.
The following 11 landscape-level metrics are calculated per postal code (PLZ) for each of the four raster layers:
| Metric | Description |
|---|---|
np |
Number of patches: overall fragmentation |
pd |
Patch density: NP normalised to 100 ha |
lpi |
Largest patch index: percentage area of largest patch |
mesh |
Effective mesh size: landscape connectivity |
split |
Splitting index: landscape split |
ed |
Edge density: Length of patch boundaries |
shdi |
Shannon Diversity Index |
pr |
Patch richness: number of unique classes |
enn_mn |
Mean nearest-neighbour distance |
cohesion |
Patch cohesion: physical connectivity |
ai |
Aggregation index |
Climate data is extracted from the CHELSA 2.1 dataset (0.5 arc-min resolution) for Germany. 13 bioclimatic varibles relevant to agricultutal practice are downloaded. The data is in a 30year timeframe, present data used reaches from 1981-2010, future data from 2041 - 2070.
| Variable | Description |
|---|---|
bio01 |
Mean Annual Near-Surface Air Temperature |
bio02 |
Mean Diurnal Near-Surface Air Temperature Range |
bio04 |
Temperature Seasonality |
bio05 |
Mean Daily Maximum Near-Surface Air Temperature of the Warmest Month |
bio06 |
Mean Daily Minimum Near-Surface Air Temperature of the Coldest Month |
bio10 |
Mean Daily Mean Near-Surface Air Temperature of the Warmest Quarter |
bio11 |
Mean Daily Mean Near-Surface Air Temperature of the Quarter |
bio12 |
Annual Precipitation |
bio13 |
Precipitation of the Wettest Month |
bio14 |
Precipitation of the Driest Month |
bio15 |
Precipitation Seasonality |
bio16 |
Mean Monthly Precipitation of the Wettest Quarter |
bio18 |
Mean Monthly Precipitation of the Warmest Quarter |
Specified landscape metrics are defined.
lsm_metrics <- c(
"lsm_l_np", "lsm_l_pd", "lsm_l_lpi", "lsm_l_mesh",
"lsm_l_split", "lsm_l_ed", "lsm_l_shdi", "lsm_l_pr",
"lsm_l_enn_mn", "lsm_l_cohesion", "lsm_l_ai"
)Metrics are computed sequentially with a for loop that
checkpoints each postal code as an .rds file.
# PLZ shapefile
plz_de <- st_read("/Users/jonasschreiber/Downloads/PLZ_Gebiete_9143106783908117499-2.gpkg") %>%
st_transform(3035)
poll_friendly_path <- "/Volumes/SSK SSD/agri4pol_jonas/poll_classes/combined_final_10m.tif"
for (i in 1:nrow(plz_de)) {
out_path <- paste0("/Volumes/SSK SSD/agri4pol_jonas/poll_classes/metrics/", i, ".rds")
if (file.exists(out_path)) next
tryCatch({
r_worker <- rast(poll_friendly_path)
plz_single <- vect(plz_de[i, ])
poll_plz <- crop(r_worker, plz_single) %>% mask(plz_single)
metrics <- calculate_lsm(poll_plz, what = lsm_metrics)
metrics$plz <- plz_de$plz[i]
saveRDS(metrics, out_path)
}, error = function(e) saveRDS(NULL, out_path))
if (i %% 100 == 0)
cat("\r", i, "/ 8170 (", round(i / 8170 * 100), "%) ")
}# Collect all .rds files into a single long-format data frame
poll_metrics <- list.files(
"/Volumes/SSK SSD/agri4pol_jonas/poll_classes/metrics/",
full.names = TRUE
) %>%
purrr::map(readRDS) %>%
bind_rows()
write_csv(poll_metrics, "/Volumes/SSK SSD/agri4pol_jonas/poll_classes/lsm_de.csv")
#for wide format:
poll_metrics_wide <- pivot_wider(poll_metrics,
names_from = "metric",
values_from = "value") %>%
dplyr::select(-class, -id, -level, -layer) %>%
rename(
patch_density = pd,
largest_patch_index = lpi,
edge_density = ed,
aggregation_index = ai,
number_patches = np,
patch_richness = pr
)
write_csv(poll_metrics_wide,
"/Volumes/SSK SSD/agri4pol_jonas/poll_classes/lsm_de_wide.csv")Note – Thüringen update: Thüringen biotope data was obtained later and reclassified separately. The PLZ-level metrics were recalculated on the Thüringen-specific raster and merged back into the Germany-wide table.
Metrics from the three individual source layers (CLC, Thünen, biotopes) and the combined layer are compared to check for systematic differences. Variation between maps shows that including the three data sources is needed for the highest spatial resolution.
clc_metrics <- read_csv("/Users/jonasschreiber/Documents/AGRI4POL/geodaten/landschaftsstruktur/landscape_metrics/long/clc_lm_long.csv") %>%
mutate(source = "Landcover (Copernicus)")
thuenen_metrics <- read_csv("/Volumes/SSK SSD/agri4pol_jonas/landscape/crops/crops_metrics.csv") %>%
mutate(source = "Agricultural Landuse (Thünen)")
bt_metrics <- read_csv("/Volumes/SSK SSD/agri4pol_jonas/landscape/gesch_biotope/bt_metrics.csv") %>%
mutate(source = "Protected Biotopes")
alle_metrics <- read_csv("/Volumes/SSK SSD/agri4pol_jonas/poll_classes/lsm_de.csv") %>%
mutate(source = "Combined")
all_metrics <- bind_rows(clc_metrics, thuenen_metrics, bt_metrics, alle_metrics)# Trimmed to 1.5×IQR per group for visualisation
all_metrics_trimmed <- all_metrics %>%
group_by(metric, source) %>%
mutate(q1 = quantile(value, 0.25, na.rm = TRUE),
q3 = quantile(value, 0.75, na.rm = TRUE),
iqr = q3 - q1) %>%
filter(value >= q1 - 1.5 * iqr, value <= q3 + 1.5 * iqr) %>%
dplyr::select(-q1, -q3, -iqr) %>%
ungroup() %>%
mutate(source = factor(source, levels = c(
"Landcover (Copernicus)", "Agricultural Landuse (Thünen)",
"Protected Biotopes", "Combined"
)))
metric_labels <- c(
ai = "Aggregation Index (AI)", cohesion = "Patch Cohesion Index",
ed = "Edge Density (ED)", enn_mn = "Mean ENN Distance",
lpi = "Largest Patch Index (LPI)", mesh = "Effective Mesh Size",
np = "Number of Patches (NP)", pd = "Patch Density (PD)",
pr = "Patch Richness (PR)", shdi = "Shannon Diversity Index",
split = "Splitting Index"
)
#boxplot
ggplot(all_metrics_trimmed, aes(x = source, y = value, fill = source)) +
geom_boxplot(outlier.shape = NA, width = 0.6, staplewidth = 0.5) +
facet_wrap(~ metric, scales = "free_y", labeller = as_labeller(metric_labels)) +
scale_fill_manual(values = c(
"Landcover (Copernicus)" = "#E69F00",
"Agricultural Landuse (Thünen)" = "#56B4E9",
"Protected Biotopes" = "#009E73",
"Combined" = "#CC79A7"
)) +
theme_minimal() +
theme(legend.position = "none", strip.text = element_text(face = "bold", size = 8),
axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank()) +
labs(title = "Landscape Metrics per PLZ",
subtitle = "CLC, Thünen, Protected Biotopes and Combined",
y = "Value",
caption = "Values outside 1.5×IQR per group excluded for visualisation.")Climatic variables are the clipped to postal code boundaries and a mean is calculated per postal code. After that they are joined to landscape metrics by postal code.
set_data_path("/Volumes/SSK SSD/agri4pol_jonas/temp")
# Download with pastclim
download_dataset(
dataset = "CHELSA_2.1_0.5m",
bio_variables = c("bio01", "bio02", "bio04", "bio05", "bio06",
"bio10", "bio11", "bio12", "bio13", "bio14",
"bio15", "bio16", "bio18")
)#climate per postal code
metrics_climate_de <- poll_metrics_wide %>%
left_join(
current_climate_de %>%
dplyr::select(plz, bio01, bio02, bio04, bio05, bio06,
bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio18),
by = "plz"
)
current_files <- c(
bio01 = "CHELSA_bio1_1981-2010_V.2.1.tif",
bio02 = "CHELSA_bio2_1981-2010_V.2.1.tif",
bio04 = "CHELSA_bio4_1981-2010_V.2.1.tif",
bio05 = "CHELSA_bio5_1981-2010_V.2.1.tif",
bio06 = "CHELSA_bio6_1981-2010_V.2.1.tif",
bio10 = "CHELSA_bio10_1981-2010_V.2.1.tif",
bio11 = "CHELSA_bio11_1981-2010_V.2.1.tif",
bio12 = "CHELSA_bio12_1981-2010_V.2.1.tif",
bio13 = "CHELSA_bio13_1981-2010_V.2.1.tif",
bio14 = "CHELSA_bio14_1981-2010_V.2.1.tif",
bio15 = "CHELSA_bio15_1981-2010_V.2.1.tif",
bio16 = "CHELSA_bio16_1981-2010_V.2.1.tif",
bio18 = "CHELSA_bio18_1981-2010_V.2.1.tif"
)
current_climate_de <- data.frame(plz = plz_de$plz)
for (var in names(current_files)) {
r <- rast(file.path(chelsa_dir, current_files[var]))
r_de <- crop(r, plz_vect)
zonal_mean <- terra::extract(r_de, plz_vect, fun = mean, na.rm = TRUE, ID = FALSE)
current_climate_de[[var]] <- zonal_mean[, 1]
cat(var, "done\n")
}
write.csv(current_climate_de,
file.path(output_dir, "current_climate_de_zonal.csv"),
row.names = FALSE)
#join with landscape metrics
metrics_climate_de <- poll_metrics_wide %>%
left_join(
current_climate_de %>%
dplyr::select(plz, bio01, bio02, bio04, bio05, bio06,
bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio18),
by = "plz"
)Farmer acceptance scores (1–6 scale) from the survey are joined to the combined landscape-climate table.
survey_raw <- read_csv("/Users/jonasschreiber/Downloads/ergebnis_verknuepft.csv")
survey_clean <- survey_raw %>%
dplyr::select(PLZ, Akzeptanz_Mean, Bereitschaft) %>%
mutate(PLZ = as.character(PLZ))
survey_data_combined <- survey_clean %>%
left_join(metrics_climate_de, by = c("PLZ" = "plz"))
cat("Survey observations:", nrow(survey_data_combined), "\n")Variables are tested for correlation.
# Landscape metrics only
landscape_metrics <- metrics_climate_de %>%
dplyr::select(aggregation_index, cohesion, edge_density, mesh,
number_patches, patch_density, patch_richness, shdi,
split, largest_patch_index, enn_mn)
cor_landscape <- cor(landscape_metrics, use = "complete.obs", method = "spearman")
corrplot(cor_landscape, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black",
number.cex = 0.7, tl.cex = 0.8,
title = "Spearman correlation – landscape metrics")# Correlation of climate variables with acceptance
cor_climate <- survey_data_combined %>%
dplyr::select(Akzeptanz_Mean, bio01, bio02, bio04, bio05,
bio06, bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio18) %>%
cor(use = "complete.obs") %>%
as.data.frame() %>%
dplyr::select(Akzeptanz_Mean) %>%
arrange(desc(abs(Akzeptanz_Mean)))
print(cor_climate)Predictors are specified on the hypothesis that landscape strructure and climate influence farmers acceptance of pollinator friendly meassures. Before selection they are tested for collinearity.
VIF is run on the full theoretical candidate pool (all landscape
metrics + all bioclim variables available in survey_model)
before finalizing which variables to keep, rather than only
checking the final 4 after the fact.
exclude_cols <- c("Mittelwert_Bewertung", "Gewichteter_Index", "PLZ", "X", "Y",
"resid", "...1") # adjust to match actual non-predictor columns
candidate_vars <- c(
"aggregation_index", "cohesion", "edge_density", "enn_mn",
"largest_patch_index", "mesh", "number_patches", "patch_density",
"patch_richness", "shdi", "split",
"bio01", "bio02", "bio04", "bio05", "bio06", "bio10", "bio11",
"bio12", "bio13", "bio14", "bio15", "bio16", "bio18"
)
formula_candidate <- as.formula(
paste("Mittelwert_Bewertung ~", paste(candidate_vars, collapse = " + "))
)
m_candidate <- glm(formula_candidate, data = survey_model)
vif(m_candidate)Decision rule (pre-specified): variables with VIF > 10 are dropped iteratively — removing the variable with the highest VIF first, then re-checking remaining variables, until all are below the threshold. Where a choice exists between equally collinear variables, the one with stronger theoretical justification is retained (e.g. bio01 over bio12 for direct heat-stress relevance; shdi/edge_density retained together as they capture distinct compositional vs. structural landscape dimensions).
Note: with 18 candidate predictors and n = 91, this candidate model is underpowered (~5 obs/predictor) and is used purely as a VIF diagnostic — not as an analysis model in its own right.
After testing, the following variables are selected: -
Climate block: bio01 (mean annual
temperature), bio12 (annual precipitation) -
Landscape block: shdi (Shannon diversity
index — compositional diversity), edge_density
(configurational/structural complexity)
model_formula <- Mittelwert_Bewertung ~ bio01 + bio12 + edge_density + shdiSo the model is:
m_full <- glm(model_formula, data = survey_model)
summary(m_full)set.seed(42) # fixed seed for reproducibility
train_control <- trainControl(method = "repeatedcv", number = 10, repeats = 10)
cv_model <- train(
model_formula,
data = survey_model,
method = "glm",
trControl = train_control
)
print(cv_model)Reported performance metrics (RMSE, R-squared, MAE) are taken as the honest, locked-in estimate of predictive performance and are not used to iteratively modify the model.
# join PLZ centroid coordinates onto survey_model
plz_centroids <- st_centroid(plz_de)
coords_df <- plz_centroids %>%
st_coordinates() %>%
as.data.frame() %>%
mutate(plz = plz_de$plz)
# ensure consistent types before joining (avoids leading-zero / casing mismatches)
coords_df$plz <- as.character(coords_df$plz)
survey_model$PLZ <- as.character(survey_model$PLZ)
survey_model <- survey_model %>%
left_join(coords_df, by = c("PLZ" = "plz"))survey_model$resid <- residuals(m_full)
coords <- cbind(survey_model$X, survey_model$Y)
knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)
weights <- nb2listw(nb, style = "W")
moran.test(survey_model$resid, weights)Result: Moran’s I close to zero and non-significant indicates no evidence of spatial autocorrelation in residuals. The independence assumption underlying the GLM and CV results holds reasonably well for this sample.
Method follows the logic of MESS-style extrapolation checks used in species distribution modeling (Elith, Kearney & Phillips, 2010, Methods in Ecology and Evolution): extrapolation occurs when an environmental value falls outside the range present in the training data. This implementation checks each variable’s range independently (a simplified, univariate version of the full multivariate MESS approach).
training_ranges <- survey_model %>%
summarise(
bio01_min = min(bio01), bio01_max = max(bio01),
bio12_min = min(bio12), bio12_max = max(bio12)
)
training_rangesmetrics_climate_de <- metrics_climate_de %>%
mutate(
bio01_extrap_baseline = bio01 < training_ranges$bio01_min | bio01 > training_ranges$bio01_max,
bio12_extrap_baseline = bio12 < training_ranges$bio12_min | bio12 > training_ranges$bio12_max
)
mean(metrics_climate_de$bio01_extrap_baseline) * 100 # % of PLZs extrapolated, bio01
mean(metrics_climate_de$bio12_extrap_baseline) * 100 # % of PLZs extrapolated, bio12metrics_climate_de <- metrics_climate_de %>%
mutate(
bio01_gfdl_126_extrap = bio01_gfdl_126 < training_ranges$bio01_min | bio01_gfdl_126 > training_ranges$bio01_max,
bio01_gfdl_585_extrap = bio01_gfdl_585 < training_ranges$bio01_min | bio01_gfdl_585 > training_ranges$bio01_max,
bio01_ipsl_585_extrap = bio01_ipsl_585 < training_ranges$bio01_min | bio01_ipsl_585 > training_ranges$bio01_max,
bio12_gfdl_126_extrap = bio12_gfdl_126 < training_ranges$bio12_min | bio12_gfdl_126 > training_ranges$bio12_max,
bio12_gfdl_585_extrap = bio12_gfdl_585 < training_ranges$bio12_min | bio12_gfdl_585 > training_ranges$bio12_max,
bio12_ipsl_585_extrap = bio12_ipsl_585 < training_ranges$bio12_min | bio12_ipsl_585 > training_ranges$bio12_max
)
extrapolation_summary <- data.frame(
scenario = c("gfdl_126", "gfdl_585", "ipsl_585"),
pct_bio01_extrap = c(
mean(metrics_climate_de$bio01_gfdl_126_extrap) * 100,
mean(metrics_climate_de$bio01_gfdl_585_extrap) * 100,
mean(metrics_climate_de$bio01_ipsl_585_extrap) * 100
),
pct_bio12_extrap = c(
mean(metrics_climate_de$bio12_gfdl_126_extrap) * 100,
mean(metrics_climate_de$bio12_gfdl_585_extrap) * 100,
mean(metrics_climate_de$bio12_ipsl_585_extrap) * 100
)
)
extrapolation_summaryKey finding: bio01 extrapolation is already ~100% at baseline (current climate), before any future scenario is applied. This indicates a sampling/generalizability limitation. Survey respondents occupy a narrow temperature band relative to Germany as a whole, rather than an artifact of future climate change. bio12 (precipitation) extrapolation is low (~1-2%) both at baseline and under future scenarios, indicating the survey sample reasonably represents the national precipitation range.
For future surveys, more representative survey participants in wider climate zones should be chosen. In this case study we still show a German-wide map, giving a blueprint at german wide prediction, with the disclaimer that more surveys are necessary.
metrics_climate_de$predicted_acceptance <- predict(m_full, newdata = metrics_climate_de, type = "response")Landscape variables (edge_density, shdi)
are held constant; only climate variables are swapped in per
scenario.
predict_scenario <- function(data, bio01_col, bio12_col, model) {
temp_data <- data %>%
mutate(
bio01 = .data[[bio01_col]],
bio12 = .data[[bio12_col]]
)
predict(model, newdata = temp_data, type = "response")
}
metrics_climate_de$predicted_gfdl_126 <- predict_scenario(metrics_climate_de, "bio01_gfdl_126", "bio12_gfdl_126", m_full)
metrics_climate_de$predicted_gfdl_585 <- predict_scenario(metrics_climate_de, "bio01_gfdl_585", "bio12_gfdl_585", m_full)
metrics_climate_de$predicted_ipsl_585 <- predict_scenario(metrics_climate_de, "bio01_ipsl_585", "bio12_ipsl_585", m_full)This documment described the workflow to link farmers acceptance of pollinator friendly meassures to climate data and landscape metrics, based on a pollinator friendly map. The model was chosen on the Hypothesis that both climate, and landscape influence farmers acceptance. As predictors, Annual temperature, annual preciaption, Shannon-Diversity-Index and Edge Density were selected. The model was checked for collinearitry, spatial autocorrelation and extrapolation. After showing that more surveys in different climate zones are needed for a more generalized spatial prediction, it is still succesfull in linking farmers acceptance to climate and landscape data.
sessionInfo()