---
title: "Script_2_Spatial_data_integration_fungalradar"
author: "Hylke Kortenbosch"
date: "2024-03-01"
output: html_document
---

Load packages
```{r}
# List of required packages
packages <- c("sf", "gstat", "stars", "ggplot2", "tidyverse", "terra", "tmap", "ncdf4", "raster")

# Function to check, install (if needed), and load packages
install_and_load <- function(pkg) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}

# Apply the function to each package
lapply(packages, install_and_load)
```

# Key session info
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] svglite_2.1.3      performance_0.12.4 gridExtra_2.3      terra_1.7-71       lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
[12] tibble_3.2.1       tidyverse_2.0.0    stars_0.6-5        abind_1.4-5        MASS_7.3-60.2      gstat_2.1-1        tmap_3.3-4         ggspatial_1.1.9    ggplot2_3.5.1      MuMIn_1.47.5       glmmTMB_1.1.9     
[23] sf_1.0-16 

#Loading and cleaning up data, generate Figure S1 (histogram of sample exposure times) 
```{r}
# read data
Schimmelradar360 <- read.csv("Schimmelradar_360_submission.csv")
#subset data to remove samples that had well over 300 or below 25 counts for RF estimations
Schimmelradar360 <- Schimmelradar360[Schimmelradar360$Note != "Exclude",]

# add median spore CFUs per trap
Schimmelradar360 <- Schimmelradar360 %>%
  rowwise() %>%
  mutate(median_CFU = median(c(Tot.ITR, Tot.VOR, Tot.CON)))

#Make dates out of relevant variables in df:
Schimmelradar360$Arrived..back.on  <- as.Date(Schimmelradar360$Arrived..back.on , format = "%d/%m/%Y")
Schimmelradar360$Start.sampling  <- as.Date(Schimmelradar360$Start.sampling , format = "%d/%m/%Y")
Schimmelradar360$End.sampling  <- as.Date(Schimmelradar360$End.sampling , format = "%d/%m/%Y")

# Normalise and round median CFU counts to 28 days exposure
Schimmelradar360$Days_exposed <- as.numeric(Schimmelradar360$End.sampling - Schimmelradar360$Start.sampling)
#Plot histogram of days exposed for 
hist(Schimmelradar360$Days_exposed, main = "Histogram of number of days participants exposed the traps", xlab = "Days Exposed")
Schimmelradar360$norm_median_CFU <-round(Schimmelradar360$median_CFU/(Schimmelradar360$Days_exposed/28))

hist(Schimmelradar360$norm_median_CFU)
```

Make summary table of the Dutch population data (Table 1)
```{r}
Population_summary <- data.frame(Schimmelradar360$norm_median_CFU, Schimmelradar360$Tot.VOR, Schimmelradar360$RES.VOR, Schimmelradar360$RF.VOR, Schimmelradar360$Tot.ITR, Schimmelradar360$RES.ITR,  Schimmelradar360$RF.ITR)


totals <- colSums(Population_summary) # ignore totals of resistance fractions here
means <- colMeans(Population_summary)
medians <- apply(Population_summary, 2, median)


# Standard error = Standard deviation / sqrt(n)
std_error.norm.median.CFU <- sd(Schimmelradar360$norm_median_CFU) / sqrt(length(Schimmelradar360$norm_median_CFU)) # Standard error = Standard deviation / sqrt(n)
std_error.Tot.VOR <- sd(Schimmelradar360$Tot.VOR) / sqrt(length(Schimmelradar360$Tot.VOR))
std_error.Res.VOR <- sd(Schimmelradar360$RES.VOR) / sqrt(length(Schimmelradar360$RES.VOR))
std_error.RF.VOR <- sqrt(mean(Schimmelradar360$RF.VOR) * (1 - mean(Schimmelradar360$RF.VOR)) / length(Schimmelradar360$Tot.VOR)) # SE proportions = sqrt(p*(1-p)/n)

std_error.Tot.ITR <- sd(Schimmelradar360$Tot.ITR) / sqrt(length(Schimmelradar360$Tot.ITR))
std_error.Res.ITR <- sd(Schimmelradar360$RES.ITR) / sqrt(length(Schimmelradar360$RES.ITR))
std_error.RF.ITR <- sqrt(mean(Schimmelradar360$RF.ITR) * (1 - mean(Schimmelradar360$RF.ITR)) / length(Schimmelradar360$Tot.ITR)) 

std_errors <- c(std_error.norm.median.CFU, std_error.Tot.VOR, std_error.Res.VOR, std_error.RF.VOR, std_error.Tot.ITR, std_error.Res.ITR, std_error.RF.ITR)



# Create the summary table
summary_table <- data.frame(Totals = totals, Mean_per_Sample = means,  Standard_Error = std_errors, Median_of_dataset = medians)

# Print the summary table
print(summary_table)
```

Make data spatial
```{r}
Schimmelradar360spat_round <- st_as_sf(Schimmelradar360, coords=c("Round_long", "Round_lat"), 
                    crs=4326)

# transform to Dutch grid
Schimmelradar360spat_round <- st_transform(Schimmelradar360spat_round, crs=28992)
str(Schimmelradar360spat_round)
```


Make empty grid for spatial prediction
```{r}
# make grid to fill in with prediction
data(NLD_prov)
provincesNL <- sf::st_sf(NLD_prov)

bbox <- st_bbox(provincesNL_transformed_dummy)

# Print the bounding box
print(bbox)
provincesNL_transformed_dummy <- st_transform(provincesNL, crs = 4326)

# convert to new style CRS and Dutch grid
st_crs(provincesNL) <- st_crs(provincesNL)
provincesNL_transformed <- st_transform(provincesNL, crs = 28992)
```

Load LGN22 raster file (Land use of the Netherlands for 2022) dataset. Can be downloaded as geopackage from: https://lgn.nl/bestanden
```{r}
lgn22 <- rast("LGN2022.tif")
```

Load data from Basisregistratie Gewaspercelen (BRP) ATOM
This is an underlying source of the LGN database: https://nationaalgeoregister.nl/geonetwork/srv/dut/catalog.search#/metadata/44e6d4d3-8fc5-47d6-8712-33dd6d244eef
```{r}
#Gewassen2022 <- st_read("brpgewaspercelen_definitief_2022.gpkg") # this file is very large and can by downloaded from the link above

#legenda_gewas_2022 <- as.data.frame(cbind(as.vector(unique(Gewassen2022$gewas)), as.vector(unique(Gewassen2022$gewascode))))

#save(legenda_gewas_2022, file = "legenda_gewas_2022.Rdata")
#load("legenda_gewas_2022.Rdata")
#view(legenda_gewas_2022) # note that this legend is in Dutch, if required it can be downloaded as a csv and translated

#create a rasterfile projected the same as lgn22, maintain the codes for the different crops from the brp dataset
#rastgewas22 <- rasterize(vect(Gewassen2022), lgn22, field="gewascode", filename = "Gewas22rast.tiff")
rastgewas22 <- rast("Gewas22rast.tiff")  #simply load this to run the code below
```


