---
title: "Spatial PECs across EU countries"
output:
  html_notebook:
  toc: true
  toc_depth: 3
  toc_float: true
---

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(sf)
library(fs)
library(terra)
library(tidyterra)
library(tmap)
library(readxl)
library(scales)
library(ggprism)
library(ggrepel)
```

# Evaluation steps

## Data preparation

### Country sampling site coordinates

#### Local Administrative Units (LAUs)

```{r message=FALSE, warning=FALSE}
lau_degurba_cz <- path_wd("/Data/Shapes of NUTS/LAU_RG_01M_2024_4326.gpkg") |>
  vect() |> 
  tidyterra::filter(CNTR_CODE == "CZ")

lau_degurba_nl <- path_wd("/Data/Shapes of NUTS/LAU_RG_01M_2024_4326.gpkg") |> 
  vect() |> 
  tidyterra::filter(CNTR_CODE == "NL")

lau_degurba_dk <- path_wd("/Data/Shapes of NUTS/LAU_RG_01M_2024_4326.gpkg") |>
  vect() |> 
  tidyterra::filter(CNTR_CODE == "DK")
```

#### Sampling coordinates

```{r message=FALSE, warning=FALSE}

site_id_cz <- c("CZ_EF01_S_P", "CZ_EF02_S_P", "CZ_EF03_S_P", "CZ_EF04_S_P", "CZ_EF05_S_P", "CZ_EF06_S_P", "CZ_EF07_S_P", "CZ_EF08_S_P", "CZ_EF09_S_P", "CZ_EF10_S_P", "CZ_EF11_S_P", "CZ_EF12_S_P", "CZ_EF13_S_P", "CZ_EF14_S_P", "CZ_EF15_S_P", "CZ_EF16_S_P", "CZ_EF17_S_P", "CZ_EF18_S_P", "CZ_EF19_S_P", "CZ_EF20_S_P", "CZ_EF21_S_P", "CZ_EF22_S_P", "CZ_EF23_S_P", "CZ_EF24_S_P")

site_id_lat_cz <- c(48.75428, 48.7651997, 49.5882283, 48.9971858, 48.8974494, 48.9203217, 49.5573031, 49.7089406, 49.7190786, 50.5149308, 50.5392117, 50.4087597, 49.0598786, 49.0573903, 50.1286728, 50.2487525, 50.22116, 49.7082947, 48.9869569, 49.6619561, 49.7321864, 48.8276081, 48.8326775, 48.7364833)

site_id_lon_cz <- c(16.9401214, 16.9324394, 17.1860597, 16.9157722, 16.0395308, 15.9784944, 13.9029903, 14.8711725, 14.8851656, 16.196935, 16.1935875, 14.1494492, 15.4676047, 15.4676261, 16.2425311, 17.6263658, 15.9933244, 14.7895919, 17.0273617, 12.9795581, 12.9605036, 16.5171839, 16.4857483, 16.9331583)

site_id_nl <- c("NL_EF01_S_P", "NL_EF02_S_P", "NL_EF03_S_P", "NL_EF04_S_P", "NL_EF05_S_P", "NL_EF06_S_P", "NL_EF07_S_P", "NL_EF08_S_P", "NL_EF09_S_P", "NL_EF10_S_P", "NL_EF11_S_P", "NL_EF12_S_P", "NL_EF13_S_P", "NL_EF14_S_P", "NL_EF15_S_P", "NL_EF16_S_P")
site_id_lat_nl <- c(53.31, 53.32, 53.39, 53.38, 53.22, 53.4, 53.36, 53.39, 53.33, 53.28, 53.42, 53.35, 53.43, 53.38, 53.19, 53.32)
site_id_lon_nl <- c(6.27, 6.13, 6.37, 6.33, 5.47, 6.42, 6.41, 6.03, 6.32, 6.26, 6.6, 6.46, 6.63, 6.34, 5.46, 6.33)

site_id_dk <- c("DK_EF01_S_P", "DK_EF02_S_P", "DK_EF03_S_P", "DK_EF04_S_P", "DK_EF05_S_P", "DK_EF06_S_P", "DK_EF07_S_P", "DK_EF08_S_P", "DK_EF09_S_P", "DK_EF10_S_P", "DK_EF11_S_P", "DK_EF12_S_P", "DK_EF13_S_P", "DK_EF14_S_P", "DK_EF15_S_P", "DK_EF16_S_P", "DK_EF17_S_P", "DK_EF18_S_P", "DK_EF19_S_P", "DK_EF20_S_P")

site_id_lat_dk <- c(57, 57, 56.5, 56.51, 56.43, 56.44, 56.95, 56.97, 56.57, 56.57, 56.34, 56.35, 56.81, 56.81, 56.2, 56.23, 56.75, 56.78, 56.53,56.53)
site_id_lon_dk <- c(8.75, 8.73, 8.17, 8.28, 8.98, 8.94, 10.12, 10.16, 8.21, 8.19, 10.35, 10.37, 9.35, 9.35, 9.39, 9.4, 9.96, 9.89, 8.19, 8.21)

