R/AccDIFNSER.R

Defines functions AccDIFNSER

Documented in AccDIFNSER

#' Accroissement sur le diametre moyen en cm par an sur 5 ans
#'
#' @description La fonction renvoie les accroissements moyen sur les 5 dernieres annees
#' sur le diametre par essence et classe de diametre en utilisant la base de donnees arbres de l'IFN.
#' La fonction necessite en entree un shape correspondant au perimetre retenu : foret, massif,
#' sylvoecoregion, etc.
#'
#' @return La fonction renvoie un tableau et un graphique des accroissements sur le diametre par essence
#' par periode et classe de diametre.
#'
#' @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
#' @param periode = NULL par defaut, couvre donc l'ensemble de la periode du vecteur annee, si 1 alors resultat calcule annee par annee
#' @param tailleBuffer = Taille du buffer en metre. Par defaut elle est egale a 3000 m.
#' @param seuilCircf = 50 cm par defaut, seuil minimal de circonference en cm en dessous duquel l'arbre n'est pas retenu.
#' @param seuilNb = 10 tiges par defaut, seuil minimal du nombre d'arbres dans la base IFN pour qu'une essence soit retenue.
#' @param useSer = TRUE par defaut, argument permettant de choisir si le calcul des accroissements doit se faire
#' au sein d'une meme sylvoecoregion.
#' @param enreg = FALSE par defaut, argument permettant de choisir l'enregistrement ou pas du tableau au format csv.
#'
#' @export AccDIFNSER
#'
#' @author Pascal Obstetar
#'
#' @examples
#' \dontrun{
#' AccDIFNSER(fichier = '~/ec21.shp', enreg = T, seuilNb = 30)
#' AccDIFNSER(fichier = '~/ec21.shp', enreg = T, seuilNb = 100)
#' AccDIFNSER(fichier = '~/ec21.shp', enreg = T, seuilNb = 500)
#' test <- AccDIFNSER(fichier = '~/ec21.shp', enreg = F, seuilNb = 100, periode=3)
#' test$Tableau
#' test$Effectif
#' test$Graphe
#' AccDIFNSER(fichier = '~/ec21.shp', enreg = F, seuilNb = 500, periode=6)
#' AccDIFNSER(fichier = '~/ec21.shp', enreg = F, seuilNb = 1000, periode=6)
#' }
#'