# reclassification & indicator coding of variables of interest at High definition and subsequencet agregatation in 250m raster to make the focal computation of variables through space less intense. 
```{r}
water22HD2 <- app(lgn22, function(x) {ifelse(x %in% c(16,17), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "water22HD2.tif", overwrite = T)
water22HD2 <- rast("water22HD2.tif")

water22_250 <- aggregate(water22HD2, fact=50, fun="mean",
                            filename="water22_250_2.tif")
water22_250 <- rast("water22_250_2.tif")

saltwater22HD2 <- app(lgn22, function(x) {ifelse(x %in% c(17), 
                                             1, 
                                             ifelse(NA))},
              filename= "saltwater22HD.tif", overwrite = T)
saltwater22HD2 <- rast("saltwater22HD.tif")

saltwater22_250 <- aggregate(saltwater22HD2, fact=50, fun="mean",
                            filename="saltwater22_250.tif", overwrite = T)
saltwater22_250 <- rast("saltwater22_250.tif")
plot(saltwater22_250)

grass22HD <- app(lgn22, function(x) {ifelse(x %in% c(1, 23, 28, 45, 46, 47), 
                                                  1, 
                                                  ifelse(is.na(x), NA, 0))},
                  filename="grass22HD.tif", overwrite = T)
grass22HD <- rast("grass22HD.tif")

grass22_250 <- aggregate(grass22HD, fact=50, fun="mean",
                            filename="grass22_250.tif", overwrite = T)
grass22_250 <- rast("grass22_250.tif")

maize22HD <- app(lgn22, function(x) {ifelse(x %in% c(2), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename="maize22HD.tif", overwrite = T)
maize22HD <- rast("maize22HD.tif")

maize22_250 <- aggregate(maize22HD, fact=50, fun="mean",
                            filename="maize22_250.tif")
maize22_250 <- rast("maize22_250.tif")
plot(maize22_250)

potato22HD <- app(lgn22, function(x) {ifelse(x %in% c(3), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename="potato22HD.tif", overwrite = T)
potato22HD <- rast("potato22HD.tif")

potato22_250 <- aggregate(potato22HD, fact=50, fun="mean",
                            filename="potato22_250.tif")
potato22_250 <- rast("potato22_250.tif")
plot(potato22_250)

seed_potato_22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(2015, 2016), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "seed_potato22HD2.tif", overwrite = T)
seed_potato_22HD2 <- rast("seed_potato22HD2.tif")
seed_potato_22HD2[is.na(seed_potato_22HD2)] <- 0

seed_potato22HD_NA <- mask(seed_potato_22HD2, vect(provincesNL_transformed), filename="seed_potato22HD2_NA.tif", overwrite =T)

seed_potato22_250 <- aggregate(seed_potato_22HD2, fact=50, fun="mean",
                            filename="seed_potato22_250.tif", overwrite = T)
seed_potato22_250 <- rast("seed_potato22_250.tif")
plot(seed_potato22_250)

starch_potato22HD <- app(rastgewas22, function(x) {ifelse(x %in% c(2017), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "starch_potato22HD.tif", overwrite = T)
starch_potato22HD <- rast("starch_potato22HD.tif")
starch_potato22HD[is.na(starch_potato22HD)] <- 0

starch_potato22HD_NA <- mask(starch_potato22HD, vect(provincesNL_transformed), filename="starch_potato22HD_NA.tif", overwrite =T)

starch_potato22_250 <- aggregate(starch_potato22HD, fact=50, fun="mean",
                            filename="starch_potato22_2502.tif", overwrite = T)
starch_potato22_250 <- rast("starch_potato22_2502.tif")
plot(starch_potato22_250)

consumption_potato22HD <- app(rastgewas22, function(x) {ifelse(x %in% c(2014), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "consumption_potato22HD.tif", overwrite = T)
consumption_potato22HD <- rast("consumption_potato22HD.tif")
consumption_potato22HD[is.na(consumption_potato22HD)] <- 0

consumption_potato22HD_NA <- mask(consumption_potato22HD, vect(provincesNL_transformed), filename="consumption_potato22HD_NA.tif", overwrite =T)

consumption_potato22_250 <- aggregate(consumption_potato22HD, fact=50, fun="mean",
                            filename="consumption_potato22_2502.tif", overwrite = T)
consumption_potato22_250 <- rast("consumption_potato22_2502.tif")
plot(consumption_potato22_250)

beetroot22HD <- app(lgn22, function(x) {ifelse(x %in% c(4), 
                                                  1, 
                                                  ifelse(is.na(x), NA, 0))},
                   filename="beetroot22HD.tif", overwrite = T)
beetroot22HD <- rast("beetroot22HD.tif")

beetroot22_250 <- aggregate(beetroot22HD, fact=50, fun="mean",
                            filename="beetroot22_250.tif")
beetroot22_250 <- rast("beetroot22_250.tif")

beetroot_potato22HD <- app(lgn22, function(x) {ifelse(x %in% c(3,4), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename="beetroot_potato22HD.tif", overwrite = T)
beetroot_potato22HD <- rast("beetroot_potato22HD.tif")

beetroot_potato22_250 <- aggregate(beetroot_potato22HD, fact=50, fun="mean",
                            filename="beetroot_potato22_250.tif")
beetroot_potato22_250 <- rast("beetroot_potato22_250.tif")
plot(beetroot_potato22_250)

grain22HD <- app(lgn22, function(x) {ifelse(x %in% c(5), 
                                               1, 
                                               ifelse(is.na(x), NA, 0))},
                filename="grain22HD.tif", overwrite = T)
grain22HD <- rast("grain22HD.tif")

grain22_250 <- aggregate(grain22HD, fact=50, fun="mean",
                            filename="grain22_250.tif")
grain22_250 <- rast("grain22_250.tif")

greenhouse22HD <- app(lgn22, function(x) {ifelse(x %in% c(8), 
                                               1, 
                                               ifelse(is.na(x), NA, 0))},
                filename="greenhouse22HD.tif", overwrite = T)
greenhouse22HD <- rast("greenhouse22HD.tif")

greenhouse22_250 <- aggregate(greenhouse22HD, fact=50, fun="mean",
                            filename="greenhouse22_250.tif")
greenhouse22_250 <- rast("greenhouse22_250.tif")

buildings22HD <- app(lgn22, function(x) {ifelse(x %in% c(18, 19, 24, 26, 251), 
                                            1, 
                                            ifelse(is.na(x), NA, 0))},
             filename="buildings22HD.tif", overwrite = T)
buildings22HD <- rast("buildings22HD.tif")

buildings22_250 <- aggregate(buildings22HD, fact=50, fun="mean",
                            filename="buildings22_250.tif")
buildings22_250 <- rast("buildings22_250.tif")


heather22HD <- app(lgn22, function(x) {ifelse(x %in% c(34, 36, 37, 38), 
                                            1, 
                                            ifelse(is.na(x), NA, 0))},
             filename="heather22HD.tif", overwrite = T)
heather22HD <- rast("heather22HD.tif")

heather22_250 <- aggregate(heather22HD, fact=50, fun="mean",
                            filename="heather22_250.tif")
heather22_250 <- rast("heather22_250.tif")

Peatland22HD <- app(lgn22, function(x) {ifelse(x %in% c(39), 
                                            1, 
                                            ifelse(is.na(x), NA, 0))},
             filename="Peatland22HD.tif", overwrite = T)
Peatland22HD <- rast("Peatland22HD.tif")

Peatland22_250 <- aggregate(Peatland22HD, fact=50, fun="mean",
                            filename="Peatland22_250.tif")
Peatland22_250 <- rast("Peatland22_250.tif")

forest22HD <- app(lgn22, function(x) {ifelse(x %in% c(11, 12, 20, 22, 40, 43), 
                                    1, 
                                    ifelse(is.na(x), NA, 0))},
           filename="forest22HD.tif", overwrite = T)
forest22HD <- rast("forest22HD.tif")

forest22_250 <- aggregate(forest22HD, fact=50, fun="mean",
                            filename="forest22_250.tif")
forest22_250 <- rast("forest22_250.tif")

Othercrops22HD <- app(lgn22, function(x) {ifelse(x %in% c(6), 
                                      1, 
                                      ifelse(is.na(x), NA, 0))},
           filename="Othercrops22HD.tif", overwrite = T)
Othercrops22HD <- rast("Othercrops22HD.tif")

Othercrops22_250 <- aggregate(Othercrops22HD, fact=50, fun="mean",
                            filename="Othercrops22_250.tif")
Othercrops22_250 <- rast("Othercrops22_250.tif")
plot(Othercrops22_250)

onion22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(6664, 1934, 1933, 6660, 263), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "onion22HD2.tif", overwrite = T)
onion22HD2 <- rast("onion22HD2.tif")
onion22HD2[is.na(onion22HD2)] <- 0

onion22HD_NA <- mask(onion22HD2, vect(provincesNL_transformed), filename="onion22HD2_NA.tif", overwrite =T)

onion22_250 <- aggregate(onion22HD2, fact=50, fun="mean",
                            filename="onion22_250.tif", overwrite = T)
onion22_250 <- rast("onion22_250.tif")
plot(onion22_250)

seed_onion22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(6664, 6660), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "seed_onion22HD2.tif", overwrite = T)
seed_onion22HD2 <- rast("seed_onion22HD2.tif")
seed_onion22HD2[is.na(seed_onion22HD2)] <- 0

seed_onion22HD_NA <- mask(seed_onion22HD2, vect(provincesNL_transformed), filename="seed_onion22HD2_NA.tif", overwrite =T)

seed_onion22_250 <- aggregate(seed_onion22HD2, fact=50, fun="mean",
                            filename="seed_onion22_250.tif", overwrite = T)
seed_onion22_250 <- rast("seed_onion22_250.tif")
plot(seed_onion22_250)

set_onion22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(1934, 1933, 263), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "set_onion22HD2.tif", overwrite = T)
set_onion22HD2 <- rast("set_onion22HD2.tif")
set_onion22HD2[is.na(set_onion22HD2)] <- 0

set_onion22HD_NA <- mask(set_onion22HD2, vect(provincesNL_transformed), filename="set_onion22HD2_NA.tif", overwrite =T)

set_onion22_250 <- aggregate(set_onion22HD2, fact=50, fun="mean",
                            filename="set_onion22_250.tif", overwrite = T)
set_onion22_250 <- rast("set_onion22_250.tif")
plot(set_onion22_250)

fruits22HD <- app(lgn22, function(x) {ifelse(x %in% c(9, 62), 
                                            1, 
                                            ifelse(is.na(x), NA, 0))},
             filename="fruits22HD.tif", overwrite = T)
fruits22HD <- rast("fruits22HD.tif")

fruit22_250 <- aggregate(fruits22HD, fact=50, fun="mean",
                            filename="fruit22_250.tif")
fruit22_250 <- rast("fruit22_250.tif")
plot(fruit22_250)

apple22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(1095, 1096), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "apple22HD2.tif", overwrite = T)
apple22HD2 <- rast("apple22HD2.tif")
apple22HD2[is.na(apple22HD2)] <- 0

apple22HD_NA <- mask(apple22HD2, vect(provincesNL_transformed), filename="apple22HD2_NA.tif", overwrite =T)

apple22_250 <- aggregate(apple22HD2, fact=50, fun="mean",
                            filename="apple22_250.tif", overwrite = T)
apple22_250 <- rast("apple22_250.tif")
plot(apple22_250)

apple_pear22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(1095, 1096, 1097, 1098), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "apple_pear22HD2.tif", overwrite = T)
apple_pear22HD2 <- rast("apple_pear22HD2.tif")
apple_pear22HD2[is.na(apple_pear22HD2)] <- 0

apple_pear22HD_NA <- mask(apple_pear22HD2, vect(provincesNL_transformed), filename="apple_pear22HD2_NA.tif", overwrite =T)

apple_pear22_250 <- aggregate(apple_pear22HD2, fact=50, fun="mean",
                            filename="apple_pear22_250.tif", overwrite = T)
apple_pear22_250 <- rast("apple_pear22_250.tif")
plot(apple_pear22_250)

pear22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(1097, 1098), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "pear22HD2.tif", overwrite = T)
pear22HD2 <- rast("pear22HD2.tif")
pear22HD2[is.na(pear22HD2)] <- 0

pear22HD_NA <- mask(pear22HD2, vect(provincesNL_transformed), filename="pear22HD2_NA.tif", overwrite =T)

pear22_250 <- aggregate(pear22HD2, fact=50, fun="mean",
                            filename="pear22_250.tif", overwrite = T)
pear22_250 <- rast("pear22_250.tif")
plot(pear22_250)

tree_farming22HD <- app(lgn22, function(x) {ifelse(x %in% c(61), 
                                            1, 
                                            ifelse(is.na(x), NA, 0))},
             filename="tree_farming22HD.tif", overwrite = T)
tree_farming22HD <- rast("tree_farming22HD.tif")

tree_farming22_250 <- aggregate(tree_farming22HD, fact=50, fun="mean",
                            filename="tree_farming22_250.tif")
tree_farming22_250 <- rast("tree_farming22_250.tif")



bulbs22HD <- app(lgn22, function(x) {ifelse(x %in% c(10), 
                                           1, 
                                           ifelse(is.na(x), NA, 0))},
                filename="bulbs22HD.tif", overwrite = T)
bulbs22HD <- rast("bulbs22HD.tif")

bulbs22_250 <- aggregate(bulbs22HD, fact=50, fun="mean",
                            filename="bulbs22_250.tif")

bulbs22_250 <- rast("bulbs22_250.tif")

tulip22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(985, 1004), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "tulip22HD2.tif", overwrite = T)
tulip22HD2 <- rast("tulip22HD2.tif")
tulip22HD2[is.na(tulip22HD2)] <- 0

tulip22HD_NA <- mask(tulip22HD2, vect(provincesNL_transformed), filename="tulip22HD2_NA.tif", overwrite =T)

tulip22_250 <- aggregate(tulip22HD2, fact=50, fun="mean",
                            filename="tulip22_250.tif", overwrite = T)
tulip22_250 <- rast("tulip22_250.tif")
plot(tulip22_250)


lily22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(979, 980, 1002), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "lily22HD2.tif", overwrite = T)
lily22HD2 <- rast("lily22HD2.tif")
lily22HD2[is.na(lily22HD2)] <- 0

lily22HD_NA <- mask(lily22HD2, vect(provincesNL_transformed), filename="lily22HD2_NA.tif", overwrite =T)

lily22_250 <- aggregate(lily22HD2, fact=50, fun="mean",
                            filename="lily22_250.tif", overwrite = T)
lily22_250 <- rast("lily22_250.tif")
plot(lily22_250)

gladiolus22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(967, 968, 998), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "gladiolus22HD2.tif", overwrite = T)
gladiolus22HD2 <- rast("gladiolus22HD2.tif")
gladiolus22HD2[is.na(gladiolus22HD2)] <- 0

gladiolus22HD_NA <- mask(gladiolus22HD2, vect(provincesNL_transformed), filename="gladiolus22HD2_NA.tif", overwrite =T)

gladiolus22_250 <- aggregate(gladiolus22HD2, fact=50, fun="mean",
                            filename="gladiolus22_250.tif", overwrite = T)
gladiolus22_250 <- rast("gladiolus22_250.tif")
plot(gladiolus22_250)

daffodil22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(982, 983, 1003), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "daffodil22HD2.tif", overwrite = T)
daffodil22HD2 <- rast("daffodil22HD2.tif")
daffodil22HD2[is.na(daffodil22HD2)] <- 0

daffodil22HD_NA <- mask(daffodil22HD2, vect(provincesNL_transformed), filename="daffodil22HD2_NA.tif", overwrite =T)

daffodil22_250 <- aggregate(daffodil22HD2, fact=50, fun="mean",
                            filename="daffodil22_250.tif", overwrite = T)
daffodil22_250 <- rast("daffodil22_250.tif")
plot(daffodil22_250)

hyacinth22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(970, 999, 1018), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "hyacinth22HD2.tif", overwrite = T)
hyacinth22HD2 <- rast("hyacinth22HD2.tif")
hyacinth22HD2[is.na(hyacinth22HD2)] <- 0

hyacinth22HD_NA <- mask(hyacinth22HD2, vect(provincesNL_transformed), filename="hyacinth22HD2_NA.tif", overwrite =T)

hyacinth22_250 <- aggregate(hyacinth22HD2, fact=50, fun="mean",
                            filename="hyacinth22_250.tif", overwrite = T)
hyacinth22_250 <- rast("hyacinth22_250.tif")
plot(hyacinth22_250)

daffodil_hyacinth22HD2 <- app(rastgewas22, function(x) {ifelse(x %in% c(970, 999, 1018, 982, 983, 1003), 
                                             1, 
                                             ifelse(is.na(x), NA, 0))},
              filename= "daffodil_hyacinth22HD2.tif", overwrite = T)
daffodil_hyacinth22HD2 <- rast("daffodil_hyacinth22HD2.tif")
daffodil_hyacinth22HD2[is.na(daffodil_hyacinth22HD2)] <- 0

daffodil_hyacinth22HD_NA <- mask(daffodil_hyacinth22HD2, vect(provincesNL_transformed), filename="daffodil_hyacinth22HD2_NA.tif", overwrite =T)

daffodil_hyacinth22_250 <- aggregate(daffodil_hyacinth22HD2, fact=50, fun="mean",
                            filename="daffodil_hyacinth22_250.tif", overwrite = T)
daffodil_hyacinth22_250 <- rast("daffodil_hyacinth22_250.tif")
plot(daffodil_hyacinth22_250)
```