site_coord_cz <- cbind(tibble(
  site_id = site_id_cz,
  long = site_id_lon_cz,
  lat = site_id_lat_cz)) |>
  sf::st_as_sf(coords = c("long", "lat"), crs = st_crs(lau_degurba_cz)) |>
  vect() |>
  terra::intersect(lau_degurba_cz)

site_coord_nl <- cbind(tibble(
  site_id = site_id_nl,
  long = site_id_lon_nl,
  lat = site_id_lat_nl)) |> 
  sf::st_as_sf(coords = c("long", "lat"), crs = st_crs(lau_degurba_nl)) |>
  vect() |>
  terra::intersect(lau_degurba_nl)

site_coord_dk <- cbind(tibble(
  site_id = site_id_dk,
  long = site_id_lon_dk,
  lat = site_id_lat_dk)) |>
  sf::st_as_sf(coords = c("long", "lat"), crs = st_crs(lau_degurba_dk)) |>
  vect() |>
  terra::intersect(lau_degurba_dk)
```

##### LAUs with sampling sites
```{r message=FALSE, warning=FALSE}
lau_sample_cz <- lau_degurba_cz |> filter(LAU_NAME %in% site_coord_cz$LAU_NAME) |> distinct(LAU_NAME)
lau_sample_nl <- lau_degurba_nl |> filter(LAU_NAME %in% site_coord_nl$LAU_NAME) |> distinct(LAU_NAME)
lau_sample_dk <- lau_degurba_dk |> filter(LAU_NAME %in% site_coord_dk$LAU_NAME) |> distinct(LAU_NAME)
```

##### GEMUPs
```{r message=FALSE, warning=FALSE}
gemup_lau_cz <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "gemup100_slim_cz.gpkg", recursive = T)) |> 
      vect(extent = ext(lau_sample_cz)) |> 
      intersect(lau_sample_cz) |>
      mutate(District = str_to_title((gsub("_", " ", District, fixed = TRUE))),
             Active = str_to_title(Active),
             Crop = str_to_sentence(Crop)) |>
      # filter(!is.na(Crop)) |>
      rename(infactor_av = IFav,
             infactor_min = IFmin,
             infactor_max = IFmax,
             acsubst = Active,
             acsubst_mass_kg = ASmass,
             acsubst_area_ha = ASarea,
             field_area = FieldAr,
             aprate_farm_g.ha = ARfarm,
             aprate_min = ARmin,
             aprate_mmax = ARmax,
             apfreq = ApFreq) |>   
  makeValid()

crop_nl <- paste0(path_home_r(),
                  "/",
                  list.files(path_home_r(),
                             pattern = "brpgewaspercelen_definitief_2021",
                             recursive = T)) |> vect() |> project(crs(lau_degurba_nl))

crop_nl |> mask(lau_sample_nl[5]) |> plot("gewas")


gemup_lau_nl <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "gemap_nl", recursive = T)) |> 
      vect() |> 
      makeValid() |> 
      terra::intersect(lau_sample_nl) |> 
      rename(infactor_av = inter_frac_bbch.ave,
             infactor_min = inter_frac_bbch.min,
             infactor_max = inter_frac_bbch.max,
             acsubst = EU_name,
             field_area = area,
             aprate_farm_g.ha = Dose_kg_ha,
             Crop = EC_hcat_n) |> 
      mutate(field_area = field_area * 0.0001,
             Crop = str_replace_all(Crop, "_", " "),
             Crop = str_to_title(Crop))
 

 gemup_lau_dk <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "gemap_dk", recursive = T)) |>
      vect(extent = ext(lau_sample_dk)) |>
   makeValid() |> 
   terra::intersect(lau_sample_dk) |> 
      mutate(Active = str_to_sentence(Active)) |> 
      rename(infactor_av = inter_frac_bbch.ave,
             infactor_min = inter_frac_bbch.min,
             infactor_max = inter_frac_bbch.max,
             acsubst = Active,
             field_area = field_size,
             aprate_farm_kg.ha = App_dose_kg.ha,
             Crop = EC_hcat_n) |> 
    mutate(Crop = str_replace_all(Crop, "_", " "),
           Crop = str_to_title(Crop))
```
###### Crops data

```{r}
crop_as_cz <- gemup_lau_cz |> values() |> group_by(LAU_NAME) |> reframe(n_distinct(Crop), n_distinct(acsubst))
crop_as_nl <- gemup_lau_nl |> values() |> group_by(LAU_NAME) |> reframe(n_distinct(Crop), n_distinct(acsubst))
crop_as_dk <- gemup_lau_dk |> values() |> group_by(LAU_NAME) |> reframe(n_distinct(Crop), n_distinct(acsubst))