AccDIFNSER <- function(fichier = NULL, annee = c(2005:2016), periode = NULL, tailleBuffer = 1000, seuilCircf = 50, seuilNb = 1000, useSer = T, 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
  periode <- dplyr::if_else(is.null(periode), max(annee) - min(annee) + 1, periode)
  dossier <- system.file("extdata/IFN", package = "gftools")
  CodesEssIFN <- gftools::getData("IFNCODE") %>% dplyr::filter(donnee %in% "ESPAR")
  # Extract 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 = dossier)
    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 = dossier)
    arbre <- data.frame(yrs = yrs, ifn$Arb[[1]][, c("idp", "espar", "c13", "w", "ir5", "veget")])
    arb <- rbind(arb, arbre)
  }
  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)
  zone <- sf::st_buffer(perimetre, dist = tailleBuffer)
  placettes <- sf::st_intersection(sf::st_as_sf(plac, coords = c("xl93", "yl93"), crs = 2154, agr = "constant"), zone)

  AcctDs <- arb %>%
    dplyr::filter(idp %in% placettes$idp) %>%
    dplyr::filter(veget == "0") %>%
    dplyr::filter(!is.na(ir5) & c13 >= seuilCircf) %>%
    dplyr::select(espar, c13, ir5, yrs) %>%
    dplyr::mutate(Diam = round(c13 / pi, 0)) %>%
    dplyr::mutate(Classe = floor(c13 / pi / 5 + 0.5) * 5) %>%
    dplyr::mutate(espar = as.character(espar)) %>%
    dplyr::mutate(Population = 1 + (yrs - min(yrs)) %/% periode) %>%
    dplyr::mutate(Periode = paste(min(arb$yrs) + (Population - 1) * periode, min(arb$yrs) + (Population * periode) - 1, sep = " - ")) %>%
    dplyr::left_join(CodesEssIFN[, 2:3], by = c(espar = "code")) %>%
    dplyr::left_join(Ecorces[, c(1, 4)], by = c(espar = "codeIFN")) %>%
    dplyr::mutate(AcctD = round(ir5 * 2 / 50 / (1 - 2 * pi * b / 1000), 3))

  if (enreg) {
    readr::write_excel_csv(AcctDs, paste0(dirname(fichier), "/AccDs.csv"))
    message(paste0("Result is saved in: ", paste0(dirname(fichier), "/AccDs.csv")))
  }

  df <- AcctDs %>%
    dplyr::group_by(Periode, libelle) %>%
    dplyr::summarise(Nidp = n()) %>%
    dplyr::filter(Nidp >= seuilNb) %>%
    dplyr::arrange(libelle, Periode)

  t1 <- AcctDs %>%
    dplyr::filter(libelle %in% unique(df$libelle))

  t2 <- AcctDs
  t2$libelle <- " Toutes essences"
  t2$espar <- "99"

  t3 <- rbind(t1, t2)

  m1 <- aggregate(AcctD ~ Diam + libelle + Periode, t3, mean)
  m2 <- aggregate(Diam ~ libelle + Periode, t3, max)
  colnames(m2) <- c("libelle", "Periode", "DMax")
  mm <- m1 %>%
    dplyr::group_by(libelle, Periode) %>%
    dplyr::summarise_all(mean)
  mm <- merge(mm, m2)
  mm$sd_AcctD <- aggregate(AcctD ~ libelle + Periode, m1, sd)[, 3]
  mm <- mm %>%
    dplyr::mutate(m_plus_sd = AcctD + 1.96 * sd_AcctD) %>%
    dplyr::mutate(m_moins_sd = AcctD - 1.96 * sd_AcctD) %>%
    dplyr::mutate(cv = sd_AcctD / AcctD)

  options(digits = 1)

  p <- ggplot(t3, aes(x = Diam, y = AcctD)) +
    geom_point(alpha = 0.5, col = "gray") +
    geom_boxplot(
      outlier.colour = "green", outlier.size = 2, notch = TRUE, alpha = 0.3,
      fill = "green", show.legend = FALSE
    ) +
    geom_smooth(span = 1, method = "loess", se = TRUE, color = "red", fill = "blue", alpha = 0.3, size = 0.3) +
    facet_wrap(libelle ~ Periode) +
    geom_label(data = mm, aes(x = DMax + 20, label = round(AcctD, 1), y = AcctD), col = "red") +
    geom_label(data = mm, aes(x = DMax + 15, label = round(m_plus_sd, 1), y = m_plus_sd), col = "blue", alpha = 0.1) +
    geom_label(data = mm, aes(x = DMax + 15, label = round(m_moins_sd, 1), y = m_moins_sd), col = "blue", alpha = 0.1) +
    geom_errorbar(data = mm, mapping = aes(x = DMax + 10, ymin = m_moins_sd, ymax = m_plus_sd), size = 0.5, color = "blue", width = 2) +
    geom_point(data = mm, aes(x = DMax + 10), col = "red") + xlab("Diametre") + ylab("Acct.D en cm/an/5 ans")

  ListeEss <- unique(t1[, c("espar", "libelle")])
  Xnew <- data.frame(Diam = seq(20, 100, 5))

  options(digits = 2)

  tab <- Xnew

  for (i in 1:dim(ListeEss)[1]) {
    t2 <- subset(t1, espar == ListeEss[i, 1])
    try(model <- stats::loess(AcctD ~ Diam, data = t2, span = 1, na.action = na.exclude))
    res <- data.frame(round(predict(model, newdata = Xnew), 2))
    names(res) <- ListeEss[i, 2]
    tab <- cbind(tab, res)
  }

  if (enreg) {
    readr::write_excel_csv(tab, paste0(dirname(fichier), "/tabAccDIFNSER.csv"))
    message(paste0("Result is saved in: ", paste0(dirname(fichier), "/tabAccDIFNSER.csv")))
  }
  out <- list(tab, df, p)
  names(out) <- c("Tableau", "Effectif", "Graphe")
  message("Calculation realized!")
  options(digits = 7)
  return(out)
}
pobsteta/gftools documentation built on March 28, 2020, 8:25 p.m.