update RPG calculation

Author: rfrelat

README.md

@@ -45,6 +45,12 @@ quarto::quarto_render("analysis/06_haie.qmd")
 quarto::quarto_render("analysis/07_fields_Toulouse.qmd")
 ```
 
+#### 8. exploration of the Swiss dataset
+
+```r
+quarto::quarto_render("analysis/08_Swiss_data.qmd")
+```
+
 
 #### Mise à jour de l'index:
 ```r

analysis/05_summary.qmd

@@ -67,6 +67,8 @@ indr <- read.csv(here("data", "derived-data", "raster_indicators.csv"))
 
 # rpg and oso classes
 ref <- read.csv(here("data", "derived-data", "RPG-OSO_classes.csv"))
+ref_rpg <- read.csv(here("data", "raw-data", "RPG", "RPG_Explorer_Annexe2.csv"))
+
 
 # vector indicators
 indv <- read.csv(here("data", "derived-data", "vector_indicators.csv"))

analysis/06_haie.qmd

@@ -189,7 +189,14 @@ mapview(
 ```
 
 
+### IGN
 
+```{r}
+library(happign)
+meta <- get_layers_metadata("wfs")
+meta <- get_layers_metadata("wmr")
+meta$Name[grep("haie", meta$Name)]
+```
 
 
 ## Comparaison des indicateurs

analysis/F01_get_plot_CH.R

@@ -0,0 +1,147 @@
+# Get all plots around the fields coordinates in a buffer
+# considering the year of the measurement
+
+# For the Swiss dataset, we use data received from Selma Cadot
+# Run in ~1min
+
+library(terra)
+library(sf)
+library(here)
+
+# Load home made functions
+devtools::load_all()
+
+# The data was split per year using S01_prep_data.R
+# data are available for the period 2019 - 2023
+period <- 2019:2023
+nutz_layer <- "nutz_XXXX.gpkg"
+labRPG <- c("id_parcel", "code_cultu", "code_group")
+colNUTZ <- c("nutzungsidentifikator", "nutzung_fr", "Hauptkategorie_FR")
+
+buffer_fields <- c(500, 1000, 1500) #in m
+years <- 5 # for crop rotation
+
+datafolder <- here("data", "raw-data")
+outfolder <- here("data", "derived-data")
+
+# get the coordinates from the points
+# from shinyFunbiodiv/analysis/03_update_data.R
+df <- read.csv(here(datafolder, "coordinates_year_crop.csv"))
+pts <- vect(df, geom = c("Long", "Lat"), crs = "EPSG:4326")
+# project in EPSG 2056
+pts <- project(pts, "EPSG:2056")
+
+# select only observation in Switzerland
+keep <- pts$Study_ID %in% "PestiRed" & !is.na(df$Lat)
+# table(keep) # 203 points
+# table(pts$Year[keep])
+
+# table(nutz$Hauptkategorie_FR, useNA="ifany")
+rmCat <- c(
+  "Forêt",
+  "Haies, bosquets et berges boisées",
+  "Surfaces en dehors de la SAU"
+)
+
+nutz_out <- c() # save the nutz information
+df_out <- c() # save the parameters
+# for testing: sample(which(keep), 10)
+for (i in which(keep)) {
+  cat(i)
+
+  # select point i
+  pti <- pts[i]
+
+  # create maximum buffer
+  buffi <- buffer(pti, max(buffer_fields))
+
+  # get rpg around buffer
+  nutzi <- vect(
+    file.path(datafolder, "CH", gsub("XXXX", pti$Year, nutz_layer)),
+    ext = buffi
+  )
+
+  # remove non agricultural fields
+  nutzi <- nutzi[!nutzi$Hauptkategorie_FR %in% rmCat, ]
+  # continue only if some fields in nutz
+  if (nrow(nutzi) > 0) {
+    # calculate geometrical characteristics
+    nutzi$Perim_m <- perim(nutzi)
+    nutzi$Area_ha <- expanse(nutzi) * 0.0001
+
+    ##
+    # ponctual information on the site
+    cat(".")
+    rei <- terra::relate(nutzi, pti, "intersects")
+
+    if (sum(rei) > 0) {
+      exi <- nutzi[which(rei)[1], ]
+      # add funbiodiv ID
+      exi$Funbiodiv_ID <- pti$ID
+      # keep in nutz_out
+      nutz_out <- c(nutz_out, exi)
+      out_i <- data.frame(
+        data.frame(pti),
+        data.frame(exi)[, c(colNUTZ, "Perim_m", "Area_ha")]
+      )
+    } else {
+      out_i <- data.frame(
+        data.frame(pti),
+        matrix(NA, ncol = length(colNUTZ) + 2)
+      )
+      #fmt: skip
+      names(out_i)[(ncol(pti)+1):ncol(out_i)] <- c(colNUTZ, "Perim_m", "Area_ha")
+    }
+
+    ##
+    # average field size in buffer
+    cat(".")
+    for (f in buffer_fields) {
+      lab <- paste0("Mean_fieldsize_", f, "m_ha")
+      buf_pts <- buffer(pti, f)
+      nutz_buf <- relate(nutzi, buf_pts, "intersects")
+      out_i[, lab] <- mean(nutzi$Area_ha[nutz_buf], na.rm = TRUE)
+    }
+
+    ##
+    # crop rotation
+    cat(".")
+    if (sum(rei) > 0) {
+      nutz_time <- exi[, colNUTZ]
+      nutz_time$Year <- pti$Year
+    } else {
+      nutz_time <- vect()
+    }
+    for (y in 1:years) {
+      out_i[paste0("N-", y)] <- pti$Year - y
+      lab <- paste0(colNUTZ, "_N-", y)
+      out_i[, lab] <- NA
+      if ((pti$Year - y) %in% period) {
+        nutzy <- vect(
+          file.path(datafolder, "CH", gsub("XXXX", pti$Year - y, nutz_layer)),
+          ext = pti
+        )
+
+        if (nrow(nutzy) > 0) {
+          rey <- relate(nutzy, pti, "intersects")
+          if (sum(rey) > 1) {
+            nutzy <- nutzi[which(rey)[1], ]
+            # add information in out_i
+            out_i[, lab] <- data.frame(nutzy)[, colNUTZ]
+            # keep it for save
+            nutzy$Year <- pti$Year - y
+            nutz_time <- rbind(nutz_time, nutzy[c(colNUTZ, "Year")])
+          }
+        }
+      }
+    }
+    df_out <- rbind(df_out, out_i)
+  }
+}
+
+# export indicators
+write.csv(df_out, file.path(outfolder, "metrics_nutz.csv"), row.names = FALSE)
+
+# save rpg fields
+nutz_all <- do.call(rbind, nutz_out)
+writeVector(nutz_all, file.path(outfolder, "nutz_fields.gpkg"))

analysis/F01_get_plot_FR.R

@@ -0,0 +1,161 @@
+# Get all plots around the fields coordinates in a buffer
+# considering the year of the measurement
+
+# For the French dataset, we will download the interesting part from the API directly
+# using happign package: https://paul-carteron.github.io/happign
+library(happign)
+
+# Run in ~2h
+
+library(terra)
+library(sf)
+library(here)
+
+# Load home made functions
+devtools::load_all()
+
+# check out which RPG layer is available
+# meta_vect <- get_layers_metadata("wfs") # all layers for altimetrie wms
+# look_up <- "RPG"
+# #fmt: skip
+# found <- grepl(tolower(look_up), tolower(meta_vect$Name)) | grepl(tolower(look_up), tolower(meta_vect$Name))
+# meta_vect$Name[found]
+period <- 2015:2024
+rpg_layer <- "RPG.XXXX:parcelles_graphiques"
+colRPG <- c("id_parcel", "code_cultu", "code_group")
+
+buffer_fields <- c(500, 1000, 1500) #in m
+years <- 5 # for crop rotation
+
+datafolder <- here("data", "raw-data")
+outfolder <- here("data", "derived-data")
+
+# remove jachere? J5M, J6P, J6S, JAC, JNO
+rmCat <- c("BFP", "BFS", "BOR", "BTA", "SNA")
+
+# get the coordinates from the points
+# from shinyFunbiodiv/analysis/03_update_data.R
+df <- read.csv(here(datafolder, "coordinates_year_crop.csv"))
+pts <- vect(df, geom = c("Long", "Lat"), crs = "EPSG:4326")
+# table(pts$Year)
+
+# select only observation in France between 2015 and 2024
+keep <- pts$Year %in% period & !pts$Study_ID %in% "PestiRed"
+# table(keep) # 1692 points
+
+rpg_out <- c() # save the rpg information
+df_out <- c() # save the parameters
+# for testing: sample(which(keep), 10)
+for (i in which(keep)) {
+  cat(i)
+
+  # select point i
+  pti <- pts[i]
+
+  # create maximum buffer
+  buffi <- buffer(pti, max(buffer_fields))
+
+  # get rpg around buffer
+  rpgi <- get_wfs(
+    x = st_as_sf(buffi),
+    layer = gsub("XXXX", pti$Year, rpg_layer)
+  )
+  # remove non agricultural fields
+  rpgi <- rpgi[!rpgi$code_cultu %in% rmCat, ]
+
+  # continue only if some fields in rpg
+  if (nrow(rpgi) > 0) {
+    # transform as terra SpatVect object
+    rpgi <- vect(rpgi)
+
+    # save (in case api doesn't work anymore)
+    outi <- file.path(datafolder, "happign", paste0("RPGbuf_", pti$ID, ".gpkg"))
+    writeVector(rpgi, outi, overwrite = TRUE)
+
+    # calculate geometrical characteristics
+    rpgi$Perim_m <- perim(rpgi)
+    rpgi$Area_ha <- expanse(rpgi) * 0.0001
+
+    ##
+    # ponctual information on the site
+    cat(".")
+    rei <- relate(rpgi, pti, "intersects")
+    if (sum(rei) == 1) {
+      exi <- rpgi[rei]
+      # keep in rpg_out
+      rpg_out <- c(rpg_out, exi)
+      out_i <- data.frame(
+        data.frame(pti),
+        data.frame(exi)[, c(colRPG, "Perim_m", "Area_ha")]
+      )
+    } else {
+      out_i <- data.frame(
+        data.frame(pti),
+        matrix(NA, ncol = length(colRPG) + 2)
+      )
+      #fmt: skip
+      names(out_i)[(ncol(pti)+1):ncol(out_i)] <- c(colRPG, "Perim_m", "Area_ha")
+    }
+
+    ##
+    # average field size in buffer
+    cat(".")
+    for (f in buffer_fields) {
+      lab <- paste0("Mean_fieldsize_", f, "m_ha")
+      buf_pts <- buffer(pti, f)
+      rpg_buf <- relate(rpgi, buf_pts, "intersects")
+      out_i[, lab] <- mean(rpgi$Area_ha[rpg_buf], na.rm = TRUE)
+    }
+
+    ##
+    # crop rotation
+    cat(".")
+    if (sum(rei) == 1) {
+      rpg_time <- exi[, colRPG]
+      rpg_time$Year <- pti$Year
+    } else {
+      rpg_time <- vect()
+    }
+    for (y in 1:years) {
+      out_i[paste0("N-", y)] <- pti$Year - y
+      lab <- paste0(colRPG, "_N-", y)
+      out_i[, lab] <- NA
+      if ((pti$Year - y) %in% period) {
+        rpgy <- get_wfs(
+          x = st_as_sf(pti),
+          layer = gsub("XXXX", pti$Year - y, rpg_layer)
+        )
+
+        if (nrow(rpgy) > 0) {
+          rpgy <- vect(rpgy)
+          rey <- relate(rpgy, pti, "intersects")
+          if (sum(rey) == 1) {
+            rpgy <- rpgy[rey]
+            # add information in out_i
+            out_i[, lab] <- data.frame(rpgy)[, colRPG]
+            # keep it for save
+            rpgy$Year <- pti$Year - y
+            rpg_time <- rbind(rpg_time, rpgy[c(colRPG, "Year")])
+          }
+        }
+      }
+    }
+    if (nrow(rpg_time) > 1) {
+      outi <- file.path(
+        datafolder,
+        "happign",
+        paste0("RPGrot_", pti$ID, ".gpkg")
+      )
+      writeVector(rpg_time, outi, overwrite = TRUE)
+    }
+
+    df_out <- rbind(df_out, out_i)
+  }
+}
+
+# export indicators
+write.csv(df_out, file.path(outfolder, "metrics_rpg.csv"), row.names = FALSE)
+
+# save rpg fields
+rpg_all <- do.call(rbind, rpg_out)
+writeVector(rpg_all, file.path(outfolder, "rpg_fields.gpkg"))

analysis/P01_get_points.R

@@ -98,9 +98,9 @@ st_write(
 )
 
 ## Get classes of RPG and OSO
-rpg <- read.csv2("data/raw-data/RPG/REF_CULTURES_GROUPES_CULTURES_2023.csv", )
+# rpg <- read.csv2("data/raw-data/RPG/REF_CULTURES_GROUPES_CULTURES_2023.csv", )
 
-rpg2 <- read.csv2("data/raw-data/RPG/REF_CULTURES_2023.csv")
-keep <- which(rpg2$CAMPAGNE_DEBUT > 2015 | rpg2$CAMPAGNE_FIN > 2015)
-rpg2$LIBELLE_CULTURE[keep][!rpg2$CODE[keep] %in% rpg$CODE_CULTURE]
-rpg$CODE_CULTURE %in% rpg2$CODE[keep]
+# rpg2 <- read.csv2("data/raw-data/RPG/REF_CULTURES_2023.csv")
+# keep <- which(rpg2$CAMPAGNE_DEBUT > 2015 | rpg2$CAMPAGNE_FIN > 2015)
+# rpg2$LIBELLE_CULTURE[keep][!rpg2$CODE[keep] %in% rpg$CODE_CULTURE]
+# rpg$CODE_CULTURE %in% rpg2$CODE[keep]

analysis/S01_prep_data.R

@@ -0,0 +1,20 @@
+library(terra)
+library(here)
+
+dirfile <- here("data", "raw-data", "CH")
+file <- file.path(dirfile, "buffer_4000_lnf_hk_area_distance_nutz_60.gpkg")
+
+# super heavy (2.4Gb)
+rpg_ch <- vect(file)
+
+# split the file per year (400Mb)
+for (i in 2019:2023) {
+  chi <- rpg_ch[rpg_ch$year == i, ]
+  # could also select the relevant column to make the data smaller
+  writeVector(chi, file.path(dirfile, paste0("nutz_", i, ".gpkg")))
+}
+
+zoomin <- ext(c(2508265, 2509047, 1119836, 1121222)) + 2000
+v <- vect(file.path(dirfile, "nutz_2019.gpkg"), extent = zoomin)
+
+mapview::mapview(v)

make.R

@@ -0,0 +1,7 @@
+system.time({
+  source("analysis/F01_get_plot_FR.R")
+})
+
+system.time({
+  source("analysis/F01_get_plot_CH.R")
+})