bind_rows(crop_as_cz, crop_as_nl, crop_as_dk)
```


## Application rates evaluation data
```{r message=FALSE, warning=FALSE}
aprat_nl <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "_PEC_Soil_NL\\.[csv]", recursive = T)) |>
  map(\(x) read_csv(x)) |>
  bind_rows() |> 
  add_column(Country = "NL")

aprat_cz <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "_PEC_Soil_CZ\\.[csv]", recursive = T)) |>
  map(\(x) read_csv(x)) |>
  bind_rows() |> 
  add_column(Country = "CZ")

aprat_dk <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "_PEC_Soil_DK\\.[csv]", recursive = T)) |>
  map(\(x) read_csv(x)) |>
  bind_rows() |> 
  add_column(Country = "DK")

aprat_gemup <- bind_rows(aprat_nl, aprat_cz, aprat_dk)

aprat_mark <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "Mark\\s\\w+\\s\\w+\\s\\d+\\w+", recursive = T)) |> 
  read_xlsx()

aprate_df <- right_join(
aprat_mark[c("Active Substance", "kg.ha.app")],
aprat_gemup[c("acsubst",
              "aprate_farm_kg.ha",
              "Country")] |>
  distinct(acsubst, Country, .keep_all = T),
by = c("Active Substance" = "acsubst")) |> 
  rename(AS = "Active Substance",
        ar_gemup_kg.ha = aprate_farm_kg.ha,
        ar_sprint_hg.ha = kg.ha.app)
```

## PEC soil evaluation data
### Observed dataset
```{r message=FALSE, warning=FALSE}
# MEC soil Knuth et al. 2024
soil_mec_cz <- read_excel(paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "Knuth soil sampling", recursive = T)), range = "A1:GO65") |> 
  mutate(across(6:197, ~as.numeric(.))) |> 
  rename(Sample_code = `SPRINT Sample code`) |> 
  pivot_longer(cols = -c(Sample_code, Country, Management, Region, Crop),
               names_to = "Active",
               values_to = "MEC_µg.kg") |> 
  filter(!is.na(MEC_µg.kg), Country == "CZ")

soil_mec_nl <- read_excel(paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "Knuth soil sampling", recursive = T)), range = "A1:GO65") |> 
  mutate(across(6:197, ~as.numeric(.))) |> 
  rename(Sample_code = `SPRINT Sample code`) |> 
  pivot_longer(cols = -c(Sample_code, Country, Management, Region, Crop),
               names_to = "Active",
               values_to = "MEC_µg.kg") |> 
  filter(!is.na(MEC_µg.kg), Country == "NL") |> 
  mutate(Sample_code = str_remove(Sample_code, "_Field\\s\\d"))

soil_mec_dk <- read_excel(paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "Knuth soil sampling", recursive = T)), range = "A1:GO65") |> 
  mutate(across(6:197, ~as.numeric(.))) |> 
  rename(Sample_code = `SPRINT Sample code`) |> 
  pivot_longer(cols = -c(Sample_code, Country, Management, Region, Crop),
               names_to = "Active",
               values_to = "MEC_µg.kg") |> 
  filter(!is.na(MEC_µg.kg), Country == "DK")
```

### Modelled dataset
```{r message=FALSE, warning=FALSE}

soil_pec_cz_path <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "PEC_Soil_CZ.csv", recursive = T))
soil_pec_cz <- map(soil_pec_cz_path, read_csv) |>
  bind_rows() |>
  rename(Active = acsubst) |> 
  filter(str_detect(Crop, "rape"))

soil_pec_nl_path <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "PEC_Soil_NL.csv", recursive = T))
soil_pec_nl <- map(soil_pec_nl_path, read_csv) |>
  bind_rows() |>
  rename(Active = acsubst) |> 
  filter(str_detect(Crop, "Potat"))

soil_pec_nl |> filter(Active == "Fluopicolide")

winter_cereals <- c("Grain maize corn popcorn", "Winter barley", "Winter barley", "Winter rapeseed rape", "Winter rye", "Winter spelt", "Winter triticale")
spring_cereals <- c("Grain maize corn popcorn", "Spring barley", "Spring common soft wheat", "Spring oats", "Spring rapeseed rape", "Spring rye", "Spring triticale")

soil_pec_dk_path <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "PEC_Soil_DK.csv", recursive = T))
soil_pec_dk <- map(soil_pec_dk_path, read_csv) |>
  bind_rows() |>
  rename(Active = acsubst) |> 
  filter(Crop %in% spring_cereals)