Set focal filter for 10k radius circle
```{r}
wndw <- matrix(rep(1, 81^2), nrow=81)

wndw <- ifelse((col(wndw)-41)^2 +  (row(wndw)-41)^2 <= 1600, 1, 0)

# apply the focal filter
water22fra <- focal(water22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="water_fract2.tif")
plot(water22fra)
water22fra <- rast("water_fract2.tif")

grass22fra <- focal(grass22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="grass_fract.tif", overwrite=TRUE)
plot(grass22fra)
grass22fra <- rast("grass_fract.tif")

maize22fra <- focal(maize22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="maize_fract.tif")
plot(maize22fra)
maize22fra <- rast("maize_fract.tif")

potato22fra <- focal(potato22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="potato_fract.tif")
plot(potato22fra)
potato22fra <- rast("potato_fract.tif")

beetroot_potato22fra <- focal(beetroot_potato22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="beetroot_potato_fract.tif")
plot(beetroot_potato22fra)
beetroot_potato22fra <- rast("beetroot_potato_fract.tif")

seed_potato22fra <- focal(seed_potato22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="seed_potato_fract.tif")
seed_potato22fra <- rast("seed_potato_fract.tif")
plot(seed_potato22fra)

starch_potato22fra <- focal(starch_potato22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="starch_potato_fract.tif")
starch_potato22fra <- rast("starch_potato_fract.tif")
plot(starch_potato22fra)

consumption_potato22fra <- focal(consumption_potato22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="consumption_potato_fract.tif")
consumption_potato22fra <- rast("consumption_potato_fract.tif")
plot(consumption_potato22fra)

beetroot22fra <- focal(beetroot22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="beetroot_fract.tif")
plot(beetroot22fra)
beetroot22fra <- rast("beetroot_fract.tif")

grain22fra <- focal(grain22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="grain_fract.tif")
plot(grain22fra)
grain22fra <- rast("grain_fract.tif")

greenhouse22fra <- focal(greenhouse22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="greenhouse_fract.tif")
plot(greenhouse22fra)
greenhouse22fra <- rast("greenhouse_fract.tif")

buildings22fra <- focal(buildings22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="buildings_fract.tif", overwrite =T)
plot(buildings22fra)
buildings22fra <- rast("buildings_fract.tif")

forest22fra <- focal(forest22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="forest_fract.tif")
plot(forest22fra)
forest22fra <- rast("forest_fract.tif")

Onion22fra <- focal(onion22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="Onion_fract.tif", overwrite = T)
plot(Onion22fra)
Onion22fra <- rast("Onion_fract.tif")

seed_onion22fra <- focal(seed_onion22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="seed_onion_fract.tif", overwrite = T)
plot(seed_onion22fra)
seed_onion22fra <- rast("seed_onion_fract.tif")

set_onion22fra <- focal(set_onion22_250, w=wndw, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="set_onion_fract.tif", overwrite = T)
plot(set_onion22fra)
set_onion22fra <- rast("set_onion_fract.tif")

Othercrops22fra <- focal(Othercrops22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="Othercrops_fract.tif")
plot(Othercrops22fra)
Othercrops22fra <- rast("Othercrops_fract.tif")

fruit22fra <- focal(fruit22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="fruit_fract.tif")
plot(fruit22fra)
fruit22fra <- rast("fruit_fract.tif")

apple_pear22fra <- focal(apple_pear22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="apple_pear_fract.tif")
plot(apple_pear22fra)
apple_pear22fra <- rast("apple_pear_fract.tif")

apple22fra <- focal(apple22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="apple_fract.tif")
plot(apple22fra)
apple22fra <- rast("apple_fract.tif")

pear22fra <- focal(pear22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="pear_fract.tif")
plot(pear22fra)
pear22fra <- rast("pear_fract.tif")


tree_farming22fra <- focal(tree_farming22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="tree_farming_fract.tif")
plot(tree_farming22fra)
tree_farming22fra <- rast("tree_farming_fract.tif")

heather22fra <- focal(heather22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="heather_fract.tif", overwrite = T)
plot(heather22fra)
heather22fra <- rast("heather_fract.tif")

Peatland22fra <- focal(Peatland22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="Peatland_fract.tif", overwrite = T)
plot(Peatland22fra)
Peatland22fra <- rast("Peatland_fract.tif")

bulbs22fra <- focal(bulbs22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="bulbs_fract.tif")
plot(bulbs22fra)
bulb22fra <- rast("bulbs_fract.tif")

tulip22fra <- focal(tulip22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="tulip_fract.tif")
plot(tulip22fra)
tulip22fra <- rast("tulip_fract.tif")

lily22fra <- focal(lily22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="lily_fract.tif")
plot(lily22fra)
lily22fra <- rast("lily_fract.tif")

daffodil22fra <- focal(daffodil22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="daffodil_fract.tif")
plot(daffodil22fra)
daffodil22fra <- rast("daffodil_fract.tif")


hyacinth22fra <- focal(hyacinth22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="hyacinth_fract.tif")
plot(hyacinth22fra)
hyacinth22fra <- rast("hyacinth_fract.tif")

daffodil_hyacinth22fra <- focal(daffodil_hyacinth22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="daffodil_hyacinth_fract.tif",overwrite=TRUE)
plot(daffodil_hyacinth22fra)
daffodil_hyacinth22fra <- rast("daffodil_hyacinth_fract.tif")

gladiolus22fra <- focal(gladiolus22_250, w=wndw, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="gladiolus_fract.tif")
plot(gladiolus22fra)
gladiolus22fra <- rast("gladiolus_fract.tif")

```

