#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.