R/TarifsEmerge.R

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

#' @title Tarifs EMERGE
#'
#' @description Fonction pour calculer les volumes geometriques des especes francaises issus du projet ANR EMERGE
#'
#' @details Calcule le volume tige, houppier ou total d'un arbre selon une decoupe donnee
#' pour la quasi totalite des espece francaises ou sinon pour des groupes d'essences
#'
#' Ces tarifs sont issus du projet ANR Emerge 2009-2014 pilote par R&D ONF
#'
#' @name TarifsEmerge
#'
#' @aliases TarEmerge
#'
#' @param c130 Circonference(s) a 1,30 m (en cm) : vecteur \code{numeric}
#' @param htot Hauteur(s) totale(s) (en m) : vecteur \code{numeric}
#' @param hdec Hauteur(s) de decoupe(s) (au sens IFN) (en m).
#' Vaut \code{NULL} par defaut : pas besoin de \code{hdec} pour calculer le volume total.
#' @param dec Decoupe pour calculer le volume cumule depuis le bas de l'arbre vers
#'  cette decoupe. Vaut \code{'7'} par defaut, c'est a dire decoupe bois-fort.
#' @param ess Le code de l'essence d'apres le code de l'ONF.  Vaut \code{NULL} par defaut.
#' Ne peut etre renseigne en meme temps que \code{nomess}.
#' Le vecteur \code{ess} est recycle s'il est moins long que \code{c130}.
#' @param nomess Le nom detaille de l'essence. Vaut \code{NULL} par defaut.
#' Ne peut etre renseigne en meme temps que \code{ess}.
#' Le vecteur \code{nomess} est recycle s'il est moins long que \code{c130}.
#' @param espar Le code ifn  de l'essence. Vaut \code{NULL} par defaut.
#' Ne peut etre renseigne en meme temps que \code{ess} et \code{nomess}.
#' @param typevol Type de volume renvoye par le tarif : volume de la
#' Peut etre abrege. Vaut \code{'total'} par defaut.
#' @param perchebrin Indique si le calcul porte sur des perches et brins, par defaut egal a FALSE.
#' Si TRUE, on ne vérifie pas diam/2 > dec.
#'
#' @return Un vecteur donnant le volume (en m3) recherche
#'
#' @note La liste des parametres et des tarifs disponibles (ess, nomess, etc.) est fournie
#' dans le \code{\link[base]{data.frame}}: ListTarEmerge.
#'
#' Attention, tous les tarifs n'ont pas de code \code{ess} dans les codes de l'ONF.
#' Pour les utiliser, il faut imperativement utiliser \code{nomess}.
#'
#' Des \code{NA} sont permis pour c130, htot et hdec, auquel cas \code{v = NA}
#'
#' @author Christine Deleuze (portage sous R par Francois Morneau)
#'
#' @references Projet ANR-08-BIOE-003 EMERGE
#'
#' Rapports de C. Deleuze et publications ONF :
#'
#' Deleuze, C., Morneau, F., Renaud, J.-P., Rivoire, M., Santenoise, P., Longuetaud, F., Mothe, F.,
#'   Herve, J.-C., Vallet, P. & Vivien, Y.
#'   Estimer le volume total d'un arbre, quelles que soient l'essence, la taille, la sylviculture, la station,
#'   Rendez-Vous Techniques, 2014, 44, 22-32
#'
#' Deleuze, C., Morneau, F., Renaud, J.-P., Rivoire, M., Santenoise, P., Longuetaud, F., Mothe, F., Herv?, J.-C. & Vivien, Y.
#'   Estimation harmonisee du volume de tige a differentes decoupes,
#'   Rendez-Vous Techniques, 2014, 44, 33-42
#'
#' @family Tarifs
#'
#' @keywords function
#' @export TarEmerge
#' @examples
#' \dontrun{
#' require(graphics)
#'
#' ## Utilisation de base
#' TarEmerge(c130 = 50, htot = 15, dec = 7,  hdec = 10,  nomess = 'Fagus sylvatica', typevol = 'total')
#' TarEmerge(c130 = 50, htot = 15, dec = 7, hdec = 10,  ess = 'HET', typevol = 'tige')
#' TarEmerge(c130 = 50, htot = 15, dec = 7, hdec = 10,  espar = '09', typevol = 'tige')
#' }
#'