```

### Evaluation dataset
```{r}
eval_data_pec_cz <- soil_pec_cz |>
  group_by(LAU_NAME, Active) |> 
  summarise(soil_conc_ini_med = mean(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_ini_sd = sd(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_56d_med = mean(conc_acsubst_total_soil_56twa_ug.kg),
            soil_conc_56d_sd = sd(conc_acsubst_total_soil_56twa_ug.kg),
            organic_carbon = mean(oc_perc*100),
            organic_carbon = sd(oc_perc*100),
            density_bulk = mean(bulk_dens_kg.dm3),
            density_bulk = sd(bulk_dens_kg.dm3),
            clay = mean(clay_perc),
            clay = sd(clay_perc),
            sand = mean(sand_perc),
            sand = sd(sand_perc),
            slope = mean(slope_perc),
            slope = sd(slope_perc),
            nr_field = n()) |>
  bind_cols("Concentration units" = "µg × kg⁻¹") |> 
  rename("Location" = LAU_NAME,
         "Active substance" = Active,
         "Concentration intial" = soil_conc_ini_med,
         "Concentration intial SD" = soil_conc_ini_sd,
         "Concentration 56 days" = soil_conc_56d_med,
         "Concentration 56 days SD" = soil_conc_56d_sd,
         "# of fields" = nr_field)

eval_data_mec_cz <- soil_mec_cz |>
  left_join(site_coord_cz |> values(), by = join_by("Sample_code" == "site_id")) |> 
  summarise(soil_conc_sample = max(MEC_µg.kg),
            nr_samples = n(),
            .by = c(LAU_NAME, Active, Management, Country)) |>
  rename("Location" = LAU_NAME,
         "Concentration sample" = soil_conc_sample,
         "Active substance" = Active)

eval_data_pec_nl <- soil_pec_nl |>
  group_by(LAU_NAME, Active) |> 
  summarise(soil_conc_ini_med = mean(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_ini_sd = sd(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_56d_med = mean(conc_acsubst_total_soil_56twa_ug.kg),
            soil_conc_56d_sd = sd(conc_acsubst_total_soil_56twa_ug.kg),
            organic_carbon = mean(oc_perc*100),
            organic_carbon = sd(oc_perc*100),
            density_bulk = mean(bulk_dens_kg.dm3),
            density_bulk = sd(bulk_dens_kg.dm3),
            clay = mean(clay_perc),
            clay = sd(clay_perc),
            sand = mean(sand_perc),
            sand = sd(sand_perc),
            slope = mean(slope_perc),
            slope = sd(slope_perc),
            nr_field = n()) |>
  bind_cols("Concentration units" = "µg × kg⁻¹") |> 
  rename("Location" = LAU_NAME,
         "Active substance" = Active,
         "Concentration intial" = soil_conc_ini_med,
         "Concentration intial SD" = soil_conc_ini_sd,
         "Concentration 56 days" = soil_conc_56d_med,
         "Concentration 56 days SD" = soil_conc_56d_sd,
         "# of fields" = nr_field) 

eval_data_mec_nl <- soil_mec_nl |>
  left_join(site_coord_nl |> values(), by = join_by("Sample_code" == "site_id")) |> 
  # group_by(LAU_NAME, Active, Management) |> 
  summarise(soil_conc_sample = max(MEC_µg.kg),
            nr_samples = n(),
            .by = c(LAU_NAME, Active, Management, Country)) |>
  rename("Location" = LAU_NAME,
         "Concentration sample" = soil_conc_sample,
         "Active substance" = Active)

eval_data_pec_dk <- soil_pec_dk |>
  group_by(LAU_NAME, Active) |> 
  summarise(soil_conc_ini_med = mean(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_ini_sd = sd(conc_acsubst_total_soil_ini_ug.kg),
            soil_conc_56d_med = mean(conc_acsubst_total_soil_56twa_ug.kg),
            soil_conc_56d_sd = sd(conc_acsubst_total_soil_56twa_ug.kg),
            organic_carbon = mean(oc_perc*100),
            organic_carbon = sd(oc_perc*100),
            density_bulk = mean(bulk_dens_kg.dm3),
            density_bulk = sd(bulk_dens_kg.dm3),
            clay = mean(clay_perc),
            clay = sd(clay_perc),
            sand = mean(sand_perc),
            sand = sd(sand_perc),
            slope = mean(slope_perc),
            slope = sd(slope_perc),
            nr_field = n()) |>
  bind_cols("Concentration units" = "µg × kg⁻¹") |> 
  rename("Location" = LAU_NAME,
         "Active substance" = Active,
         "Concentration intial" = soil_conc_ini_med,
         "Concentration intial SD" = soil_conc_ini_sd,
         "Concentration 56 days" = soil_conc_56d_med,
         "Concentration 56 days SD" = soil_conc_56d_sd,
         "# of fields" = nr_field) 

eval_data_mec_dk <- soil_mec_dk |>
  left_join(site_coord_dk |> values(), by = join_by("Sample_code" == "site_id")) |> 
  summarise(soil_conc_sample = max(MEC_µg.kg),
            nr_samples = n(),
            .by = c(LAU_NAME, Active, Management, Country)) |>
  rename("Location" = LAU_NAME,
         "Concentration sample" = soil_conc_sample,
         "Active substance" = Active)

eval_data_df_cz <- left_join(eval_data_mec_cz, eval_data_pec_cz , by = c("Location", "Active substance")) |> na.omit()
write_excel_csv(eval_data_df_cz, "CZ soil concentration comparison.csv")

eval_data_df_nl <- left_join(eval_data_mec_nl, eval_data_pec_nl , by = c("Location", "Active substance")) |> na.omit()
write_excel_csv(eval_data_df_nl, "NL soil concentration comparison.csv")

eval_data_df_dk <- left_join(eval_data_mec_dk, eval_data_pec_dk, by = c("Location", "Active substance")) |> na.omit()
write_excel_csv(eval_data_df_dk, "DK soil concentration comparison.csv")
```

## PEC water evaluation data
```{r message=FALSE, warning=FALSE}
water_pec_nl_path <- paste0(path_home_r(), "/",
                            list.files(path_home_r(),
                                       pattern = "Harlingen\\S+PEC_Water_NL\\.csv",
                                       recursive = T))
  
river_class <- basename(water_pec_nl_path[1:5]) |> str_extract("\\d\\S+_") |> str_remove("_PEC_Water_")
river_buff100 <- water_pec_nl_path[1] |> str_extract("\\w+_\\d+") |> str_remove("Water_")
river_buff20 <- water_pec_nl_path[6] |> str_extract("\\w+_\\d+") |> str_remove("Water_")

data <- water_pec_nl_path[1] |> read_csv() |> add_column("class" = river_class[1], "buffer_m" = river_buff100)
data$conc_river_ug.dm3 |> hist()
```

# Plots 
## Sampling sites
```{r}
# Sample site map
sampling_sites_map <- function(lau, site_coord){

  tmap_mode("plot")
  basemap_lau <- tm_shape(lau) +
    tm_borders(col = "black",
               lty = "solid",
               lwd = 0.25) +
    tm_basemap("Esri.WorldTopoMap", alpha = 0.85) +
    # tmap_options(component.autoscale = T) +
    tm_scalebar(position = c("left", "bottom"))
  
  layout_map <- tm_layout(legend.outside = T,
                          legend.outside.position = "right",
                          legend.outside.size = 0.5,
                          legend.text.size = 0.5,
                          legend.title.size = 0.75,
                          main.title.fontface = "bold",
                          scale = 1, 
                          # main.title.size = 0.5
                          )
ss.map <- basemap_lau +
    layout_map +
    tm_shape(lau |> filter(LAU_NAME %in% site_coord$LAU_NAME)) +
    tm_polygons(fill = "LAU_NAME",
                fill.scale = tm_scale(values = "poly.alphabet2"),
                fill.legend = tm_legend("LAU (2024)", frame = F),
                fill_alpha = 0.75) +
    tm_shape(site_coord["site_id"]) +
    tm_symbols(fill = "magenta",
               size = .3,) +
    tm_add_legend(type = "symbols",
                  fill = "magenta",
                  size = 0.5,
                  labels = "Sampling site (2021)",
                  frame = F) + 
    tm_add_legend(labels = "LAU border",
                  type = "lines") +
    tm_title(paste0("Soil sampling locations in ",
                    lau$CNTR_CODE |> unique(), 
                     " (# of sites ",
                     site_coord$site_id |> unique() |> length(),
                     ")"), size = 1)
  
  tmap_save(ss.map, paste0(lau$CNTR_CODE |> unique(), "_ss.png"), width = 8, height = 4, dpi = 300)
}

sampling_sites_map(lau_degurba_nl, site_coord_nl)
```
## Crop maps LAU

```{r}
crop_map <- function(gemup, lau, site_coord){

  lau_name <- lau$LAU_NAME
  
  for(i in seq_along(lau_name)) {
  
  tmap_mode("plot")
  basemap_lau <- tm_shape(lau |> filter(LAU_NAME %in% lau_name[i])) +
    tm_borders(col = "black",
               lty = "solid",
               lwd = 0.15) +
    tm_basemap("Esri.WorldTopoMap", alpha = 0.85) +
    # tmap_options(component.autoscale = T) +
    tm_scalebar(position = c("left", "bottom"))
 
  layout_map <- tm_layout(legend.outside = T,
                          legend.outside.position = "right",
                          legend.outside.size = 0.5,
                          legend.text.size = 0.5,
                          legend.title.size = 0.75,
                          main.title.fontface = "bold",
                          scale = 1, 
                          # main.title.size = 0.5
                          )
crop.map <- basemap_lau +
    layout_map +
  tm_shape(gemup |> mask(lau |> filter(LAU_NAME %in% lau_name[i]))) +
  tm_polygons(fill = "Crop",
                fill.scale = tm_scale_categorical(values = "poly.alphabet2", n.max = 40),
                fill.legend = tm_legend("Crops (2021)", frame = F),
                fill_alpha = 0.75,
              lwd = 0.15) +
    tm_shape(lau |> filter(LAU_NAME %in% lau_name[i])) +
    tm_borders() +
  tm_shape(site_coord["site_id"] |> mask(lau|> filter(LAU_NAME %in% lau_name[i]))) +
    tm_symbols(fill = "magenta",
               size = .5,) +
    tm_add_legend(labels = "LAU border",
                  type = "lines") +
   tm_add_legend(type = "symbols",
                  fill = "magenta",
                  size = 0.5,
                  labels = "Sampling site",
                  frame = F) + 
    tm_title(lau |> values() |> filter(LAU_NAME %in% lau_name[i]),
                    size = 1)
  
  tmap_save(crop.map, paste0(lau_name[i],"_crop.png"), width = 8, height = 4, dpi = 300)
  }
}

crop_map(gemup_lau_nl, lau_sample_nl[5], site_coord_nl)
```
## Application rates
```{r}
aprate.plt <- ggplot(aprate_df) +
  geom_point(aes(x = ar_gemup_kg.ha, y = ar_sprint_hg.ha, colour = Country),
             shape = 21,
             size = 2,
             stroke = 0.9,
             alpha = 0.85) +
  geom_text_repel(aes(x = ar_gemup_kg.ha,
                      y = ar_sprint_hg.ha,
                      label = AS,
                      fontface = "bold",
                      segment.size = 0.25,
                      segment.linetype = "dotted",
                      segment.colour = "gray25"),
                  size = 2.5,
                  max.overlaps = 30) +
  geom_abline(aes(intercept = 0, slope = 1, linetype = "1:1 line"), linewidth = 0.25) +
  geom_abline(aes(intercept = 0.4043007, slope = 1, linetype = "±1sd line"), linewidth = 0.25) +
  geom_abline(aes(intercept = -0.4043007, slope = 1, linetype = "±1sd line"), linewidth = 0.25) +
  # geom_line(data = data.frame(x = c(0, 2), y = c(0, 2)), aes(x = x, y = y, linetype = "1:1 line"), lwd = 0.25) +
  # geom_line(data = data.frame(x = c(0, 2), y = c(0.5, 2.5)), aes(x = x, y = y, linetype = "±.5× line"), lwd = 0.25) +
  # geom_line(data = data.frame(x = c(0.5, 2), y = c(0, 1.5)), aes(x = x, y = y, linetype = "±.5× line"), lwd = 0.25) +
  scale_linetype_manual(
      name = "Reference lines",
      values = c("1:1 line" = "solid", 
                 # "±2 OoM (100×)" = "dotted",
                 "±1sd line" = "dashed")) +
  scale_color_manual(name = "Case Study countries",
                     labels = c("Czechia (oilseed rape)",
                                "Denmark (cereals)",
                                "Netherlands (potatos)"),
                     values = c("steelblue", "seagreen", "tomato")) +
  labs(title = "Surveyed vs estimated application rates across case studies, 2021",
       x = "Estimated application rate [kg\u00D7ha⁻¹]",
       y = "Surveyed application rate [kg\u00D7ha⁻¹]") +
  theme(
    aspect.ratio = 1,
    plot.title = element_text(hjust = 0.5, face = "bold", size = 7),
    legend.position = "right",
    legend.box = "vertical",
    legend.spacing.y = unit(0.0, "cm"),
    panel.grid = element_blank()) +
  theme_bw(base_size = 13) +
  coord_cartesian(xlim = c(-0.5, 2), ylim = c(-0.5,2))
  
ggsave("Aprate.png", 
       aprate.plt,
       width = 9,
       height = 7,
       dpi = 300)
```
## PEC soil
### Bar plot
```{r}
# Read the data
nl_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "NL soil concentration comparison.csv", recursive = T)) |> read_csv()
cz_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "CZ soil concentration comparison.csv", recursive = T)) |> read_csv()
dk_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "DK soil concentration comparison.csv", recursive = T)) |> read_csv()

