R/TarifIFNSER.R

In pobsteta/gftools: Geographic Data Analysis and Modeling for Forests

#' Tarif de cubage volume geometrique bois fort
#'
#' @description La fonction renvoie un tableau des types et numeros de tarif Schaeffer par essence
#' en utilisant la base de donnees arbres de l'IFN. Il s'agit d'un volume geometrique bois fort tige.
#' La fonction necessite en entree un shape correspondant au perimetre retenu : forêt, massif,
#' sylvoecoregion, etc.
#'
#' @return La fonction renvoie un tableau des types et numeros de tarif Schaeffer par essence.
#'
#' @param seuilCircf = 50 cm par defaut, seuil minimal de circonference en cm en dessous duquel l'arbre n'est pas retenu
#' pour le calcul du tarif.
#' @param seuilNb = 10 individus par defaut, seuil minimal du nombre d'arbres dans la base IFN pour qu'une essence soit retenue.
#' @param resolution = 500 m par defaut, resolution de la grille servant de support au krigeage
#' @param useSer = FALSE par defaut, argument permettant de choisir si le calcul des numeros de tarif doit se faire
#' au sein d'une même sylvoecoregion
#' @param distmax = 20 000 m par defaut, distance maximale a ne pas depasser pour la recherche de placette IFN
#' @param enreg = FALSE par defaut, argument permettant de choisir l'enregistrement ou pas du tableau resultat au format csv
#' @param fichier = nom du fichier (chemin du repertoire compris) de la zone sur laquelle sont calcules les resultats
#' @param annee = c(2005:2016) par defaut, vecteur indiquant les donnees IGN servant au calcul des resultats
#'
#' @author Bruciamacchie Max, Pascal Obstetar
#'
#' @export TarifIFNSER
#'
#' @examples
#' \dontrun{
#' TarifIFNSER(fichier = '~/ec21.shp', enreg = F)
#' }
#'

