R/TarifFindSch.R
In pobsteta/gftools: Geographic Data Analysis and Modeling for Forests
Defines functions
TarifFindSch
Documented in
TarifFindSch
#' Calcul des tarifs de cubage
#'
#' @description La fonction permet de trouver parmi les tarifs Schaeffer, ceux qui sont le plus adaptes au cubage
#' de la zone inventoriee.
#' La fonction necessite en entree un fichier au format csv (separateur champ Tabulation, separateur decimal point, encodage UTF-8) contenant 4 colonnes : essence,
#' diam, htot et hdec, meme si hdec est vide
#'
#' @return La fonction renvoie un tableau qui fournit par essence le numero de tarif Schaeffer rapide, ou lent
#' ainsi que les coefficients de variations associes, ainsi qu'un graphique permettant d’apprecier la variabilite
#' des numeros des tarifs retenus.
#'
#' @param enreg = FALSE par defaut, argument permettant de choisir l'enregistrement ou pas du tableau au format .csv
#' @param fichier = nom du fichier csv des data a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param mercuriale = nom du fichier csv de la mercuriale a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param clause = nom du fichier csv des clauses territoriales des decoupes appliquees a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param decemerge = 7 par defaut, indique le diametre de la decoupe fin bout pour le calcul du volume EMERGE
#' @param typvolemerge = 'total' par defaut, indique le type de volume EMERGE calculé
#' @param classearbremin = 20 par defaut, indique le diametre minimal a partir duquel les arbres sont pris en compte pour le calcul
#' @param classearbremax = 80 par defaut, indique le diametre maximal a partir duquel les arbres ne sont plus pris en compte pour le calcul
#' @param mappoint = FALSE par defaut, indique si on utilise les data IFN autour d'un point pour le calcul
#' @param tailleBuffer = NULL par defaut, taille du buffer en metre
#' @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 latitude = NULL par defaut, latitude utilisee pour les data IFN (EPSG=4326)
#' @param longitude = NULL par defaut, longitude utilisee pour les data IFN (EPSG=4326)
#' @param zonecalc = NULL par defaut, objet sf_spatialdataframe (EPSG=2154) limitant la zone de calcul a la zone SER par defaut sinon à l'objet considere
#' @param houppier = 'N' par défaut, indique si le houppier est compris dans le calcul du volume ou pas
#' @param echantillon = NULL par defaut, tableau des valeurs (htot= f(diam) de l'echantillon
#'
#' @author Bruciamacchie Max, Pascal Obstetar
#'
#' @export TarifFindSch
#'
#' @examples
#' \dontrun{
#' TarifFindSch(fichier = "~/ex_EPC_FOP_IGN.csv", enreg = F)
#' TarifFindSch(fichier = "~/EPC_SP.csv", enreg = F)
#' }
#'
TarifFindSch <- function(fichier = NULL, enreg = F, mercuriale = NULL, decemerge = 7, typvolemerge = "total",
classearbremin = 20, classearbremax = 80, mappoint = FALSE, tailleBuffer = NULL,
essence = c("02", "09"), latitude = NULL, longitude = NULL, zonecalc = NULL, houppier = "N",
clause = NULL, echantillon = NULL) {
split <- function(texte) {
strsplit(texte, " ")[[1]][1]
}
splitv <- Vectorize(split)
TarONF3v <- Vectorize(gftools::TarONF3)
if (is.null(fichier) & !mappoint) {
stop("Your fichier is NULL!")
}
# if (is.null(mercuriale) & !mappoint) {
# stop("Your mercuriale is NULL!")
# }
message("Extract mercuriale file...")
if (!is.null(mercuriale)) {
mer <- readr::read_tsv(
mercuriale,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(cdiam = readr::col_integer(), tarif = readr::col_character(), houppier = readr::col_integer(), hauteur = readr::col_double()),
col_names = T
) %>%
filter(!is.na(tarif))
} else {
mer <- data.frame(cdiam = seq(from = 10, to = 120, by = 5), tarif = rep("SR14", 23), houppier = c(rep(0, 3), rep(30, 20)), hauteur = rep(0, 23))
}
if (!is.null(clause)) {
clo <- readr::read_tsv(
clause,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(ess = readr::col_character(), dmin = readr::col_integer(), dmax = readr::col_integer(), dec = readr::col_double()),
col_names = T
)
} else {
clo <- data.frame(ess = "Defaut", dmin = 10, dmax = 200, dec = 7)
}
if (mappoint) {
vecteur_annee <- c(2008:2016)
dossier <- system.file("extdata/IFN", package = "gftools")
# 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", "htot", "hdec", "veget")])
arb <- rbind(arb, arbre)
}
if (is.null(zonecalc)) {
zonecalc <- sf::st_read(system.file("extdata/shapes/SER/ser_l93.shp", package = "gftools"), quiet = TRUE)
zonecalc <- sf::st_transform(zonecalc, crs = 2154)
point <- sf::st_as_sf(data.frame(click = "Point", lat = latitude, lng = longitude), coords = c("lng", "lat"), crs = 4326, agr = "constant")
point <- sf::st_transform(point, crs = 2154)
zone <- sf::st_join(zonecalc, point, join = st_covers) %>%
filter(click == "Point")
} else if (!is.null(tailleBuffer)) {
point <- sf::st_as_sf(data.frame(click = "Point", lat = latitude, lng = longitude), coords = c("lng", "lat"), crs = 4326, agr = "constant")
point <- sf::st_transform(point, crs = 2154)
zone <- sf::st_buffer(point, dist = tailleBuffer)
} else {
zone <- sf::st_transform(zonecalc, crs = 2154)
}
placettes <- sf::st_intersection(sf::st_as_sf(plac, coords = c("xl93", "yl93"), crs = 2154, agr = "constant"), sf::st_geometry(zone))
tab <- arb %>%
dplyr::filter(espar %in% essence) %>%
dplyr::filter(idp %in% placettes$idp) %>%
dplyr::filter(veget == "0") %>%
dplyr::select(espar, c13, htot, hdec) %>%
dplyr::mutate(diam = round(c13 / pi, 0)) %>%
dplyr::mutate(espar = as.character(espar)) %>%
dplyr::left_join(code_ifn_onf, by = c(espar = "espar")) %>%
dplyr::filter(!is.na(htot)) %>%
dplyr::mutate(esscct = ifelse(splitv(essence) %in% c("Hetre", "Chene", "Pin"), splitv(essence), fr))
} else {
message("Extract data file...")
tab <- readr::read_tsv(
fichier,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(code = readr::col_character(), diam = readr::col_integer(), htot = readr::col_double(), hdec = readr::col_double()),
col_names = T
)
tab <- tab %>%
dplyr::left_join(code_ifn_onf, by = c(code = "code")) %>%
dplyr::mutate(esscct = ifelse(splitv(essence) %in% c("Hetre", "Chene", "Pin"), splitv(essence), fr))
}
tab <- tab %>%
dplyr::mutate(Classe = as.integer(floor(diam / 5 + 0.5) * 5)) %>%
dplyr::inner_join(mer, by = c(Classe = "cdiam"))
## recherche la decoupe par defaut emerge dans les clauses ou prend decemerge
defo <- ifelse(!is.null(clause), clo[clo$ess %in% "Defaut", "dec"][[1]], decemerge)
## joint les clauses en fonction de essence et classe diam
tab1 <- tab %>%
dplyr::mutate(clo = splitv(esscct)) %>%
dplyr::left_join(clo, by = c(clo = "ess"))
# on gere les essences esscct decrites dans la mercuriale
tab2 <- tab1 %>%
dplyr::filter(!is.na(dmin)) %>%
dplyr::filter(Classe >= dmin & Classe <= dmax)
# on gere les essences esscct non decrites dans la mercuriale
tab3 <- tab1 %>%
dplyr::filter(is.na(dmin)) %>%
dplyr::select(espar, htot, hdec, diam, essence, Classe, fr, clo, tarif, houppier, hauteur) %>%
dplyr::left_join(clo, by = c(fr = "ess"))
tab4 <- dplyr::bind_rows(tab2, tab3) %>%
dplyr::mutate(defaut = defo)
tab4$decemerge <- ifelse(!is.na(tab4$dmin) & tab4$Classe >= tab4$dmin, tab4$dec, tab4$defaut)
tab <- tab4 %>%
dplyr::select(espar, essence, diam, Classe, htot, hdec, decemerge, tarif, hauteur, houppier)
if (houppier == "N") {
tab <- tab %>%
dplyr::mutate(
L_VbftigCom = as.numeric(TarONF3v(tarif = tarif, entr1 = diam, entr2 = hauteur, details = FALSE)),
L_VHouppiers = L_VbftigCom * houppier / 100,
L_Vbftot7cm = L_VbftigCom + L_VHouppiers,
L_PHouppiers = houppier / 100
)
} else {
tab <- tab %>%
dplyr::mutate(
L_Vbftot7cm = as.numeric(TarONF3v(tarif = tarif, entr1 = diam, entr2 = hauteur, details = FALSE)),
L_VHouppiers = L_Vbftot7cm * houppier / 100,
L_VbftigCom = L_Vbftot7cm - L_VHouppiers,
L_PHouppiers = houppier / 100
)
}
tab <- cbind(tab, E_VbftigCom = TarEmerge(c130 = pi * tab$diam, htot = tab$htot, hdec = tab$hdec, espar = tab$espar, typevol = "tige", dec = tab$decemerge))
tab <- cbind(tab, E_Vbftot7cm = TarEmerge(c130 = pi * tab$diam, htot = tab$htot, hdec = tab$hdec, espar = tab$espar, typevol = "total", dec = 7)) %>%
dplyr::mutate(E_PHouppiers = E_Vbftot7cm / E_VbftigCom - 1) %>%
dplyr::filter(!is.na(E_Vbftot7cm)) %>%
dplyr::mutate(E_VHouppiers = E_Vbftot7cm - E_VbftigCom)
tab1 <- tab %>%
dplyr::mutate_if(is.numeric, funs(round(., 2)))
if (sum(c("essence", "diam", "htot", "hdec") %in% names(tab)) == 4) {
if (typvolemerge == "tige") {
tab <- tab %>%
dplyr::filter(diam >= classearbremin & diam <= classearbremax) %>%
dplyr::mutate(
numSchR = E_VbftigCom / 5 * 70000 / (diam - 5) / (diam - 10) - 8,
numSchL = E_VbftigCom / 5 * 90000 / diam / (diam - 5) - 8,
numAlg = E_VbftigCom * 28000000 / (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3) - 8
)
} else {
tab <- tab %>%
dplyr::filter(diam >= classearbremin & diam <= classearbremax) %>%
dplyr::mutate(
numSchR = E_Vbftot7cm / 5 * 70000 / (diam - 5) / (diam - 10) - 8,
numSchL = E_Vbftot7cm / 5 * 90000 / diam / (diam - 5) - 8,
numAlg = E_Vbftot7cm * 28000000 / (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3) - 8
)
}
res <- tab %>%
dplyr::group_by(essence) %>%
dplyr::summarise_at(vars(numSchR, numSchL, numAlg), funs(mean, var))
res[, 5:7] <- round(res[, 5:7]^0.5 / res[, 2:4], 2)
res[, 2:4] <- round(res[, 2:4], 2)
names(res) <- c("essence", "SchR", "SchL", "Alg", "SchRcv", "SchLcv", "Algcv")
tab <- tab %>%
dplyr::left_join(res, by = "essence") %>%
dplyr::mutate(
VSchR = ((round(SchR, 0) + 8) / 14000) * (diam - 5) * (diam - 10),
VSchL = ((round(SchL, 0) + 8) / 18000) * diam * (diam - 5),
VAlg = ((round(Alg, 0) + 8) / 28000000 * (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3))
)
tab.p <- tab %>%
dplyr::select(essence, diam, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg")))
if (typvolemerge == "tige") {
tab.q <- reshape2::melt(tab, id.vars = c("diam", "essence"), measure.vars = c("L_VbftigCom", "E_VbftigCom", "VSchL", "VSchR", "VAlg"))
} else {
tab.q <- reshape2::melt(tab, id.vars = c("diam", "essence"), measure.vars = c("L_Vbftot7cm", "E_Vbftot7cm", "VSchL", "VSchR", "VAlg"))
}
names(tab.q) <- c("diam", "essence", "Type", "Vol")
###################################################################
# ça plante avec shinyproxy.io !!! bug dplyr memory mapped
# tab.r <- tab %>%
# dplyr::group_by(essence, Classe) %>%
# dplyr::summarise_at(vars(L_VbftigCom,L_VHouppiers,E_VbftigCom,E_VHouppiers), mean, na.rm=TRUE)
tab.r1 <- aggregate(L_VbftigCom ~ essence + Classe, tab, mean)
tab.r2 <- aggregate(L_VHouppiers ~ essence + Classe, tab, mean)
tab.r3 <- aggregate(E_VbftigCom ~ essence + Classe, tab, mean)
tab.r4 <- aggregate(E_VHouppiers ~ essence + Classe, tab, mean)
tab.r <- bind_cols(tab.r1, tab.r2[3], tab.r3[3], tab.r4[3]) %>%
arrange(essence, Classe)
###################################################################
tab.r <- reshape2::melt(tab.r, id.vars = c("Classe", "essence"), measure.vars = c("L_VbftigCom", "L_VHouppiers", "E_VbftigCom", "E_VHouppiers")) %>%
dplyr::mutate(Type = substr(variable, 1, 1), variable = substr(variable, 3, 12))
names(tab.r) <- c("Classe", "essence", "tarif", "Vol", "Type")
tab.r$tarif <- factor(tab.r$tarif, levels = c("VHouppiers", "VbftigCom"))
m1 <- aggregate(Num ~ essence + Type, tab.p, mean)
m1$var <- aggregate(Num ~ essence + Type, tab.p, var)[, 3]
m1$sd <- aggregate(Num ~ essence + Type, tab.p, sd)[, 3]
m1$cv <- round(m1$var^0.5 / m1$Num, 2)
m12 <- aggregate(Num ~ essence + Type, tab.p, median)
colnames(m12) <- c("essence", "Type", "Median")
m11 <- aggregate(Num ~ essence + Type, tab.p, FUN = function(x) quantile(x, probs = 0.25))
colnames(m11) <- c("essence", "Type", "Q1")
m13 <- aggregate(Num ~ essence + Type, tab.p, FUN = function(x) quantile(x, probs = 0.75))
colnames(m13) <- c("essence", "Type", "Q3")
m2 <- aggregate(diam ~ essence + Type, tab.p, max)
colnames(m2) <- c("essence", "Type", "DMax")
mm <- m1 %>%
dplyr::group_by(essence, Type) %>%
dplyr::summarise_all(mean)
mm <- merge(mm, m2)
mm <- merge(mm, m11)
mm <- merge(mm, m12)
mm <- merge(mm, m13)
mm <- mm %>%
dplyr::mutate(m_plus_sd = Num + 1.96 * sd) %>%
dplyr::mutate(m_moins_sd = Num - 1.96 * sd)
message("Plot...")
# EMERGE + LOCAL diagramme empile
p0 <- ggplot(tab.r, aes(x = Type, y = Vol, fill = Type, alpha = tarif)) +
scale_fill_manual(values = c("red", "darkgreen")) +
geom_bar(stat = "identity", position = "stack") +
facet_grid(essence ~ Classe)
# Boxplot
p1 <- ggplot(tab.p, aes(x = diam, y = Num)) + geom_point(alpha = 0.1) +
facet_grid(essence ~ Type) +
geom_boxplot(outlier.colour = "green", outlier.size = 2, notch = TRUE, alpha = 0.1, fill = "green") +
geom_label(data = mm, aes(x = DMax + 15, label = round(Num, 0), y = Num), col = "red") +
geom_label(data = mm, aes(x = DMax + 15, label = round(m_plus_sd, 0), y = m_plus_sd), col = "blue", alpha = 0.1) +
geom_label(data = mm, aes(x = DMax + 15, label = round(m_moins_sd, 0), 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_smooth(method = "lm") +
geom_point(data = mm, aes(x = DMax + 10), color = "red")
# EMERGE + LOCAL courbe lineaire
p2 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("darkgreen", "red"))
if (typvolemerge == "tige") {
p2 <- p2 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom", "L_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom", "L_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5)
} else {
p2 <- p2 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm", "L_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm", "L_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5)
}
# EMERGE + VSchR
p3 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p3 <- p3 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchR"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchR")
} else {
p3 <- p3 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchR"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchR")
}
# EMERGE + VSchL
p4 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p4 <- p4 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchL"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchL")
} else {
p4 <- p4 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchL"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchL")
}
# EMERGE + VAlg
p5 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p5 <- p5 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VAlg"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VAlg")
} else {
p5 <- p5 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VAlg"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VAlg")
}
# Htot= f(E_Vbftig | Vbftot)
if (typvolemerge == "tige") {
tab.h <- tab %>%
dplyr::select(essence, diam, Classe, htot, E_VbftigCom, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam, -Classe, -htot, -E_VbftigCom) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg"))) %>%
dplyr::mutate(Classe = factor(Classe))
} else {
tab.h <- tab %>%
dplyr::select(essence, diam, Classe, htot, E_Vbftot7cm, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam, -Classe, -htot, -E_Vbftot7cm) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg"))) %>%
dplyr::mutate(Classe = factor(Classe))
}
p <- c(0.10, 0.35, 0.65, 0.90)
p_names <- purrr::map_chr(p, ~paste0(.x*100))
p_fh <- purrr::map(p, ~partial(quantile, probs = .x, na.rm = TRUE)) %>%
purrr::set_names(nm = paste0("h",p_names))
p_fn <- purrr::map(p, ~partial(quantile, probs = .x, na.rm = TRUE)) %>%
purrr::set_names(nm = paste0("n",p_names))
tab.h1 <- tab.h %>%
group_by(essence, Classe, Type) %>%
summarize_at(vars(htot), funs(!!!p_fh))
tab.h2 <- tab.h %>%
group_by(essence, Classe, Type) %>%
summarize_at(vars(Num), funs(!!!p_fn))
tab.n <- merge(tab.h1, tab.h2)
tab.h3 <- tab.h %>%
group_by(essence, Type) %>%
summarize_at(vars(htot), funs(!!!p_fh)) %>%
mutate_if(is.numeric, round)
m3 <- aggregate(Num ~ essence + Type + Classe, tab.h, mean)
m3$var <- aggregate(Num ~ essence + Type + Classe, tab.h, var)[, 4]
m3$sd <- aggregate(Num ~ essence + Type + Classe, tab.h, sd)[, 4]
m3$cv <- round(m3$var^0.5 / m3$Num, 2)
if (typvolemerge == "tige") {
m4 <- aggregate(E_VbftigCom ~ essence + Type + Classe, tab.h, mean)
} else {
m4 <- aggregate(E_Vbftot7cm ~ essence + Type + Classe, tab.h, mean)
}
colnames(m4) <- c("essence", "Type", "Classe", "VMean")
m5 <- aggregate(htot ~ essence + Type + Classe, tab.h, mean)
colnames(m5) <- c("essence", "Type", "Classe", "HMean")
mmm <- m3 %>%
dplyr::group_by(essence, Type, Classe) %>%
dplyr::summarise_all(mean)
mmm <- merge(mmm, m4)
mmm <- merge(mmm, m5)
mmm <- mmm %>%
dplyr::mutate(m_plus_sd = Num + 1.96 * sd) %>%
dplyr::mutate(m_moins_sd = Num - 1.96 * sd)
# Schaeffer rapide
if (typvolemerge == "tige") {
p6 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchR") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p6 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchR") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
# Schaeffer lent
if (typvolemerge == "tige") {
p7 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchL") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p7 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchL") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
# Algan
if (typvolemerge == "tige") {
p8 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numAlg") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p8 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numAlg") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
if (!is.null(echantillon)) {
# Calcul les centroides
centroids <- aggregate(cbind(echantillon$diam, echantillon$htot, echantillon$Num) ~ essence + Type, echantillon, mean)
f <- function(z) qt(0.025, df = length(z) - 1, lower.tail = F) * sd(z) / sqrt(length(z)) # 95% confidence
# f <- function(z)sd(z)/sqrt(length(z)) # function to calculate std.err
se <- aggregate(cbind(se.x = echantillon$diam, se.y = echantillon$htot, se.n = echantillon$Num) ~ essence + Type, echantillon, f)
centroids <- merge(centroids, se) # add std.err column to centroids
names(centroids) <- c("essence", "Type", "diam", "htot", "Num", "se.x", "se.y", "se.n")
}
# Numero de tarif numSchR
if (is.null(echantillon)) {
p9 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)), label = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchR") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p9 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)), label = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchR") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numSchR")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numSchR")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numSchR")), ]
p12 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numSchR") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
# Numero de tarif numSchL
if (is.null(echantillon)) {
p10 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchL") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p10 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchL") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numSchL")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numSchL")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numSchL")), ]
p13 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numSchL") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
# Numero de tarif numAlg
if (is.null(echantillon)) {
p11 <- ggplot(tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numAlg") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p11 <- ggplot(tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numAlg") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numAlg")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numAlg")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numAlg")), ]
p14 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numAlg") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
if (enreg) {
readr::write_excel_csv(tab, paste0(dirname(fichier), "/numTarifSchaeffer.csv"))
message(paste0("Result is saved in: ", paste0(dirname(fichier), "/numTarifSchaeffer.csv")))
}
if (mappoint) {
out <- list(tab1, res, tab.p, mm, tab.h, tab.n, p1, p2, p3, p4, p5, p0, placettes, zone, p6, p7, p8, p9, p10, p11, p12, p13, p14)
names(out) <- c(
"Tableau1", "Tableau2", "Tableau3", "Tableau4", "Tableau5", "Tableau6",
"Graphe1", "Graphe2", "Graphe3", "Graphe4", "Graphe5", "Graphe6",
"Placettes", "Zone", "Graphe7", "Graphe8", "Graphe9", "Graphe10",
"Graphe11", "Graphe12", "Graphe13", "Graphe14", "Graphe15"
)
} else {
out <- list(tab1, res, tab.p, mm, tab.h, tab.n, p1, p2, p3, p4, p5, p0, p6, p7, p8, p9, p10, p11, p12, p13, p14)
names(out) <- c(
"Tableau1", "Tableau2", "Tableau3", "Tableau4", "Tableau5", "Tableau6",
"Graphe1", "Graphe2", "Graphe3", "Graphe4", "Graphe5", "Graphe6",
"Graphe7", "Graphe8", "Graphe9", "Graphe10", "Graphe11", "Graphe12",
"Graphe13", "Graphe14", "Graphe15"
)
}
message("Calculation realized!")
return(out)
} else {
message("The file must contain the columns 'essence', 'diam', 'htot' and 'hdec' (even empty)!")
}
}
pobsteta/gftools documentation built on March 28, 2020, 8:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions TarifFindSch
Documented in TarifFindSch
#' Calcul des tarifs de cubage
#'
#' @description La fonction permet de trouver parmi les tarifs Schaeffer, ceux qui sont le plus adaptes au cubage
#' de la zone inventoriee.
#' La fonction necessite en entree un fichier au format csv (separateur champ Tabulation, separateur decimal point, encodage UTF-8) contenant 4 colonnes : essence,
#' diam, htot et hdec, meme si hdec est vide
#'
#' @return La fonction renvoie un tableau qui fournit par essence le numero de tarif Schaeffer rapide, ou lent
#' ainsi que les coefficients de variations associes, ainsi qu'un graphique permettant d’apprecier la variabilite
#' des numeros des tarifs retenus.
#'
#' @param enreg = FALSE par defaut, argument permettant de choisir l'enregistrement ou pas du tableau au format .csv
#' @param fichier = nom du fichier csv des data a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param mercuriale = nom du fichier csv de la mercuriale a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param clause = nom du fichier csv des clauses territoriales des decoupes appliquees a importer (separateur champ Tabulation, separateur decimal point, encodage UTF-8)
#' @param decemerge = 7 par defaut, indique le diametre de la decoupe fin bout pour le calcul du volume EMERGE
#' @param typvolemerge = 'total' par defaut, indique le type de volume EMERGE calculé
#' @param classearbremin = 20 par defaut, indique le diametre minimal a partir duquel les arbres sont pris en compte pour le calcul
#' @param classearbremax = 80 par defaut, indique le diametre maximal a partir duquel les arbres ne sont plus pris en compte pour le calcul
#' @param mappoint = FALSE par defaut, indique si on utilise les data IFN autour d'un point pour le calcul
#' @param tailleBuffer = NULL par defaut, taille du buffer en metre
#' @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 latitude = NULL par defaut, latitude utilisee pour les data IFN (EPSG=4326)
#' @param longitude = NULL par defaut, longitude utilisee pour les data IFN (EPSG=4326)
#' @param zonecalc = NULL par defaut, objet sf_spatialdataframe (EPSG=2154) limitant la zone de calcul a la zone SER par defaut sinon à l'objet considere
#' @param houppier = 'N' par défaut, indique si le houppier est compris dans le calcul du volume ou pas
#' @param echantillon = NULL par defaut, tableau des valeurs (htot= f(diam) de l'echantillon
#'
#' @author Bruciamacchie Max, Pascal Obstetar
#'
#' @export TarifFindSch
#'
#' @examples
#' \dontrun{
#' TarifFindSch(fichier = "~/ex_EPC_FOP_IGN.csv", enreg = F)
#' TarifFindSch(fichier = "~/EPC_SP.csv", enreg = F)
#' }
#'
TarifFindSch <- function(fichier = NULL, enreg = F, mercuriale = NULL, decemerge = 7, typvolemerge = "total",
classearbremin = 20, classearbremax = 80, mappoint = FALSE, tailleBuffer = NULL,
essence = c("02", "09"), latitude = NULL, longitude = NULL, zonecalc = NULL, houppier = "N",
clause = NULL, echantillon = NULL) {
split <- function(texte) {
strsplit(texte, " ")[[1]][1]
}
splitv <- Vectorize(split)
TarONF3v <- Vectorize(gftools::TarONF3)
if (is.null(fichier) & !mappoint) {
stop("Your fichier is NULL!")
}
# if (is.null(mercuriale) & !mappoint) {
# stop("Your mercuriale is NULL!")
# }
message("Extract mercuriale file...")
if (!is.null(mercuriale)) {
mer <- readr::read_tsv(
mercuriale,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(cdiam = readr::col_integer(), tarif = readr::col_character(), houppier = readr::col_integer(), hauteur = readr::col_double()),
col_names = T
) %>%
filter(!is.na(tarif))
} else {
mer <- data.frame(cdiam = seq(from = 10, to = 120, by = 5), tarif = rep("SR14", 23), houppier = c(rep(0, 3), rep(30, 20)), hauteur = rep(0, 23))
}
if (!is.null(clause)) {
clo <- readr::read_tsv(
clause,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(ess = readr::col_character(), dmin = readr::col_integer(), dmax = readr::col_integer(), dec = readr::col_double()),
col_names = T
)
} else {
clo <- data.frame(ess = "Defaut", dmin = 10, dmax = 200, dec = 7)
}
if (mappoint) {
vecteur_annee <- c(2008:2016)
dossier <- system.file("extdata/IFN", package = "gftools")
# 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", "htot", "hdec", "veget")])
arb <- rbind(arb, arbre)
}
if (is.null(zonecalc)) {
zonecalc <- sf::st_read(system.file("extdata/shapes/SER/ser_l93.shp", package = "gftools"), quiet = TRUE)
zonecalc <- sf::st_transform(zonecalc, crs = 2154)
point <- sf::st_as_sf(data.frame(click = "Point", lat = latitude, lng = longitude), coords = c("lng", "lat"), crs = 4326, agr = "constant")
point <- sf::st_transform(point, crs = 2154)
zone <- sf::st_join(zonecalc, point, join = st_covers) %>%
filter(click == "Point")
} else if (!is.null(tailleBuffer)) {
point <- sf::st_as_sf(data.frame(click = "Point", lat = latitude, lng = longitude), coords = c("lng", "lat"), crs = 4326, agr = "constant")
point <- sf::st_transform(point, crs = 2154)
zone <- sf::st_buffer(point, dist = tailleBuffer)
} else {
zone <- sf::st_transform(zonecalc, crs = 2154)
}
placettes <- sf::st_intersection(sf::st_as_sf(plac, coords = c("xl93", "yl93"), crs = 2154, agr = "constant"), sf::st_geometry(zone))
tab <- arb %>%
dplyr::filter(espar %in% essence) %>%
dplyr::filter(idp %in% placettes$idp) %>%
dplyr::filter(veget == "0") %>%
dplyr::select(espar, c13, htot, hdec) %>%
dplyr::mutate(diam = round(c13 / pi, 0)) %>%
dplyr::mutate(espar = as.character(espar)) %>%
dplyr::left_join(code_ifn_onf, by = c(espar = "espar")) %>%
dplyr::filter(!is.na(htot)) %>%
dplyr::mutate(esscct = ifelse(splitv(essence) %in% c("Hetre", "Chene", "Pin"), splitv(essence), fr))
} else {
message("Extract data file...")
tab <- readr::read_tsv(
fichier,
locale = readr::locale(encoding = "UTF-8", decimal_mark = "."),
readr::cols(code = readr::col_character(), diam = readr::col_integer(), htot = readr::col_double(), hdec = readr::col_double()),
col_names = T
)
tab <- tab %>%
dplyr::left_join(code_ifn_onf, by = c(code = "code")) %>%
dplyr::mutate(esscct = ifelse(splitv(essence) %in% c("Hetre", "Chene", "Pin"), splitv(essence), fr))
}
tab <- tab %>%
dplyr::mutate(Classe = as.integer(floor(diam / 5 + 0.5) * 5)) %>%
dplyr::inner_join(mer, by = c(Classe = "cdiam"))
## recherche la decoupe par defaut emerge dans les clauses ou prend decemerge
defo <- ifelse(!is.null(clause), clo[clo$ess %in% "Defaut", "dec"][[1]], decemerge)
## joint les clauses en fonction de essence et classe diam
tab1 <- tab %>%
dplyr::mutate(clo = splitv(esscct)) %>%
dplyr::left_join(clo, by = c(clo = "ess"))
# on gere les essences esscct decrites dans la mercuriale
tab2 <- tab1 %>%
dplyr::filter(!is.na(dmin)) %>%
dplyr::filter(Classe >= dmin & Classe <= dmax)
# on gere les essences esscct non decrites dans la mercuriale
tab3 <- tab1 %>%
dplyr::filter(is.na(dmin)) %>%
dplyr::select(espar, htot, hdec, diam, essence, Classe, fr, clo, tarif, houppier, hauteur) %>%
dplyr::left_join(clo, by = c(fr = "ess"))
tab4 <- dplyr::bind_rows(tab2, tab3) %>%
dplyr::mutate(defaut = defo)
tab4$decemerge <- ifelse(!is.na(tab4$dmin) & tab4$Classe >= tab4$dmin, tab4$dec, tab4$defaut)
tab <- tab4 %>%
dplyr::select(espar, essence, diam, Classe, htot, hdec, decemerge, tarif, hauteur, houppier)
if (houppier == "N") {
tab <- tab %>%
dplyr::mutate(
L_VbftigCom = as.numeric(TarONF3v(tarif = tarif, entr1 = diam, entr2 = hauteur, details = FALSE)),
L_VHouppiers = L_VbftigCom * houppier / 100,
L_Vbftot7cm = L_VbftigCom + L_VHouppiers,
L_PHouppiers = houppier / 100
)
} else {
tab <- tab %>%
dplyr::mutate(
L_Vbftot7cm = as.numeric(TarONF3v(tarif = tarif, entr1 = diam, entr2 = hauteur, details = FALSE)),
L_VHouppiers = L_Vbftot7cm * houppier / 100,
L_VbftigCom = L_Vbftot7cm - L_VHouppiers,
L_PHouppiers = houppier / 100
)
}
tab <- cbind(tab, E_VbftigCom = TarEmerge(c130 = pi * tab$diam, htot = tab$htot, hdec = tab$hdec, espar = tab$espar, typevol = "tige", dec = tab$decemerge))
tab <- cbind(tab, E_Vbftot7cm = TarEmerge(c130 = pi * tab$diam, htot = tab$htot, hdec = tab$hdec, espar = tab$espar, typevol = "total", dec = 7)) %>%
dplyr::mutate(E_PHouppiers = E_Vbftot7cm / E_VbftigCom - 1) %>%
dplyr::filter(!is.na(E_Vbftot7cm)) %>%
dplyr::mutate(E_VHouppiers = E_Vbftot7cm - E_VbftigCom)
tab1 <- tab %>%
dplyr::mutate_if(is.numeric, funs(round(., 2)))
if (sum(c("essence", "diam", "htot", "hdec") %in% names(tab)) == 4) {
if (typvolemerge == "tige") {
tab <- tab %>%
dplyr::filter(diam >= classearbremin & diam <= classearbremax) %>%
dplyr::mutate(
numSchR = E_VbftigCom / 5 * 70000 / (diam - 5) / (diam - 10) - 8,
numSchL = E_VbftigCom / 5 * 90000 / diam / (diam - 5) - 8,
numAlg = E_VbftigCom * 28000000 / (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3) - 8
)
} else {
tab <- tab %>%
dplyr::filter(diam >= classearbremin & diam <= classearbremax) %>%
dplyr::mutate(
numSchR = E_Vbftot7cm / 5 * 70000 / (diam - 5) / (diam - 10) - 8,
numSchL = E_Vbftot7cm / 5 * 90000 / diam / (diam - 5) - 8,
numAlg = E_Vbftot7cm * 28000000 / (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3) - 8
)
}
res <- tab %>%
dplyr::group_by(essence) %>%
dplyr::summarise_at(vars(numSchR, numSchL, numAlg), funs(mean, var))
res[, 5:7] <- round(res[, 5:7]^0.5 / res[, 2:4], 2)
res[, 2:4] <- round(res[, 2:4], 2)
names(res) <- c("essence", "SchR", "SchL", "Alg", "SchRcv", "SchLcv", "Algcv")
tab <- tab %>%
dplyr::left_join(res, by = "essence") %>%
dplyr::mutate(
VSchR = ((round(SchR, 0) + 8) / 14000) * (diam - 5) * (diam - 10),
VSchL = ((round(SchL, 0) + 8) / 18000) * diam * (diam - 5),
VAlg = ((round(Alg, 0) + 8) / 28000000 * (310000 - 45200 * diam + 2390 * diam^2 - 2.9 * diam^3))
)
tab.p <- tab %>%
dplyr::select(essence, diam, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg")))
if (typvolemerge == "tige") {
tab.q <- reshape2::melt(tab, id.vars = c("diam", "essence"), measure.vars = c("L_VbftigCom", "E_VbftigCom", "VSchL", "VSchR", "VAlg"))
} else {
tab.q <- reshape2::melt(tab, id.vars = c("diam", "essence"), measure.vars = c("L_Vbftot7cm", "E_Vbftot7cm", "VSchL", "VSchR", "VAlg"))
}
names(tab.q) <- c("diam", "essence", "Type", "Vol")
###################################################################
# ça plante avec shinyproxy.io !!! bug dplyr memory mapped
# tab.r <- tab %>%
# dplyr::group_by(essence, Classe) %>%
# dplyr::summarise_at(vars(L_VbftigCom,L_VHouppiers,E_VbftigCom,E_VHouppiers), mean, na.rm=TRUE)
tab.r1 <- aggregate(L_VbftigCom ~ essence + Classe, tab, mean)
tab.r2 <- aggregate(L_VHouppiers ~ essence + Classe, tab, mean)
tab.r3 <- aggregate(E_VbftigCom ~ essence + Classe, tab, mean)
tab.r4 <- aggregate(E_VHouppiers ~ essence + Classe, tab, mean)
tab.r <- bind_cols(tab.r1, tab.r2[3], tab.r3[3], tab.r4[3]) %>%
arrange(essence, Classe)
###################################################################
tab.r <- reshape2::melt(tab.r, id.vars = c("Classe", "essence"), measure.vars = c("L_VbftigCom", "L_VHouppiers", "E_VbftigCom", "E_VHouppiers")) %>%
dplyr::mutate(Type = substr(variable, 1, 1), variable = substr(variable, 3, 12))
names(tab.r) <- c("Classe", "essence", "tarif", "Vol", "Type")
tab.r$tarif <- factor(tab.r$tarif, levels = c("VHouppiers", "VbftigCom"))
m1 <- aggregate(Num ~ essence + Type, tab.p, mean)
m1$var <- aggregate(Num ~ essence + Type, tab.p, var)[, 3]
m1$sd <- aggregate(Num ~ essence + Type, tab.p, sd)[, 3]
m1$cv <- round(m1$var^0.5 / m1$Num, 2)
m12 <- aggregate(Num ~ essence + Type, tab.p, median)
colnames(m12) <- c("essence", "Type", "Median")
m11 <- aggregate(Num ~ essence + Type, tab.p, FUN = function(x) quantile(x, probs = 0.25))
colnames(m11) <- c("essence", "Type", "Q1")
m13 <- aggregate(Num ~ essence + Type, tab.p, FUN = function(x) quantile(x, probs = 0.75))
colnames(m13) <- c("essence", "Type", "Q3")
m2 <- aggregate(diam ~ essence + Type, tab.p, max)
colnames(m2) <- c("essence", "Type", "DMax")
mm <- m1 %>%
dplyr::group_by(essence, Type) %>%
dplyr::summarise_all(mean)
mm <- merge(mm, m2)
mm <- merge(mm, m11)
mm <- merge(mm, m12)
mm <- merge(mm, m13)
mm <- mm %>%
dplyr::mutate(m_plus_sd = Num + 1.96 * sd) %>%
dplyr::mutate(m_moins_sd = Num - 1.96 * sd)
message("Plot...")
# EMERGE + LOCAL diagramme empile
p0 <- ggplot(tab.r, aes(x = Type, y = Vol, fill = Type, alpha = tarif)) +
scale_fill_manual(values = c("red", "darkgreen")) +
geom_bar(stat = "identity", position = "stack") +
facet_grid(essence ~ Classe)
# Boxplot
p1 <- ggplot(tab.p, aes(x = diam, y = Num)) + geom_point(alpha = 0.1) +
facet_grid(essence ~ Type) +
geom_boxplot(outlier.colour = "green", outlier.size = 2, notch = TRUE, alpha = 0.1, fill = "green") +
geom_label(data = mm, aes(x = DMax + 15, label = round(Num, 0), y = Num), col = "red") +
geom_label(data = mm, aes(x = DMax + 15, label = round(m_plus_sd, 0), y = m_plus_sd), col = "blue", alpha = 0.1) +
geom_label(data = mm, aes(x = DMax + 15, label = round(m_moins_sd, 0), 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_smooth(method = "lm") +
geom_point(data = mm, aes(x = DMax + 10), color = "red")
# EMERGE + LOCAL courbe lineaire
p2 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("darkgreen", "red"))
if (typvolemerge == "tige") {
p2 <- p2 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom", "L_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom", "L_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5)
} else {
p2 <- p2 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm", "L_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm", "L_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5)
}
# EMERGE + VSchR
p3 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p3 <- p3 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchR"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchR")
} else {
p3 <- p3 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchR"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchR")
}
# EMERGE + VSchL
p4 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p4 <- p4 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchL"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchL")
} else {
p4 <- p4 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VSchL"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VSchL")
}
# EMERGE + VAlg
p5 <- ggplot(tab.q, aes(x = diam, y = Vol, colour = Type, group = Type)) +
facet_grid(essence ~ .) + scale_color_manual(values = c("red", "darkgreen"))
if (typvolemerge == "tige") {
p5 <- p5 + geom_point(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VAlg"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_VbftigCom"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VAlg")
} else {
p5 <- p5 + geom_point(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(colour = Type), alpha = 0.1) +
geom_line(data = tab.q[tab.q$Type %in% c("VAlg"), ], aes(colour = Type), size = 1) +
geom_smooth(data = tab.q[tab.q$Type %in% c("E_Vbftot7cm"), ], aes(x = diam, y = Vol), span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
ggtitle("VAlg")
}
# Htot= f(E_Vbftig | Vbftot)
if (typvolemerge == "tige") {
tab.h <- tab %>%
dplyr::select(essence, diam, Classe, htot, E_VbftigCom, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam, -Classe, -htot, -E_VbftigCom) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg"))) %>%
dplyr::mutate(Classe = factor(Classe))
} else {
tab.h <- tab %>%
dplyr::select(essence, diam, Classe, htot, E_Vbftot7cm, numSchR, numSchL, numAlg) %>%
tidyr::gather(Type, Num, -essence, -diam, -Classe, -htot, -E_Vbftot7cm) %>%
dplyr::mutate(Type = factor(Type, levels = c("numSchR", "numSchL", "numAlg"))) %>%
dplyr::mutate(Classe = factor(Classe))
}
p <- c(0.10, 0.35, 0.65, 0.90)
p_names <- purrr::map_chr(p, ~paste0(.x*100))
p_fh <- purrr::map(p, ~partial(quantile, probs = .x, na.rm = TRUE)) %>%
purrr::set_names(nm = paste0("h",p_names))
p_fn <- purrr::map(p, ~partial(quantile, probs = .x, na.rm = TRUE)) %>%
purrr::set_names(nm = paste0("n",p_names))
tab.h1 <- tab.h %>%
group_by(essence, Classe, Type) %>%
summarize_at(vars(htot), funs(!!!p_fh))
tab.h2 <- tab.h %>%
group_by(essence, Classe, Type) %>%
summarize_at(vars(Num), funs(!!!p_fn))
tab.n <- merge(tab.h1, tab.h2)
tab.h3 <- tab.h %>%
group_by(essence, Type) %>%
summarize_at(vars(htot), funs(!!!p_fh)) %>%
mutate_if(is.numeric, round)
m3 <- aggregate(Num ~ essence + Type + Classe, tab.h, mean)
m3$var <- aggregate(Num ~ essence + Type + Classe, tab.h, var)[, 4]
m3$sd <- aggregate(Num ~ essence + Type + Classe, tab.h, sd)[, 4]
m3$cv <- round(m3$var^0.5 / m3$Num, 2)
if (typvolemerge == "tige") {
m4 <- aggregate(E_VbftigCom ~ essence + Type + Classe, tab.h, mean)
} else {
m4 <- aggregate(E_Vbftot7cm ~ essence + Type + Classe, tab.h, mean)
}
colnames(m4) <- c("essence", "Type", "Classe", "VMean")
m5 <- aggregate(htot ~ essence + Type + Classe, tab.h, mean)
colnames(m5) <- c("essence", "Type", "Classe", "HMean")
mmm <- m3 %>%
dplyr::group_by(essence, Type, Classe) %>%
dplyr::summarise_all(mean)
mmm <- merge(mmm, m4)
mmm <- merge(mmm, m5)
mmm <- mmm %>%
dplyr::mutate(m_plus_sd = Num + 1.96 * sd) %>%
dplyr::mutate(m_moins_sd = Num - 1.96 * sd)
# Schaeffer rapide
if (typvolemerge == "tige") {
p6 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchR") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p6 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchR"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchR") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
# Schaeffer lent
if (typvolemerge == "tige") {
p7 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchL") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p7 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numSchL"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numSchL") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
# Algan
if (typvolemerge == "tige") {
p8 <- ggplot(tab.h, aes(x = htot, y = E_VbftigCom, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numAlg") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
} else {
p8 <- ggplot(tab.h, aes(x = htot, y = E_Vbftot7cm, colour = Classe)) + geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(Num, 0), y = VMean, fill = Classe), col = "red") +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_plus_sd, 0), y = VMean + 0.5, fill = Classe), col = "blue", alpha = 0.5) +
geom_label(data = mmm[mmm$Type %in% c("numAlg"), ], aes(x = HMean, label = round(m_moins_sd, 0), y = VMean - 0.5, fill = Classe), col = "blue", alpha = 0.5) +
ggtitle("numAlg") +
scale_fill_discrete(guide = guide_legend(reverse = TRUE)) +
scale_color_discrete(guide = guide_legend(reverse = TRUE))
}
if (!is.null(echantillon)) {
# Calcul les centroides
centroids <- aggregate(cbind(echantillon$diam, echantillon$htot, echantillon$Num) ~ essence + Type, echantillon, mean)
f <- function(z) qt(0.025, df = length(z) - 1, lower.tail = F) * sd(z) / sqrt(length(z)) # 95% confidence
# f <- function(z)sd(z)/sqrt(length(z)) # function to calculate std.err
se <- aggregate(cbind(se.x = echantillon$diam, se.y = echantillon$htot, se.n = echantillon$Num) ~ essence + Type, echantillon, f)
centroids <- merge(centroids, se) # add std.err column to centroids
names(centroids) <- c("essence", "Type", "diam", "htot", "Num", "se.x", "se.y", "se.n")
}
# Numero de tarif numSchR
if (is.null(echantillon)) {
p9 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)), label = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchR") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p9 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)), label = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchR") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numSchR")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numSchR")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numSchR")), ]
p12 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numSchR")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numSchR") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
# Numero de tarif numSchL
if (is.null(echantillon)) {
p10 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchL") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p10 <- ggplot(tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numSchL") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numSchL")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numSchL")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numSchL")), ]
p13 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numSchL")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numSchL") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
# Numero de tarif numAlg
if (is.null(echantillon)) {
p11 <- ggplot(tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numAlg") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1))
} else {
p11 <- ggplot(tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot, color = factor(round(Num, 0)), fill = factor(round(Num, 0)))) +
geom_point() +
facet_grid(essence ~ .) + geom_smooth(method = "lm", linetype = 2, size = 0.5) +
ggtitle("numAlg") +
scale_color_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
scale_fill_discrete("Num", guide = guide_legend(reverse = TRUE, ncol = 1)) +
geom_point(data = echantillon[which(echantillon$Type %in% c("numAlg")), ], aes(x = diam, y = htot), shape = 3, size = 4, col = "black") +
geom_label(data = centroids[which(centroids$Type %in% c("numAlg")), ], aes(x = diam, y = htot, label = round(Num, 0)), col = "black") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h10,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h35,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h65,
color = "black",
fill="white") +
geom_label(data = tab.h3[which(tab.h3$Type %in% c("numSchR")), ],
x = 40,
y = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
label = tab.h3[which(tab.h3$Type %in% c("numSchR")), ]$h90,
color = "black",
fill="white")
}
tabo <- tab.n[which(tab.n$Type %in% c("numAlg")), ]
p14 <- ggplot() +
geom_point(data = tab.h[which(tab.h$Type %in% c("numAlg")), ], aes(x = diam, y = htot)) +
facet_grid(essence ~ .) +
ggtitle("numAlg") +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe)[tabo$Classe]), y = h10),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h10, label = floor(n10), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h35, label = floor(n35))) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5, col = "red") +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h65, label = floor(n65), col = "red")) +
geom_smooth(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90),
span = 1, method = "loess", se = TRUE, alpha = 0.2, size = 0.5) +
geom_label_repel(data = tabo,
aes(x = as.numeric(levels(tabo$Classe))[tabo$Classe], y = h90, label = floor(n90))) +
theme(legend.position="none")
if (enreg) {
readr::write_excel_csv(tab, paste0(dirname(fichier), "/numTarifSchaeffer.csv"))
message(paste0("Result is saved in: ", paste0(dirname(fichier), "/numTarifSchaeffer.csv")))
}
if (mappoint) {
out <- list(tab1, res, tab.p, mm, tab.h, tab.n, p1, p2, p3, p4, p5, p0, placettes, zone, p6, p7, p8, p9, p10, p11, p12, p13, p14)
names(out) <- c(
"Tableau1", "Tableau2", "Tableau3", "Tableau4", "Tableau5", "Tableau6",
"Graphe1", "Graphe2", "Graphe3", "Graphe4", "Graphe5", "Graphe6",
"Placettes", "Zone", "Graphe7", "Graphe8", "Graphe9", "Graphe10",
"Graphe11", "Graphe12", "Graphe13", "Graphe14", "Graphe15"
)
} else {
out <- list(tab1, res, tab.p, mm, tab.h, tab.n, p1, p2, p3, p4, p5, p0, p6, p7, p8, p9, p10, p11, p12, p13, p14)
names(out) <- c(
"Tableau1", "Tableau2", "Tableau3", "Tableau4", "Tableau5", "Tableau6",
"Graphe1", "Graphe2", "Graphe3", "Graphe4", "Graphe5", "Graphe6",
"Graphe7", "Graphe8", "Graphe9", "Graphe10", "Graphe11", "Graphe12",
"Graphe13", "Graphe14", "Graphe15"
)
}
message("Calculation realized!")
return(out)
} else {
message("The file must contain the columns 'essence', 'diam', 'htot' and 'hdec' (even empty)!")
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.