R/MapSpeciesSER.R

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

#' Cartographie du volume d'une essence ou des essences par sylvoecoregion.
#'
#' @description La fonction renvoie une carte de la France decoupee en sylvoecoregion. La couleur des sylvoecoregions
#' est proportionnelle au volume moyen a l'hectare de l'essence retenue (rouge d'autant plus intense que
#' le volume a l'hectare est eleve). Les points bleus permettent de localiser les placettes IFN.
#' Leur densite permet d'estimer la precision du volume.
#'
#' @return La fonction renvoie une carte des sylvoecoregions avec un niveau de rouge proportionnel
#' au volume de l'essence retenue et un graphique boîte a moustaches (boxplot) de ces valeurs.
#' Les sylvoecoregions qui ne contiennent pas l'essence retenue sont en gris fonce.
#' Les placettes IFN utilisees sont en bleu.
#'
#' @param essence = c('02', '09') par defaut. Vecteur code(s) essence(s) de l'IFN. Par defaut ce code est egal
#' aux codes EPSAR de l'IFN '02' et '09', ils correspondent au chêne pedoncule et au hêtre.
#' @param annee = c(2005:2016) par defaut. Vecteur annee des mesures IFN consideree.
#' Par defaut ce vecteur represente les annees 2005 a 2016.
#'
#' @export MapSpeciesSER
#'
#' @author Pascal Obstetar
#'
#' @examples
#' \dontrun{
#' test <- MapSpeciesSER()
#' test$Carte
#' test$Box
#' MapSpeciesSER(essence = c('02','03','09','11','52','61','62', '54','64'))
#' }
#'

MapSpeciesSER <- function(essence = c("02", "09"), annee = c(2005:2016)) {
  vecteur_annee <- annee
  donnees <- system.file("extdata/IFN", package = "gftools")
  CodesEssIFN <- gftools::getData("IFNCODE") %>% dplyr::filter(donnee %in% "ESPAR")
  # 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", "v", "w")])
    arb <- rbind(arb, arbre)
  }
  ser <- st_read(system.file("extdata/shapes/SER/ser_l93.shp", package = "gftools"), quiet = TRUE)

  t1 <- arb %>%
    dplyr::filter(espar %in% essence) %>%
    dplyr::group_by(idp, espar) %>%
    dplyr::summarise(Vol = sum(as.numeric(v) * as.numeric(w)))

  suppressWarnings(
    suppressMessages(
      t2 <- t1 %>%
        dplyr::left_join(plac[, c("idp", "ser")], by = "idp") %>%
        dplyr::mutate(Vol = replace(Vol, which(is.na(Vol)), 0)) %>%
        dplyr::group_by(ser, espar) %>%
        dplyr::summarise(Vol = mean(Vol)) %>%
        dplyr::left_join(CodesEssIFN[, 2:3], by = c(espar = "code"))
    )
  )

  shp <- merge(plac, t1, by = "idp")
  t3 <- as.data.frame(shp)
  t3 <- subset(t3, !is.na(t3$Vol))

  t4 <- reshape2::dcast(t2, ser ~ espar, value.var = "Vol")

  suppressWarnings(
    suppressMessages(
      ser.f <- ser %>%
        dplyr::left_join(t4, by = c(codeser = "ser")) %>%
        dplyr::select(essence, geometry) %>%
        tidyr::gather(Species, Volume, -geometry) %>%
        dplyr::left_join(CodesEssIFN[, 2:3], by = c(Species = "code"))
    )
  )

  g <- ggplot() + geom_sf(data = ser.f, aes(fill = Volume), color = NA) +
    facet_wrap(~libelle) +
    scale_fill_gradient(low = "yellow", high = "red") + xlab("") + ylab("") +
    geom_point(data = t3, aes(x = xl93, y = yl93), size = 0.1, alpha = 0.01, color = "blue") +
    guides(fill = guide_legend(title = "Vol. (m3/ha)"))

  h <- ggplot(t2, aes(libelle, Vol)) +
    xlab("") + ylab("Volume (m3/ha)") +
    geom_boxplot(outlier.colour = "red", outlier.size = 1, notch = TRUE) +
    stat_summary(
      fun.data = "mean_sdl",
      geom = "errorbar", width = 0.3, colour = "green", alpha = 0.5
    ) +
    stat_summary(fun.y = mean, geom = "point", colour = "green") +
    stat_summary(
      fun.y = mean, fun.ymax = max,
      fun.ymin = min, geom = "label", color = "darkgreen", hjust = -0.5, aes(alpha = 10, label = round(..y.., digits = 0)), show.legend = FALSE
    ) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  out <- list(g, h)
  names(out) <- c("Carte", "Boxplot")
  message("Graph are realized!")
  return(out)
}