pec_mec_bar.plt <- function(df){
# Prepare data for plotting
soil_conc_long <- df |> 
  # Create location label with number of fields
  mutate(Location_label = paste0(Location, "\n(n=", `# of fields`, ")")) |> 
  pivot_longer(
    cols = c(`Concentration sample`, `Concentration intial`, `Concentration 56 days`),
    names_to = "Measurement_type",
    values_to = "Concentration"
  ) |> 
  mutate(
  SD = case_when(
    Measurement_type == "Concentration intial" ~ "Concentration intial SD",
    Measurement_type == "Concentration 56 days" ~ "Concentration 56 days SD"),
  # Clean up measurement type labels
  Measurement_type = case_when(
    Measurement_type == "Concentration sample" ~ "Sample (measured)",
    Measurement_type == "Concentration intial" ~ "Initial (modeled)",
    Measurement_type == "Concentration 56 days" ~ "56-day TWA (modeled)"
  )) |> 
  mutate(
  # Reorder factor levels for plotting
    Measurement_type = factor(Measurement_type, 
                              levels = c("Sample (measured)", 
                                         "Initial (modeled)", 
                                         "56-day TWA (modeled)")),
    Location_label = factor(Location_label, 
                            levels = unique(Location_label[order(Location, `Active substance`)]))) |> 
  mutate(conc_ini_lsd = Concentration - `Concentration intial SD`,
         conc_ini_lsd = case_when(Measurement_type == "56-day TWA (modeled)" ~ NA,
                                  Measurement_type == "Sample (measured)" ~ NA,
                                  .default = conc_ini_lsd),
         conc_ini_usd = Concentration + `Concentration intial SD`,
         conc_ini_usd = case_when(Measurement_type == "56-day TWA (modeled)" ~ NA,
                                  Measurement_type == "Sample (measured)" ~ NA,
                                  .default = conc_ini_usd),
         conc_twa_lsd = Concentration - `Concentration 56 days SD`,
         conc_twa_lsd = case_when(Measurement_type == "Initial (modeled)" ~ NA,
                                  Measurement_type == "Sample (measured)" ~ NA,
                                  .default = conc_twa_lsd),
         conc_twa_usd = Concentration + `Concentration 56 days SD`,
         conc_twa_usd = case_when(Measurement_type == "Initial (modeled)" ~ NA,
                                  Measurement_type == "Sample (measured)" ~ NA,
                                  .default = conc_twa_usd))

  unique_locations <- soil_conc_long %>%
    filter(Management == "Conventional") %>%
    distinct(Location, `# of fields`) |> 
    mutate(location_num = row_number()) %>%
    ungroup()
  
  vline_data <- unique_locations %>%
    group_by(Location) %>%
    mutate(xintercept = location_num + 0.5) %>%
    ungroup()
 

plt <- ggplot(soil_conc_long, 
         aes(x = Location_label, y = Concentration, fill = Measurement_type)) +
    geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8) +
    geom_errorbar(aes(ymin = conc_ini_lsd , ymax = conc_ini_usd),
                  position = position_dodge(width = 0.9),
                  width = 0.25,
                  na.rm = TRUE) +
    geom_errorbar(aes(ymin = conc_twa_lsd , ymax = conc_twa_usd),
                  position = position_dodge(width = 0.9),
                  width = 0.25,
                  na.rm = TRUE) +
    facet_wrap(~ `Active substance`, scales = "free_y", ncol = 1) +
    scale_fill_manual(values = c("Sample (measured)" = "#E69F00", 
                                 "Initial (modeled)" = "#56B4E9", 
                                 "56-day TWA (modeled)" = "#009E73")) +
    geom_vline(data = vline_data[-14,], 
               aes(xintercept = xintercept), 
               linetype = "dashed", 
               color = "gray50", 
               linewidth = 0.5) +
    labs(
      title = paste0("Comparison of Measured and Modeled Pesticide Concentrations in Soil in ", df$Country |> unique()),
      subtitle = "Error bars represent standard deviation (SD) for modeled values",
      x = "Location (n = number of fields)",
      y = expression(paste("Concentration (", mu, "g ", " kg"^{-1}, ")")),
      fill = "Measurement Type"
    ) +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
      legend.position = "bottom",
      legend.title = element_text(face = "bold"),
      strip.text = element_text(face = "bold", size = 12),
      strip.background = element_rect(fill = "lightgray", color = NA),
      plot.title = element_text(face = "bold", size = 14),
      plot.subtitle = element_text(size = 10, color = "gray40")
    )

  ggsave(paste0(df$Country |> unique(),"_MPEC_soil_bar.png"), 
         plt, width = 12, height = 8, dpi = 300)
}
pec_mec_bar.plt(dk_conc)
```
### Scatter plot
```{r message=FALSE, warning=FALSE}
# Read the data
nl_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "NL soil concentration comparison.csv", recursive = T)) |> read_csv()
cz_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "CZ soil concentration comparison.csv", recursive = T)) |> read_csv()
dk_conc <- paste0(path_home_r(), "/", list.files(path_home_r(), pattern = "DK soil concentration comparison.csv", recursive = T)) |> read_csv()