TarEmerge <- function(c130, htot, hdec = NULL, dec = 7, ess = NULL, espar = NULL, nomess = NULL, typevol = "total", perchebrin = FALSE) {
  if (!is.numeric(c130) || !is.numeric(htot) || any(c130 < 0, na.rm = TRUE) || any(htot < 0, na.rm = TRUE)) {
    stop("c130 and htot must be positive numeric (in cm and m respectively)")
  }
  nobs <- length(c130)
  if (max(c130, na.rm = TRUE) > 1000) {
    warning("'c 130' important, check the units (in cm)")
  }
  if (max(htot, na.rm = TRUE) > 50) {
    warning("'hot' important, check the units (in m)")
  }
  typevol <- match.arg(tolower(typevol), c("tige", "total", "houp"), several.ok = FALSE)
  if (!is.null(ess) & !is.null(nomess)) {
    stop("One argument for species: ess OR nomess")
  }
  if (!is.null(espar) & !is.null(nomess)) {
    stop("One argument for species: espar OR nomess")
  }
  if (!is.null(espar) & !is.null(ess)) {
    stop("One argument for species: espar OR ess")
  }
  if (is.null(espar) & is.null(nomess) & !is.null(ess)) {
    if (sum(is.na(ess)) > 0) {
      stop("Missing values of 'ess'")
    }
    if (length(ess) < nobs) {
      ess <- rep_len(ess, length.out = nobs)
    }
    ## ATTENTION ---- faut retenir que les valeurs uniques ---------
    liste <- unique(ListTarEmerge[!is.na(ListTarEmerge$ess), c("ess", "a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
    matcoefs <- as.matrix(liste[charmatch(ess, liste$ess), c("a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
  } else {
    if (is.null(espar) & !is.null(nomess) & is.null(ess)) {
      if (sum(is.na(nomess)) > 0) {
        stop("Missing values of 'nomess'")
      }
      if (length(nomess) < nobs) {
        nomess <- rep_len(nomess, length.out = nobs)
      }
      liste <- unique(ListTarEmerge[!is.na(ListTarEmerge$nomess), c("nomess", "a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
      matcoefs <- as.matrix(liste[charmatch(nomess, liste$nomess), c("a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
    } else {
      if (!is.null(espar) & is.null(nomess) & is.null(ess)) {
        if (sum(is.na(espar)) > 0) {
          stop("Missing values of 'espar'")
        }
        if (length(espar) < nobs) {
          espar <- rep_len(espar, length.out = nobs)
        }
        liste <- unique(ListTarEmerge[!is.na(ListTarEmerge$ESPAR), c("ESPAR", "a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
        matcoefs <- as.matrix(liste[charmatch(espar, liste$ESPAR), c("a1", "b1", "c1", "a2", "b2", "c2", "d2", "a3", "b3")])
      }
    }
  }
  if (any(rownames(matcoefs) == "NA")) {
    warning("Some levels of 'ess' are not valid:\n     ", paste(ess[rownames(matcoefs) == "NA"], collapse = "\n"), "\n  A zero volume will be returned.")
  }
  # calcul du volume total decoupe 0
  logis_hdec <- log(hdec / (htot - hdec))
  c13 <- c130 / 100
  cDec <- dec * pi / 100
  hdn <- sqrt(c13) / htot
  hsurd <- htot / (c13)
  corr130 <- (1 - 1.3 / htot)
  # calcul du volume de cylindre equivalent a même hauteur et base surface terriere, en tenant compte de la correction bas de tige
  volCyl <- (c13) ^ 2 * htot / (4 * pi * corr130 ^ 2)
  # on propose le calcul de tout au debut et on conditionne dans un second temps
  .Xform <- matcoefs[, c("a1", "b1", "c1")] * matrix(c(rep(1, nobs), hdn, hsurd), ncol = 3, byrow = FALSE)
  formTot <- .rowSums(.Xform, n = ncol(.Xform), m = nrow(.Xform), na.rm = TRUE) # plus rapide
  # formTige <- diag(matcoefs[, c('a2', 'b2', 'c2', 'd2')] %*% matrix(c(rep(1, length(c130)), logis_hdec, hdn, 1/c13), nrow = 4, byrow = TRUE))
  .XFCT <- matcoefs[, c("a2", "b2", "c2", "d2")] * matrix(c(rep(1, nobs), logis_hdec, hdn, 1 / c13), ncol = 4, byrow = FALSE)
  FCT <- .rowSums(.XFCT, n = ncol(.XFCT), m = nrow(.XFCT), na.rm = TRUE) # plus rapide
  formTige <- formTot * FCT
  alpha <- matcoefs[, c("b3")] # c est inverse dans la listeEmerge!!
  beta <- matcoefs[, c("a3")]
  Cmax <- (c13 / corr130) * ((3 - beta) * (1 / FCT - 1) / (3 * alpha)) ^ (1 / 3)
  ## comme c'est un vecteur, il faut un ifelse
  decHoup <- ifelse(Cmax > cDec, (1 - ((cDec / c13 * corr130) * (((3 - beta) / (3 * alpha)) * (1 / FCT - 1)) ^ (-1 / 3)) ^ (3 - beta)), 0)
  volHoup <- (formTot - formTige) * volCyl * decHoup
  # Rq Christine ajouter un test ici sur dec / c130 : il faudrait que dec < diam/2 if(any((c130/pi)/2 < dec) ) stop('Au moins une tige est trop petite pour la decoupe
  # indiquee')
  # Modif Pascal Obstetar
  # volTige <- ifelse((c130/pi)/2 > dec, (formTige * volCyl * (1 - (cDec * corr130/c13)^3)), NA)
  if (!perchebrin) {
    if (any(c130 / pi / 2 > dec)) {
      volTige <- formTige * volCyl * (1 - (cDec * corr130 / c13) ^ 3)
    } else {
      volTige <- NA
      warning(paste(" This species:", ess, ", Htot =", htot, ", Hdec =", hdec, ", with Diam/2 =", round(c130 / pi / 2, 0), "<", dec, "= dec is too little for calculation!"))
    }
  } else {
    volTige <- formTige * volCyl * (1 - (cDec * corr130 / c13) ^ 3)
  }
  # ensuite on renvoie les infos, avec juste une interro sur la fonction if ? a-t-on une fonction case ?
  if (typevol == "tige") {
    return(volTige)
  } else {
    if (typevol == "houp") {
      return(volHoup)
    } else {
      volTot <- volTige + volHoup # on est bien a la même decoupe dec
      return(volTot)
    }
  }
}