Specify 10 km buffer around each fungal radar air sample location
```{r}
# buffer the data points @ 10 km radius
Schimmelradar360spat_focal <- Schimmelradar360spat_round
datbuf_focal <- st_buffer(Schimmelradar360spat_focal, dist=10000)

# check coordinate system; this LGN version is in WGS84, therefore we also
# unproject the disks to WGS84 (faster than transforming the raster)
bufWGS84_focal <- st_transform(datbuf_focal, crs= 4326)
```

Set focal filter for 5k radius circle
```{r}
wndw5 <- matrix(rep(1, 41^2), nrow=41)

wndw5 <- ifelse((col(wndw5)-21)^2 +  (row(wndw5)-21)^2 <= 800, 1, 0)

# apply the focal filter
water22fra5k <- focal(water22_250, w=wndw5, fun="mean", na.policy="omit", expand=TRUE, na.rm=TRUE, filename="water_fract5k.tif")
plot(water22fra5k)
water22fra5k <- rast("water_fract5k.tif")

grass22fra5k <- focal(grass22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="grass_fract5k.tif", overwrite =T)
plot(grass22fra5k)
grass22fra5k <- rast("grass_fract5k.tif")

maize22fra5k <- focal(maize22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="maize_fract5k.tif")
plot(maize22fra5k)
maize22fra5k <- rast("maize_fract5k.tif")

beetroot_potato22fra5k <- focal(beetroot_potato22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="beetroot_potato_fract5k.tif")
plot(beetroot_potato22fra5k)
beetroot_potato22fra5k <- rast("beetroot_potato_fract5k.tif")

potato22fra5k <- focal(potato22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="potato_fract5k.tif")
plot(potato22fra5k)
potato22fra5k <- rast("potato_fract5k.tif")

seed_potato22fra5k <- focal(seed_potato22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="seed_potato_fract5k.tif")
plot(seed_potato22fra5k)
seed_potato22fra5k <- rast("seed_potato_fract5k.tif")

starch_potato22fra5k <- focal(starch_potato22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="starch_potato_fract5k.tif")
plot(starch_potato22fra5k)
starch_potato22fra5k <- rast("starch_potato_fract5k.tif")

consumption_potato22fra5k <- focal(consumption_potato22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="consumption_potato_fract5k.tif")
plot(consumption_potato22fra5k)
consumption_potato22fra5k <- rast("consumption_potato_fract5k.tif")

beetroot22fra5k <- focal(beetroot22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="beetroot_fract5k.tif")
plot(beetroot22fra5k)
beetroot22fra5k <- rast("beetroot_fract5k.tif")

grain22fra5k <- focal(grain22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="grain_fract5k.tif")
plot(grain22fra5k)
grain22fra5k <- rast("grain_fract5k.tif")

greenhouse22fra5k <- focal(greenhouse22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="greenhouse_fract5k.tif")
plot(greenhouse22fra5k)
greenhouse22fra5k <- rast("greenhouse_fract5k.tif")

buildings22fra5k <- focal(buildings22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="buildings_fract5k.tif", overwrite =T)
plot(buildings22fra5k)
buildings22fra5k <- rast("buildings_fract5k.tif")

forest22fra5k <- focal(forest22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="forest_fract5k.tif")
plot(forest22fra5k)
forest22fra5k <- rast("forest_fract5k.tif")

Onion22fra5k <- focal(onion22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="Onion_fract5k.tif", overwrite = T)
plot(Onion22fra5k)
Onion22fra5k <- rast("Onion_fract5k.tif")

seed_onion22fra5k <- focal(seed_onion22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="seed_onion_fract5k.tif", overwrite = T)
plot(seed_onion22fra5k)
seed_onion22fra5k <- rast("seed_onion_fract5k.tif")

set_onion22fra5k <- focal(set_onion22_250, w=wndw5, fun="mean", na.policy="omit",

                     expand=TRUE, na.rm=TRUE, filename="set_onion_fract5k.tif", overwrite = T)
plot(set_onion22fra5k)
set_onion22fra5k <- rast("set_onion_fract5k.tif")

Othercrops22fra5k <- focal(Othercrops22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="Othercrops_fract5k.tif")
plot(Othercrops22fra5k)
Othercrops22fra5k <- rast("Othercrops_fract5k.tif")

fruit22fra5k <- focal(fruit22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="fruit_fract5k.tif", overwrite =T)
plot(fruit22fra5k)
fruit22fra5k <- rast("fruit_fract5k.tif")

apple_pear22fra5k <- focal(apple_pear22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="apple_pear_fract5k.tif")
plot(apple_pear22fra5k)

apple_pear22fra5k <- rast("apple_pear_fract5k.tif")

apple22fra5k <- focal(apple22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="apple_fract5k.tif")
plot(apple22fra5k)
apple22fra5k <- rast("apple_fract5k.tif")

pear22fra5k <- focal(pear22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="pear_fract5k.tif")
plot(pear22fra5k)
pear22fra5k <- rast("pear_fract5k.tif")

tree_farming22fra5k <- focal(tree_farming22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="tree_farming_fract5k.tif")
plot(tree_farming22fra5k)
tree_farming22fra5k <- rast("tree_farming_fract5k.tif")

heather22fra5k <- focal(heather22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="heather_fract5k.tif", overwrite = T)
plot(heather22fra5k)
heather22fra5k <- rast("heather_fract5k.tif")

Peatland22fra5k <- focal(Peatland22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="Peatland_fract5k.tif", overwrite = T)
plot(Peatland22fra5k)
Peatland22fra5k <- rast("Peatland_fract5k.tif")

bulbs22fra5k <- focal(bulbs22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="bulbs_fract5k.tif")
plot(bulbs22fra5k)
bulb22fra5k <- rast("bulbs_fract5k.tif")

tulip22fra5k <- focal(tulip22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="tulip_fract5k.tif")
plot(tulip22fra5k)
tulip22fra5k <- rast("tulip_fract5k.tif")

lily22fra5k <- focal(lily22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="lily_fract5k.tif")
plot(lily22fra5k)
lily22fra5k <- rast("lily_fract5k.tif")

daffodil22fra5k <- focal(daffodil22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="daffodil_fract5k.tif")
plot(daffodil22fra5k)
daffodil22fra5k <- rast("daffodil_fract5k.tif")


hyacinth22fra5k <- focal(hyacinth22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="hyacinth_fract5k.tif")
plot(hyacinth22fra5k)
hyacinth22fra5k <- rast("hyacinth_fract5k.tif")

gladiolus22fra5k <- focal(gladiolus22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="gladiolus_fract5k.tif")
plot(gladiolus22fra5k)
gladiolus22fra5k <- rast("gladiolus_fract5k.tif")

daffodil_hyacinth22fra5k <- focal(daffodil_hyacinth22_250, w=wndw5, fun="mean", na.policy="omit",

                      expand=TRUE, na.rm=TRUE, filename="daffodil_hyacinth_fract5k.tif")
plot(daffodil_hyacinth22fra5k)
daffodil_hyacinth22fra5k <- rast("daffodil_hyacinth_fract5k.tif")
```