# IMPORTANT: These colors are fixed to ensure consistency across all plots
# Each substance always gets the same color regardless of which dataset it appears in
# Colors are organized by pesticide type for better interpretability

substance_colors <- c(
  'Azoxystrobin' = '#4d4d4d',
  'Dimoxystrobin' = '#de77ae',   
  'Difenoconazole' = '#b2df8a',    
  'Fluopicolide' = '#33a02c',     
  'Fluopyram' = '#fb9a99',      
  'Mandipropamid' = '#e31a1c',  
  'Tebuconazole' = '#b2182b',   
  'Glyphosate' = '#ff7f00',       
  'Metobromuron' = '#cab2d6',      
  'Pendimethalin' = '#5ab4ac',    
  'Prosulfocarb' = '#6a3d9a',      
  'Deltamethrin' = '#636363',
  "Boscalid" = "#b15928",
  "Diflufenican" = "#2166ac"
)

# ============================================================================
# Define markers (pch values) for locations
# ============================================================================
location_markers_cz <- c(
  'Kostice' = 9,           
  'Kravsko' = 18,           
  'Krumvíř' = 17,            
  'Olomouc' = 16,             
  'Postupice' = 15,        
  'Velké Petrovice' = 6,   
  'Vlašim' = 5,         
  'Vševily' = 2,        
  'Znojmo' = 1,             
  'Česká Metuje' = 0      
)