TarifIFNSER <- function(fichier = NULL, annee = c(2005:2016), seuilCircf = 50, seuilNb = 10, resolution = 500, useSer = F, distmax = 20000, enreg = F) {
  if (is.null(fichier)) {
    stop("Your fichier is NULL!")
  }
  if (file.access(fichier) == -1) {
    stop(paste("You have no access to file:", fichier))
  }
  vecteur_annee <- annee
  donnees <- system.file("extdata/IFN", package = "gftools")
  ser <- sf::st_read(system.file("extdata/shapes/SER/ser_l93.shp", package = "gftools"), quiet = TRUE)
  perimetre <- sf::st_read(fichier, quiet = TRUE)
  perimetre <- sf::st_transform(perimetre, crs = 2154)
  ser <- sf::st_transform(ser, crs = 2154)
  perimetre <- sf::st_union(perimetre)
  surf <- as.numeric(sf::st_area(perimetre))
  zone <- sf::st_buffer(perimetre, dist = distmax)
  CodesEssIFN <- gftools::getData("IFNCODE") %>% dplyr::filter(donnee %in% "ESPAR")
  if (useSer) {
    zone <- sf::st_intersection(zone, ser)
    zone <- sf::st_union(zone)
  }
  # Extraction placettes et arbres
  message("Extract data placettes...")
  plac <- data.frame()
  for (i in 1:length(vecteur_annee)) {
    yrs <- as.numeric(vecteur_annee[i])
    ifn <- gftools::getFich_IFN(obj = c("Pla"), Peup = FALSE, Morts = FALSE, Doc = TRUE, ans = yrs, doss = donnees)
    placet <- data.frame(yrs = yrs, ifn$Pla[[1]][, c("idp", "ser", "xl93", "yl93")])
    plac <- rbind(plac, placet)
  }
  message("Extract data arbres...")
  arb <- data.frame()
  for (i in 1:length(vecteur_annee)) {
    yrs <- as.numeric(vecteur_annee[i])
    ifn <- gftools::getFich_IFN(obj = c("Arb"), Peup = FALSE, Morts = FALSE, Doc = TRUE, ans = yrs, doss = donnees)
    arbre <- data.frame(yrs = yrs, ifn$Arb[[1]][, c("idp", "espar", "c13", "w", "v")])
    arb <- rbind(arb, arbre)
  }
  placettes <- sf::st_intersection(sf::st_as_sf(plac, coords = c("xl93", "yl93"), crs = 2154, agr = "constant"), zone)
  if (length(placettes$idp) == 0) {
    stop(paste("Increase 'distmax' =", distmax, "to have placettes in your zone!"))
  }
  Tarifs <- arb %>%
    dplyr::filter(idp %in% placettes$idp) %>%
    dplyr::mutate(Diam = c13 / pi, numSchR = v / 5 * 70000 / (Diam - 5) / (Diam - 10) - 8, numSchL = v / 5 * 90000 / Diam / (Diam - 5) - 8) %>%
    dplyr::left_join(placettes, by = "idp") %>%
    dplyr::select(espar, c13, numSchR, numSchL, geometry) %>%
    na.omit()
  message("Extraction of species according to the threshold...")
  t2 <- Tarifs %>%
    dplyr::filter(c13 >= seuilCircf)
  df <- data.frame(table(t2$espar))
  names(df) <- c("code", "Freq")
  df <- df %>%
    dplyr::filter(df$Freq >= seuilNb) %>%
    dplyr::left_join(CodesEssIFN[, 2:3], by = c("code", "code")) %>%
    dplyr::arrange(desc(Freq))
  # Creation du grid
  Tab <- data.frame()
  BB <- sf::st_bbox(zone)
  BB <- resolution * (round(BB / resolution)) ## Arrondi
  grd <- expand.grid(
    x = seq(from = BB["xmin"] - 2 * resolution, to = BB["xmax"] + 2 * resolution, by = resolution),
    y = seq(from = BB["ymin"] - 2 * resolution, to = BB["ymax"] + 2 * resolution, by = resolution)
  )
  sp::coordinates(grd) <- ~x + y
  sp::gridded(grd) <- TRUE
  grd@proj4string <- sp::CRS("+init=epsg:2154")
  message("Grid is created...")
  # Extraction des essences
  for (i in 1:dim(df)[1]) {
    message(paste0("Species treatment (", i, "/", dim(df)[1], ") : ", df[i, 3]))
    b <- subset(Tarifs, espar == paste0(df[i, 1], "") & c13 >= seuilCircf, select = c("numSchR", "numSchL", "c13", "geometry"))
    if (surf < 1e+06) {
      t1 <- data.frame(Moyenne = apply(b[, 1:2], 2, mean, na.rm = T), Variance = apply(b[, 1:2], 2, sd, na.rm = T))
      row.names(t1) <- c("SchR", "SchL")
      t2 <- data.frame(
        Code = df[i, 1], Nbre = df[i, 2], Essence = df[i, 3], Tarif = rownames(t1)[which.min(t1[, 2])],
        Num = floor(t1[which.min(t1[, 2]), 1] / 0.5 + 0.5) * 0.5
      )
      Tab <- rbind(Tab, t2)
    } else {
      # Moyenne sur un pas de 500 m
      t <- b %>%
        dplyr::mutate(x = floor(sf::st_coordinates(sf::st_sf(b))[, 1] / 500 + 0.5) * 500) %>%
        dplyr::mutate(y = floor(sf::st_coordinates(sf::st_sf(b))[, 2] / 500 + 0.5) * 500)
      t <- doBy::summaryBy(numSchR + numSchL ~ x + y, data = t, FUN = mean, keep.names = T)
      sp::coordinates(t) <- ~x + y
      t@proj4string <- sp::CRS("+init=epsg:2154")
      # Krigeage pour chacun des tarifs
      for (j in 1:2) {
        v <- try(gstat::variogram(t@data[, j] ~ 1, data = t, cutoff = distmax, width = 2 * resolution), silent = FALSE)
        if (class(v) != "try-error") {
          vmf <- try(gstat::fit.variogram(v, gstat::vgm(c("Exp", "Sph", "Mat", "Ste"))), silent = FALSE)
          if (class(vmf) != "try-error") {
            k <- try(gstat::krige(t@data[, j] ~ 1, locations = t, newdata = grd, model = vmf, debug.level = 0), silent = FALSE)
            if (class(k) != "try-error") {
              if (j == 1) {
                List1 <- list(raster::raster(k[, 1]))
                List2 <- list(raster::raster(k[, 2]))
              } else {
                List1 <- c(List1, list(raster::raster(k[, 1])))
                List2 <- c(List2, list(raster::raster(k[, 2])))
              }
            }
          }
        }
      }
      b1 <- raster::brick(List1)
      names(b1) <- c("SchR", "SchL")
      b2 <- raster::brick(List2)
      names(b2) <- c("SchR", "SchL")
      # Extraction par essence
      v1 <- raster::mask(b1, as(perimetre, "Spatial"))
      v2 <- raster::mask(b2, as(perimetre, "Spatial"))
      t1 <- data.frame(Moyenne = apply(v1@data@values, 2, mean, na.rm = T), Variance = apply(v2@data@values, 2, sd, na.rm = T))
      t2 <- data.frame(
        Code = df[i, 1], Nbre = df[i, 2], Essence = df[i, 3], Tarif = rownames(t1)[which.min(t1[, 2])],
        Num = round(floor(t1[which.min(t1[, 2]), 1] / 0.5 + 0.5) * 0.5, 0)
      )
      Tab <- rbind(Tab, t2)
    }
  }
  message("Calculations done finally !!!")
  if (enreg) {
    readr::write_excel_csv(Tab, paste0(dirname(fichier), "/tarifIFNSER.csv"))
    message(paste0("Result is saved in: ", paste0(dirname(fichier), "/tarifIFNSER.csv")))
  }
  out <- list(Tab)
  names(out) <- c("Tableau")
  return(out)
}