Specify 5 km buffer around schimmelradar sample location
```{r}
# buffer the data points @ 5 km radius
Schimmelradar360spat_focal <- Schimmelradar360spat_round
datbuf_focal5k <- st_buffer(Schimmelradar360spat_focal, dist=5000)

# check coordinate system; this LGN version is in WGS84, therefore we also
# unproject the disks to WGS84 (faster than transforming the raster)
bufWGS84_focal5k <- st_transform(datbuf_focal5k, crs=4326)
```


Adding distance to the sea (salt water)
```{r}
DistanceSeaRast <- distance(saltwater22_250, filename="Distance_to_Sea.tif", overwrite =T)
plot(saltwater22_250)
# set NA beyond boundaries with NL mask
DistanceSeaRastNA <- mask(DistanceSeaRast, vect(provincesNL_transformed), filename="DistanceSeaRastNA.tif", overwrite =T)

DistanceSeaRastNA <- rast("DistanceSeaRastNA.tif")

plot(DistanceSeaRastNA) # Looks good!
```

Extracting Weatherdata 1-10-23 till 12-11-23 for schimmelradar analysis
source: https://cds.climate.copernicus.eu/datasets/sis-agrometeorological-indicators?tab=download
Subset values for NL data: West: 3.2, East: 7.4, South: 50.5, North: 53.7
# Windspeed data extraction and integration
```{r}
#Extracting windspeed data from ERA5
wind.nc <- list.files("Weather_schimmelradar/Wind_speed_schimmelradar", full.names = T)
wind_stack <- raster::stack(wind.nc)
mean_wind_schimmelrader <- mean(wind_stack)
str(mean_wind_schimmelrader)
mean_wind_schimmelrader@data@values[c(104, 143, 144, 145, 146)] <- mean_wind_schimmelrader@data@values[c(105, 185, 186, 187, mean(147, 188))] # plugging the gap of Vlieland missing values by copying values from neighboring cells, in the case of the one empty cell with two neighboring cells we take the mean of the two

sum(is.na(mean_wind_schimmelrader@data@values))

writeRaster(mean_wind_schimmelrader, "mean_wind_schimmelrader.nc", format = "CDF",overwrite =T)
spat_mean_wind_schimmelrader <- rast(mean_wind_schimmelrader)

spat_mean_wind_schimmelrader_250m <- terra::disagg(spat_mean_wind_schimmelrader, 40)
spat_mean_wind_schimmelrader_250m
plot(spat_mean_wind_schimmelrader_250m)

spat_mean_wind_schimmelrader_250m <- project(spat_mean_wind_schimmelrader_250m, DistanceSeaRastNA)
writeRaster(spat_mean_wind_schimmelrader_250m, "spat_mean_wind_schimmelrader_250m.tif", overwrite = TRUE)

spat_mean_wind_schimmelrader_250m <- rast("spat_mean_wind_schimmelrader_250m.tif")

Schimmelradar360spat_focal$wind_speed <- terra::extract(spat_mean_wind_schimmelrader_250m, Schimmelradar360spat_focal)[,2]

mock_raster_wind_NL <- mask(spat_mean_wind_schimmelrader_250m, vect(provincesNL_transformed), filename="Mock_wind_raster_NL.tif", overwrite = T)
plot(mock_raster_wind_NL)
```

