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)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.
### Comparison Across Data Sources
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.
``` r
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"
)
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)sessionInfo()