location_markers_nl <- c(
  "Noardeast-Fryslân" = 5,
  'Harlingen' = 2,       
  'Het Hogeland' = 1,       
  'Westerkwartier' = 0    
)

location_markers_dk <- c(
  "Syddjurs" = 18,
  "Silkeborg" = 17,
  "Thisted" = 5,
  "Vesthimmerlands" = 6,
  "Mariagerfjord" = 15,
  "Aalborg" = 2,
  "Holstebro" = 1,
  "Lemvig" = 0   
)
# ============================================================================
# Function to create scatterplot for a specific timepoint
# ============================================================================

pec_mec_scatter.plt <- function(data, country, timepoint = "initial") {
  
  # Select which concentration to use
  if (timepoint == "initial") {
    data <- data %>%
      mutate(
        Model = `Concentration intial`,
        Model_SD = `Concentration intial SD`
      )
  } else {
    data <- data %>%
      mutate(
        Model = `Concentration 56 days`,
        Model_SD = `Concentration 56 days SD`
      )
  }
  
plt <- ggplot(data = data, aes(x = `Concentration sample`, y = Model)) +
    geom_point(
      aes(color = `Active substance`, shape = Location, fill = `Active substance`), 
      size = 3,
      alpha = 0.7,
      stroke = 0.8
    ) +
  geom_abline(aes(intercept = 0, slope = 1, linetype = "1:1 line"), linewidth = 0.25) +
  geom_abline(aes(intercept = log10(10) , slope = 1, linetype = "±1 OoM (10×)"), linewidth = 0.25) +
  geom_abline(aes(intercept = -log10(10), slope = 1, linetype = "±1 OoM (10×)"), linewidth = 0.25) +
  # geom_line(data = data.frame(x = c(1, 100000), y = c(1, 100000)), aes(x = x, y = y, linetype = "1:1 line")) +
  # geom_line(data = data.frame(x = c(1, 10000), y = c(10, 100000)), aes(x = x, y = y, linetype = "±1 OoM (10×)")) +
  # geom_line(data = data.frame(x = c(10, 100000), y = c(1, 10000)), aes(x = x, y = y, linetype = "±1 OoM (10×)")) +
  # geom_line(data = data.frame(x = c(1, 100000), y = c(100, 100000)), aes(x = x, y = y, linetype = "±2 OoM (100×)")) +
  #   geom_line(data = data.frame(x = c(100, 100000), y = c(1, 1000)), aes(x = x, y = y, linetype = "±2 OoM (100×)")) +
    scale_linetype_manual(
      name = "Reference lines",
      values = c("1:1 line" = "solid", 
                 # "±2 OoM (100×)" = "dotted",
                 "±1 OoM (10×)" = "dashed")) +
    # geom_errorbar(aes(ymin = Model - Model_SD,
    #                   ymax = Model + Model_SD),
    #               color = "black",
    #               width = 0.075,
    #               linewidth = 0.3,
    #               alpha = 0.8) +
    # Scales
    # scale_color_manual(values = substance_colors, aesthetics = c("colour", "fill")) +
  scale_colour_brewer(palette = "Dark2", type = "div") +
  scale_fill_brewer(palette = "Dark2", type = "div") +
    scale_shape_manual(values = if(country == "Czechia"){
      location_markers_cz
      } else if (country == "Netherlands") {
        location_markers_nl
        } else {
          location_markers_dk
        }) +
    labs(x = expression(paste("Measured concentration (µg\u00D7", kg^{-1}, ")")),
         y = expression(paste("Modelled concentration (µg\u00D7", kg^{-1}, ")")),
         title = paste0(country, ": ", timepoint, " soil concentration"),
         color = "Active substance",
         shape = "Location"
    ) +
    scale_x_log10(breaks = breaks_log(n = 6, 10),
                  limits = c(1e-1, 1e4),
                  minor_breaks = rep(1:9, 5)*(10^rep(0:4, each = 9)),
                  guide = "prism_minor",
                  labels = label_log()) +
    scale_y_log10(breaks = breaks_log(n = 6, 10),
                  limits = c(1e-1, 1e4),
                  minor_breaks = rep(1:9, 5)*(10^rep(0:4, each = 9)),
                  guide = "prism_minor",
                  labels = label_log()) +
    theme_bw(base_size = 13) +
    theme(
      aspect.ratio = 1,
      plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
      legend.position = "right",
      legend.box = "vertical",
      legend.spacing.y = unit(0.0, "cm"),
      panel.grid = element_blank())
  
  print(plt)
  # Save plot
  ggsave(paste0(data$Country |> unique(), "_MPEC_soil_", timepoint ,"_scatter.png"), 
         plt, width = 7, height = 7, dpi = 300)

}

pec_mec_scatter.plt(dk_conc, "Denmark", "initial")

```