# Temperature data extraction and integration
```{r}
#Extracting temperature data from AGERA5
temp.nc <- list.files("Weather_schimmelradar/2m_temp_schimmelradar", full.names = T)
temp_stack <- raster::stack(temp.nc)
mean_temp_schimmelrader <- mean(temp_stack)
str(mean_temp_schimmelrader)
mean_temp_schimmelrader@data@values[c(104, 143, 144, 145, 146)] <- mean_temp_schimmelrader@data@values[c(105, 185, 186, 187, mean(147, 188))] # plugging the gap of Vlieland missing values by copying values from neighboring cells, in the case of the one empty cell with two neighboring cells we take the mean of the two
writeRaster(mean_temp_schimmelrader, "mean_temp_schimmelrader.nc", format = "CDF",overwrite =T)
spat_mean_temp_schimmelrader <- rast(mean_temp_schimmelrader)


spat_mean_temp_schimmelrader_250m <- terra::disagg(spat_mean_temp_schimmelrader, 40)
spat_mean_wind_schimmelrader_250m
plot(spat_mean_temp_schimmelrader_250m)
spat_mean_temp_schimmelrader_250m <- project(spat_mean_temp_schimmelrader_250m, DistanceSeaRastNA)
writeRaster(spat_mean_temp_schimmelrader_250m, "spat_mean_temp_schimmelrader_250m.tif", overwrite = TRUE)

spat_mean_temp_schimmelrader_250m <- rast("spat_mean_temp_schimmelrader_250m.tif")

Schimmelradar360spat_focal$temp <- terra::extract(spat_mean_temp_schimmelrader_250m, Schimmelradar360spat_focal)[,2]

mock_raster_temp_NL <- mask(spat_mean_temp_schimmelrader_250m, vect(provincesNL_transformed), filename="Mock_temp_raster_NL.tif", overwrite = T)
plot(mock_raster_temp_NL)
```

# Humidity data extraction and integration
```{r}
#Extracting humidity data from AGERA5
humid.nc <- list.files("Weather_schimmelradar/Relative_humidity", full.names = T)
humid_stack <- raster::stack(humid.nc)
mean_humid_schimmelrader <- mean(humid_stack)
str(mean_humid_schimmelrader)
mean_humid_schimmelrader@data@values[c(104, 143, 144, 145, 146)]<- mean_humid_schimmelrader@data@values[c(105, 185, 186, 187, mean(147, 188))] # plugging the gap of Vlieland missing values by copying values from neighboring cells, in the case of the one empty cell with two neighboring cells we take the mean of the two
writeRaster(mean_humid_schimmelrader, "mean_humid_schimmelrader.nc", format = "CDF",overwrite =T)
spat_mean_humid_schimmelrader <- rast(mean_humid_schimmelrader)


spat_mean_humid_schimmelrader_250m <- terra::disagg(spat_mean_humid_schimmelrader, 40)
spat_mean_humid_schimmelrader_250m
plot(spat_mean_humid_schimmelrader_250m)
spat_mean_humid_schimmelrader_250m <- project(spat_mean_humid_schimmelrader_250m, DistanceSeaRastNA)
writeRaster(spat_mean_humid_schimmelrader_250m, "spat_mean_humid_schimmelrader_250m.tif", overwrite = TRUE)

spat_mean_humid_schimmelrader_250m <- rast("spat_mean_humid_schimmelrader_250m.tif")

Schimmelradar360spat_focal$relative_humidity <- terra::extract(spat_mean_humid_schimmelrader_250m, Schimmelradar360spat_focal)[,2]

mock_raster_humid_NL <- mask(spat_mean_humid_schimmelrader_250m, vect(provincesNL_transformed), filename="Mock_humid_raster_NL.tif", overwrite = T)
plot(mock_raster_humid_NL)
```

# Precipitation data extraction and integration
```{r}
#Extracting precipitation data from AGERA5
precip.nc <- list.files("Weather_schimmelradar/Precipitation", full.names = T)
precip_stack <- raster::stack(precip.nc)
mean_precip_schimmelrader <- mean(precip_stack)
str(mean_precip_schimmelrader)
mean_precip_schimmelrader@data@values[c(104, 143, 144, 145, 146)]<- mean_precip_schimmelrader@data@values[c(105, 185, 186, 187, mean(147, 188))] # plugging the gap of Vlieland missing values by copying values from neighboring cells, in the case of the one empty cell with two neighboring cells we take the mean of the two
writeRaster(mean_precip_schimmelrader, "mean_precip_schimmelrader.nc", format = "CDF",overwrite =T)
spat_mean_precip_schimmelrader <- rast(mean_precip_schimmelrader)


spat_mean_precip_schimmelrader_250m <- terra::disagg(spat_mean_precip_schimmelrader, 40)
spat_mean_precip_schimmelrader_250m
plot(spat_mean_precip_schimmelrader_250m)
spat_mean_precip_schimmelrader_250m <- project(spat_mean_precip_schimmelrader_250m, DistanceSeaRastNA)
writeRaster(spat_mean_precip_schimmelrader_250m, "spat_mean_precip_schimmelrader_250m.tif", overwrite = TRUE)

spat_mean_precip_schimmelrader_250m <- rast("spat_mean_precip_schimmelrader_250m.tif")

Schimmelradar360spat_focal$precipitation <- terra::extract(spat_mean_precip_schimmelrader_250m, Schimmelradar360spat_focal)[,2]

mock_raster_precip_NL <- mask(spat_mean_precip_schimmelrader_250m, vect(provincesNL_transformed), filename="Mock_precip_raster_NL.tif", overwrite = T)
plot(mock_raster_precip_NL)
```

#Make raster with predictor variables covering the Netherlands
```{r}
Predictor_raster <- c(grass22fra, buildings22fra, maize22fra, bulb22fra, tulip22fra, lily22fra, daffodil22fra, gladiolus22fra, hyacinth22fra, daffodil_hyacinth22fra, greenhouse22fra, potato22fra, beetroot_potato22fra, seed_potato22fra, starch_potato22fra, consumption_potato22fra, beetroot22fra, forest22fra, Onion22fra, seed_onion22fra, set_onion22fra, Othercrops22fra, water22fra, fruit22fra, apple22fra, pear22fra, apple_pear22fra, tree_farming22fra, heather22fra, Peatland22fra, grain22fra, grass22fra5k, buildings22fra5k, maize22fra5k, bulb22fra5k, tulip22fra5k, lily22fra5k, daffodil22fra5k, gladiolus22fra5k, hyacinth22fra5k, daffodil_hyacinth22fra5k, greenhouse22fra5k, potato22fra5k, beetroot_potato22fra5k, seed_potato22fra5k, starch_potato22fra5k, consumption_potato22fra5k, beetroot22fra5k, forest22fra5k, Onion22fra5k, seed_onion22fra5k, set_onion22fra5k, Othercrops22fra5k, water22fra5k, fruit22fra5k, apple22fra5k, pear22fra5k, apple_pear22fra5k, tree_farming22fra5k, heather22fra5k, Peatland22fra5k, grain22fra5k, spat_mean_wind_schimmelrader_250m, spat_mean_temp_schimmelrader_250m, spat_mean_humid_schimmelrader_250m, spat_mean_precip_schimmelrader_250m, DistanceSeaRastNA)


names(Predictor_raster) <- c("grass22fra", "buildings22fra", "maize22fra", "bulb22fra", "tulip22fra", "lily22fra", "daffodil22fra", "gladiolus22fra", "hyacinth22fra", "daffodil_hyacinth22fra", "greenhouse22fra", "potato22fra", "beetroot_potato22fra", "seed_potato22fra", "starch_potato22fra", "consumption_potato22fra", "beetroot22fra", "forest22fra", "Onion22fra", "seed_onion22fra", "set_onion22fra", "Othercrops22fra", "water22fra", "fruit22fra", "apple22fra", "pear22fra", "apple_pear22fra","tree_farming22fra", "heather22fra", "Peatland22fra", "grain22fra", "grass22fra5k", "buildings22fra5k", "maize22fra5k", "bulb22fra5k", "tulip22fra5k", "lily22fra5k", "daffodil22fra5k", "gladiolus22fra5k", "hyacinth22fra5k", "daffodil_hyacinth22fra5k", "greenhouse22fra5k", "potato22fra5k", "beetroot_potato22fra5k", "seed_potato22fra5k", "starch_potato22fra5k", "consumption_potato22fra5k", "beetroot22fra5k", "forest22fra5k", "Onion22fra5k", "seed_onion22fra5k", "set_onion22fra5k", "Othercrops22fra5k", "water22fra5k", "fruit22fra5k", "apple22fra5k", "pear22fra5k", "apple_pear22fra5k", "tree_farming22fra5k", "heather22fra5k", "Peatland22fra5k", "grain22fra5k", "wind_speed", "temp", "relative_humidity", "precipitation", "distance_sea")

Predictor_raster_NL <- mask(Predictor_raster, vect(provincesNL_transformed), filename="Predictor_raster_NL.tif", overwrite = T)
Predictor_raster_NL <- rast("Predictor_raster_NL.tif")

plot(Predictor_raster_NL)
```


# Extract fraction of landuse within the disks of the variables of interst (10k radius)
```{r}
# WARNING: these take ~2 hours to run altogether load "Schimmelradar360spat_round.RData"
Schimmelradar360spat_focal$distance_sea <- extract(DistanceSeaRastNA, Schimmelradar360spat_focal, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$forest22fra <- extract(forest22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$bulb22fra <- extract(bulb22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$tulip22fra <- extract(tulip22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$lily22fra <- extract(lily22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$daffodil22fra <- extract(daffodil22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$gladiolus22fra <- extract(gladiolus22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$hyacinth22fra <- extract(hyacinth22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$daffodil_hyacinth22fra <- extract(daffodil_hyacinth22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Onion22fra <- extract(Onion22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$seed_onion22fra <- extract(seed_onion22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$set_onion22fra <- extract(set_onion22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Othercrops22fra <- extract(Othercrops22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$maize22fra <- extract(maize22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$apple_pear22fra <- extract(apple_pear22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$apple22fra <- extract(apple22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$pear22fra <- extract(pear22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$fruit22fra <- extract(fruit22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$tree_farming22fra <- extract(tree_farming22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$heather22fra <- extract(heather22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Peatland22fra <- extract(Peatland22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$buildings22fra <- extract(buildings22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$grass22fra <- extract(grass22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$potato22fra <- extract(potato22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$beetroot_potato22fra <- extract(beetroot_potato22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$seed_potato22fra <- extract(seed_potato22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$starch_potato22fra <- extract(starch_potato22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$consumption_potato22fra <- extract(consumption_potato22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$beetroot22fra <- extract(beetroot22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$grain22fra <- extract(grain22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$greenhouse22fra <- extract(greenhouse22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$water22fra <- extract(water22fra, datbuf_focal, fun=mean, na.rm=TRUE)[,2]
#save(Schimmelradar360spat_focal,file =  "Schimmelradar360spat_focal.RData")
load("Schimmelradar360spat_focal.RData")
```

# Extract fraction of land within the disks if the variables of interst (5k radius)
```{r}
# WARNING: these take ~2 hours to run altogether load "Schimmelradar360spat_round.RData"
Schimmelradar360spat_focal$forest22fra5k <- extract(forest22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$bulb22fra5k <- extract(bulb22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$tulip22fra5k <- extract(tulip22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$lily22fra5k <- extract(lily22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$daffodil22fra5k <- extract(daffodil22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$gladiolus22fra5k <- extract(gladiolus22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$hyacinth22fra5k <- extract(hyacinth22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$daffodil_hyacinth22fra5k <- extract(daffodil_hyacinth22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Onion22fra5k <- extract(Onion22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$seed_onion22fra5k <- extract(seed_onion22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$set_onion22fra5k <- extract(set_onion22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Othercrops22fra5k <- extract(Othercrops22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$maize22fra5k <- extract(maize22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$apple_pear22fra5k <- extract(apple_pear22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$apple22fra5k <- extract(apple22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$pear22fra5k <- extract(pear22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$fruit22fra5k <- extract(fruit22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$tree_farming22fra5k <- extract(tree_farming22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$heather22fra5k <- extract(heather22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$Peatland22fra5k <- extract(Peatland22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$buildings22fra5k <- extract(buildings22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$grass22fra5k <- extract(grass22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$potato22fra5k <- extract(potato22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$beetroot_potato22fra5k <- extract(beetroot_potato22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$seed_potato22fra5k <- extract(seed_potato22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$starch_potato22fra5k <- extract(starch_potato22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$consumption_potato22fra5k <- extract(consumption_potato22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$beetroot22fra5k <- extract(beetroot22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$grain22fra5k <- extract(grain22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$greenhouse22fra5k <- extract(greenhouse22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
Schimmelradar360spat_focal$water22fra5k <- extract(water22fra5k, datbuf_focal5k, fun=mean, na.rm=TRUE)[,2]
#save(Schimmelradar360spat_focal,file =  "Schimmelradar360spat_focal.RData")
load("Schimmelradar360spat_focal.RData")
```

#Plotting predictor raster for all 10km radius land use data extractions + weather data (Figures S20 and S21)
```{r}
plot(Predictor_raster_NL[[2:9]])
par(bg = "black", fg = "white", col.axis = "white", col.lab = "white", col.main = "white")
plot(Predictor_raster_NL[[6]], colNA = "black", main = "")
plot(Predictor_raster_NL[[c(17:21,23:28,31,64:67)]])
plot(Predictor_raster_NL[[64:67]])
Urbanistation_heatmap <- plot(Predictor_raster_NL[[2]], axes = FALSE, main = "Urbanization")

# Export the plot to a PDF
pdf("Urbanization_heat.pdf", width = 2.33, height = 1.66)  # Specify dimensions (in inches)
plot(Predictor_raster_NL[[2]], axes = FALSE, main = "", legend = FALSE) # Plot without axis labels and with a title
dev.off()  # Close the PDF device

# Export the plot to a PDF
pdf("Maize_heat.pdf", width = 2.33, height = 1.66)  # Specify dimensions (in inches)
plot(Predictor_raster_NL[[3]], axes = FALSE, main = "", legend = FALSE) # Plot without axis labels and with a title
dev.off()  # Close the PDF device

# Export the plot to a PDF
pdf("Bulb_heat.pdf", width = 2.33, height = 1.66)  # Specify dimensions (in inches)
plot(Predictor_raster_NL[[4]], axes = FALSE, main = "", legend = FALSE) # Plot without axis labels and with a title
dev.off()  # Close the PDF device
```

#Adding in the TR haplotype data to the phenotypic resistance data and generate table S2
```{r}
#load data
Schimmelradar_TR <- read.csv("Analysis_genotyping_schimmelradar_TR_types.csv")
# turn short format

# Replace NA with 0
Schimmelradar_TR[is.na(Schimmelradar_TR)] <- 0

# Group by Sample and sum the values (Table S2)
Schimmelradar_TR_short <- Schimmelradar_TR %>%
  group_by(Sample, azole) %>%
  summarise(across(wt:TR46_4, sum, na.rm = TRUE))
#write.csv(Schimmelradar_TR_short, "Schimmelradar_TR_per_sample.csv")

Schimmelradar_TR_VORI <- Schimmelradar_TR_short[Schimmelradar_TR_short$azole == "VORI", ]
Schimmelradar_TR_ITRA <- Schimmelradar_TR_short[Schimmelradar_TR_short$azole == "ITRA", ]


# Group by Azole and sum the values for summary table
Schimmelradar_TR_azole <- Schimmelradar_TR %>%
  group_by(azole) %>%
  summarise(across(wt:TR46_4, sum, na.rm = TRUE))
#write.csv(Schimmelradar_TR_azole, "Schimmelradar_TR_per_azole.csv")

Schimmelradar_TR_azole.df <- as.data.frame(Schimmelradar_TR_azole)



#Merge into spatial object but only samples already present in the phenotypic dataset
Schimmelradar360spat_focal_TR <- Schimmelradar360spat_focal %>%
  left_join(Schimmelradar_TR_short[,-2], by = "Sample")

#Merge into spatial object but only samples already present in the phenotypic dataset
Schimmelradar360spat_focal_TR_VORI <- Schimmelradar360spat_focal %>%
  left_join(Schimmelradar_TR_VORI[,-2], by = "Sample")

#Merge into spatial object but only samples already present in the phenotypic dataset
Schimmelradar360spat_focal_TR_ITRA <- Schimmelradar360spat_focal %>%
  left_join(Schimmelradar_TR_ITRA[,-2], by = "Sample")


# Identify the columns added during the merge (i.e., those in sample_df but not in sample_sf)
new_columns <- setdiff(names(Schimmelradar_TR_short[,-2]), names(Schimmelradar360spat_focal))

# Identify the columns added during the merge (i.e., those in sample_df but not in sample_sf)
new_columns_VORI <- setdiff(names(Schimmelradar_TR_VORI[,-2]), names(Schimmelradar360spat_focal))

# Identify the columns added during the merge (i.e., those in sample_df but not in sample_sf)
new_columns_ITRA <- setdiff(names(Schimmelradar_TR_ITRA[,-2]), names(Schimmelradar360spat_focal))

#Replace NA values with 0, as we have 0 
Schimmelradar360spat_focal_TR <- Schimmelradar360spat_focal_TR %>%
  mutate(across(all_of(new_columns), ~replace_na(., 0)))

#Replace NA values with 0, as we have 0 
Schimmelradar360spat_focal_TR_VORI <- Schimmelradar360spat_focal_TR_VORI %>%
  mutate(across(all_of(new_columns_VORI), ~replace_na(., 0)))

#Replace NA values with 0, as we have 0 
Schimmelradar360spat_focal_TR_ITRA <- Schimmelradar360spat_focal_TR_ITRA %>%
  mutate(across(all_of(new_columns_ITRA), ~replace_na(., 0)))

Schimmelradar360spat_focal_TR$Triazole_tot <- Schimmelradar360spat_focal_TR$Tot.VOR + Schimmelradar360spat_focal_TR$Tot.ITR

Schimmelradar360spat_focal_TR$WT_Rfrac <- Schimmelradar360spat_focal_TR$wt / Schimmelradar360spat_focal_TR$Triazole_tot


Schimmelradar360spat_focal_TR$TR34frac <- Schimmelradar360spat_focal_TR$TR34 / Schimmelradar360spat_focal_TR$Triazole_tot


#Itraconazole resistant TR34 fraction
Schimmelradar360spat_focal_TR_ITRA$TR34frac <- Schimmelradar360spat_focal_TR_ITRA$TR34 / Schimmelradar360spat_focal_TR_ITRA$Tot.ITR
#Itraconazole resistant TR46 fraction
Schimmelradar360spat_focal_TR_ITRA$TR46ITRfrac <- Schimmelradar360spat_focal_TR_ITRA$TR46 / Schimmelradar360spat_focal_TR_ITRA$Tot.ITR
#Itraconazole resistant WT fraction
Schimmelradar360spat_focal_TR_ITRA$WTITRfrac <- Schimmelradar360spat_focal_TR_ITRA$wt / Schimmelradar360spat_focal_TR_ITRA$Tot.ITR

#Voriconazole resistant TR46 fraction
Schimmelradar360spat_focal_TR_VORI$TR46frac <- Schimmelradar360spat_focal_TR_VORI$TR46 / Schimmelradar360spat_focal_TR_VORI$Tot.VOR
#Voriconazole resistant TR34 fraction
Schimmelradar360spat_focal_TR_VORI$TR34VORIfrac <- Schimmelradar360spat_focal_TR_VORI$TR34 / Schimmelradar360spat_focal_TR_VORI$Tot.VOR
save(Schimmelradar360spat_focal_TR_VORI, file = "Schimmelradar360spat_focal_TR_VORI.RData")
save(Schimmelradar360spat_focal_TR_ITRA, file = "Schimmelradar360spat_focal_TR_ITRA.RData")
```

#Statistical test for differing genotype frequencies TR34 vs TR46 per triazole treatment
```{r}
#Test if the ratio of TR34 to TR46 differs significantly between 
rownames(Schimmelradar_TR_azole.df) <- Schimmelradar_TR_azole.df[,1]
Schimmelradar_TR_azole.df <- Schimmelradar_TR_azole.df[,-1]
as.matrix(Schimmelradar_TR_azole.df[,2:3])
mcnemar.test(as.matrix(Schimmelradar_TR_azole.df[,2:3]))
```

