In thomasgaquiere/Maria: Simulate Different Forms of Logging on a Forest
# si on avait une ancienne version de withr pour rstan
remove.packages("withr")
install.packages("withr")
library(EcoFoG)
library(Maria)
library(readr)
library(dplyr)
library(tibble)
library(tidyr)
library(schoolmath)
library(ggplot2)
library(reshape)
library(sf)
library(nngeo)
library(raster)
library(topoDistance)
library(BIOMASS)
Paracou6_2016 <- Paracou6_2016 %>%
dplyr::select(-TreeFieldNum,-Project, -idVern, -CommercialSp, -CensusDate)
usethis::use_data(Paracou6_2016, overwrite = TRUE)
ParamCrownDiameterAllometry <- read_delim("D:/VSC Kourou/DATA/ParamCrownDiameterAllometry.csv",
";", escape_double = FALSE, col_types = cols(alpha = col_number(),
beta = col_number()), locale = locale(decimal_mark = ","),
trim_ws = TRUE) %>%
separate(ScientificName, sep = "_", into = c("Genus", "Species"), remove = F) %>%
mutate(Family = ifelse(Taxo == "fam", Genus, NA)) %>%
mutate(Genus = ifelse(is.na(Family), Genus, NA))
usethis::use_data(ParamCrownDiameterAllometry)
SpeciesCriteria <- read_csv("D:/VSC Kourou/DATA/SpeciesCriteria.csv",
col_types = cols(Commercial = col_factor(levels = c("0",
"1", "2")), MinFD = col_number(),
UpMinFD = col_number(), MaxFD = col_number()),
locale = locale(encoding = "Latin1")) %>%
separate(ScientificName, sep = " ", into = c("Genus", "Species"), remove = T) %>%
mutate(Aggregative = ifelse(VernName == "Angélique" | VernName == "Gonfolo rose",
TRUE, FALSE))
# unite(Genus, Species, col = "ScientificName", sep = "_", remove = F)
SpeciesCriteria$VernName <- tolower(SpeciesCriteria$VernName)
# usethis::use_data(SpeciesCriteria, overwrite = TRUE)
ScenariosTable <- read_csv("D:/VSC Kourou/DATA/ScenariosTable.csv",
col_types = cols(Winching = col_factor(levels = c("0",
"1", "2")), DirectionalFelling = col_factor(levels = c("0",
"1", "2"))), trim_ws = TRUE)
# usethis::use_data(ScenariosTable,overwrite = TRUE)
Data loading
data(Paracou6_2016) #data inventory in the global environment
data(ForestZoneVolumeParametersTable) # The volume parameters data in the global environment
# load(system.file("extdata", "BrokenParacou6_2016.rda", package = "Maria")) # the test inventory in the global environment
data(ParamCrownDiameterAllometry)
data(SpeciesCriteria)
lapply(Paracou6_2016, class) #give the class of all the data variables
inventorycheckformat
INPUT: inventory = Inventaires botaniques Guyafor/ONF (data frame)
OUTPUT: inventory = Inventaires botaniques type “Guyafor”(data frame)
vérification du type d’inventaire :
Test présence des variables.
Test class des variables.
Si les tests précédents ont tous reçu =TRUE, STOP de la fonction inventorycheckformat
Sinon mettre en forme “inventaire Guyafor” (nom des var, leurs classes, unités (mettre des limites avec range),...)
re-vérification du type d’inventaire :
Test présence des variables
Test class des variables
Si les tests précédent ont tous reçu =TRUE, STOP de la fonction inventorycheckformat
Sinon Message : “Votre format de données renseignés à l’argument “inventory”, n’est pas compatible. Voir la vignette du package pour y trouver le format requis.”
'pouet' %in% names(Paracou6_2016)
if (!is.null(Paracou6_2016$pouet))
stop("pouet variable should be an integer.")
# Variables presence check
if(!("Plot" %in% names(Paracou6_2016)))
stop("Plot variable is not found")
if(!("CensusYear" %in% names(Paracou6_2016)))
stop("CensusYear variable is not found")
if(!("idTree" %in% names(Paracou6_2016)))
stop("idTree variable is not found")
if(!("Family" %in% names(Paracou6_2016)))
stop("Family variable is not found")
if(!("Genus" %in% names(Paracou6_2016)))
stop("Genus variable is not found")
if(!("Species" %in% names(Paracou6_2016)))
stop("Species variable is not found")
if(!("CircCorr" %in% names(Paracou6_2016)))
stop("CircCorr variable is not found")
if(!("CodeAlive" %in% names(Paracou6_2016)))
stop("CodeAlive variable is not found")
if(!("CommercialSp" %in% names(Paracou6_2016)))
stop("CommercialSp variable is not found")
if(!("UTMZone" %in% names(Paracou6_2016)))
stop("UTMZone variable is not found")
if(!('Lat' %in% names(Paracou6_2016)))
stop("Lat variable is not found")
if(!('Lon' %in% names(Paracou6_2016)))
stop("Lon variable is not found")
if(!('VernName' %in% names(Paracou6_2016)))
stop("VernName variable is not found")
if(!('Xfield' %in% names(Paracou6_2016)))
stop("Xfield variable is not found")
if(!('Yfield' %in% names(Paracou6_2016)))
stop("Yfield variable is not found")
if(!('Xutm' %in% names(Paracou6_2016)))
stop("Xutm variable is not found")
if(!('Yutm' %in% names(Paracou6_2016)))
stop("Yutm variable is not found")
# Variables class check
if (!inherits(Paracou6_2016$idTree, "integer") && !is.null(idTree)) #majuscule ?
stop("idTree variable should be an integer.")
if (!inherits(Paracou6_2016$Plot, "character") && !is.null(Plot))
stop("Plot variable should be a character.")
if (!inherits(Paracou6_2016$Xfield, "numeric") && !is.null(Xfield))
stop("Xfield variable should be numeric.")
if (!inherits(Paracou6_2016$Yfield, "numeric") && !is.null(Yfield))
stop("Yfield variable should be numeric.")
if (!inherits(Paracou6_2016$Xutm, "numeric") && !is.null(Xutm))
stop("Xutm variable should be numeric.")
if (!inherits(Paracou6_2016$Yutm, "numeric") && !is.null(Yutm))
stop("Yutm variable should be numeric.")
if (!inherits(Paracou6_2016$UTMZone, "integer") && !is.null(UTMZone))
stop("UTMZone variable should be an integer.")
if (!inherits(Paracou6_2016$Lat, "numeric") && !is.null(Lat))
stop("Lat variable should be numeric.")
if (!inherits(Paracou6_2016$Lon, "numeric") && !is.null(Lon))
stop("Lon variable should be numeric.")
if (!inherits(Paracou6_2016$Family, "character") && !is.null(Family))
stop("Family variable should be a character.")
if (!inherits(Paracou6_2016$Genus, "character") && !is.null(Genus))
stop("Genus variable should be a character.")
if (!inherits(Paracou6_2016$Species, "character") && !is.null(Species))
stop("Species variable should be a character.")
if (!inherits(Paracou6_2016$VernName, "character") && !is.null(VernName))
stop("VernName variable should be a character.")
if (!inherits(Paracou6_2016$CommercialSp, "logical") && !is.null(CommercialSp))
stop("CommercialSp variable should be logical.")
if (!inherits(Paracou6_2016$CensusYear, "integer") && !is.null(CensusYear))
stop("CensusYear variable should be an integer.")
if (!inherits(Paracou6_2016$CodeAlive, "logical") && !is.null(CodeAlive))
stop("CodeAlive variable should be logical.")
if (!inherits(Paracou6_2016$CircCorr, "numeric") && !is.null(CircCorr)) #DBH
stop("CircCorr variable should be numeric.")
#
# if (!inherits(Paracou6_2016$Forest, "character"))
# stop("Forest variable should be a character.") #pas besoin
# if (!inherits(Paracou6_2016$idVern, "integer") && !is.null(idVern))
# stop("idVern variable should be an integer.") #pas besoin
#
#
#
# if (!inherits(Paracou6_2016$PlotArea, "numeric") && !is.null(PlotArea))
# stop("PlotArea variable should be numeric.") #pas besoin
#
# if (!inherits(Paracou6_2016$SubPlot, "integer") && !is.null(SubPlot))
# stop("SubPlot variable should be an integer.") #pas besoin
#
# if (!inherits(Paracou6_2016$TreeFieldNum, "numeric") && !is.null(TreeFieldNum))
# stop("TreeFieldNum variable should be numeric.") #pas besoin
#
# if (!inherits(Paracou6_2016$Project, "character") && !is.null(Project))
# stop("Project variable should be a character.") #pas besoin
#
# if (!inherits(Paracou6_2016$Protocole, "character") && !is.null(Protocole))
# stop("Protocole variable should be a character.") #pas besoin
if(!("Plot" %in% names(Paracou6_2016)))
stop("Plot variable is not found")
if(!("CensusYear" %in% names(Paracou6_2016)))
stop("CensusYear variable is not found")
if(!("idTree" %in% names(Paracou6_2016)))
stop("idTree variable is not found")
if(!("Family" %in% names(Paracou6_2016)))
stop("Family variable is not found")
if(!("Genus" %in% names(Paracou6_2016)))
stop("Genus variable is not found")
if(!("Species" %in% names(Paracou6_2016)))
stop("Species variable is not found")
if(!("CircCorr" %in% names(Paracou6_2016)))
stop("CircCorr variable is not found")
if(!("CodeAlive" %in% names(Paracou6_2016)))
stop("CodeAlive variable is not found")
if(!("CommercialSp" %in% names(Paracou6_2016)))
stop("CommercialSp variable is not found")
if(!("UTMZone" %in% names(Paracou6_2016)))
stop("UTMZone variable is not found")
if(!('Lat' %in% names(Paracou6_2016)))
stop("Lat variable is not found")
if(!('Lon' %in% names(Paracou6_2016)))
stop("Lon variable is not found")
if(!('VernName' %in% names(Paracou6_2016)))
stop("VernName variable is not found")
if(!('Xfield' %in% names(Paracou6_2016)))
stop("Xfield variable is not found")
if(!('Yfield' %in% names(Paracou6_2016)))
stop("Yfield variable is not found")
if(!('Xutm' %in% names(Paracou6_2016)))
stop("Xutm variable is not found")
if(!('Yutm' %in% names(Paracou6_2016)))
stop("Yutm variable is not found")
if (!inherits(Paracou6_2016$idTree, "integer.")) #majuscule ?
stop("idTree variable should be an integer.")
if (!inherits(Paracou6_2016$Plot, "character"))
stop("Plot variable should be a character.")
if (!inherits(Paracou6_2016$Xfield, "numeric"))
stop("Xfield variable should be numeric.")
if (!inherits(Paracou6_2016$Yfield, "numeric"))
stop("Yfield variable should be numeric.")
if (!inherits(Paracou6_2016$Xutm, "numeric"))
stop("Xutm variable should be numeric.")
if (!inherits(Paracou6_2016$Yutm, "numeric"))
stop("Yutm variable should be numeric.")
if (!inherits(Paracou6_2016$UTMZone, "integer"))
stop("UTMZone variable should be an integer.")
if (!inherits(Paracou6_2016$Lat, "numeric"))
stop("Lat variable should be numeric.")
if (!inherits(Paracou6_2016$Lon, "numeric"))
stop("Lon variable should be numeric.")
if (!inherits(Paracou6_2016$Family, "character"))
stop("Family variable should be a character.")
if (!inherits(Paracou6_2016$Genus, "character"))
stop("Genus variable should be a character.")
if (!inherits(Paracou6_2016$Species, "character"))
stop("Species variable should be a character.")
if (!inherits(Paracou6_2016$VernName, "character"))
stop("VernName variable should be a character.")
if (!inherits(Paracou6_2016$CommercialSp, "logical"))
stop("CommercialSp variable should be logical.")
if (!inherits(Paracou6_2016$CensusYear, "integer"))
stop("CensusYear variable should be an integer.")
if (!inherits(Paracou6_2016$CodeAlive, "logical"))
stop("CodeAlive variable should be logical.")
if (!inherits(Paracou6_2016$CircCorr, "numeric")) #DBH
stop("CircCorr variable should be numeric.")
cleaninventory
INPUT: inventory = Inventaires botaniques type “Guyafor” (Paracou6_2016)(data frame)
OUTPUT: inventory = Inventaires botaniques exploitable (Paracou6_2016Clean) (data frame)
Ne prendre que les arbres vivants : lignes ayant la valeur “TRUE” à la colonne “CodeAlive” de inventory.
Check si les valeurs de ID sont uniques. Si FALSE message d’erreur “Les identifiants des arbres ne sont pas uniques”
Enlever les arbres hors de la parcelle s’il y en a : coordonnées de l’arbre non comprises dans les coordonnées de la parcelle.
Ne prendre que les lignes pour lesquelles “DBH”> 10.
duplicated(Paracou6_2016)
unique(Paracou6_2016) #same rows number
#duplicated() gives "FALSE" when the row is unique
# any() : at least one of the values true?
any(duplicated(Paracou6_2016)) # at least one row is duplicated?
all(Paracou6_2016$Plot==Paracou6_2016$Plot[1]) #all the Plot values are equal?
if (!("DBH" %in% names(Paracou6_2016))) {add_column(Paracou6_2016, DBH = NA) #if DBH (cm) doesn't exist create it
Paracou6_2016$DBH = Paracou6_2016$CircCorr/pi} # and compute it
InventoryClean <- Paracou6_2016 %>%
filter(CodeAlive == "TRUE") %>% #only alive trees
filter(DBH >= 10) # DBH >= 10
if (any(duplicated(Paracou6_2016)==TRUE))
stop ("Tree identifiers (idTree) are not unique.") # stop function if the tree identifiers (idTree) are not unique
if (all(Paracou6_2016$Plot==Paracou6_2016$Plot[1])==FALSE) #all the Plot values are equal?
stop ("Input inventory concerns different plots (Plot). Our function simulates logging at the plot level.")
# Remove out-plot trees (A faire)
addtreedim
Fonction d’ajout de variables pour l’INPUT “inventory” :
INPUT: inventory = InventoryClean (data frame)
OUTPUT: inventory = Inventory_addtreedim (data frame)
calculer pour tous les arbres :
+ TreeBiomass : biomasse de tous les arbres (BIOMASS package)
+ TreeHeight : hauteur de l’arbre (H) (BIOMASS package)
+ TreeHarvestableVolume : volume exploitable à partir du tarif de cubage (= a + bDBH2) que représente chaque arbre (a et b dépendent de la localisation)
+ TrunkHeight : hauteur de fût (TreeHarvestableVolume = π(DBH/2)² x TrunkHeight)
+ CrownDiameter : diamètre de couronne (CD) (ln(D) = 𝜶+ 𝜷 ln(HCD) + 𝜺 (allométries de Mélaine)
+ CrownHeight : TreeHeight - TrunkHeight
library(BIOMASS)
library(knitr)
TreeheightGuyafordata <- read_delim("D:/VSC Kourou/DATA/TreeheightGuyafordata.csv",
";", escape_double = FALSE, col_types = cols(Circ = col_number(),
CircCorr = col_number(), Hauteur = col_number()),
locale = locale(encoding = "Latin1"),
trim_ws = TRUE) %>%
mutate(HeightHow = recode(HeightHow, 'Visual estimate' = 'Visual Estimate')) %>%
# filter(!(HeightHow == "Visual Estimate" | HeightHow == "Visual estimate")) %>%
unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>%
mutate(DBH = ifelse(is.na(CircCorr), Circ/pi, CircCorr/pi)) # add DBH column
any(is.na(TreeheightGuyafordata$DBH)) # all the DBH have been computed?
SineTelemeter <- TreeheightGuyafordata %>% #data to fit model
filter(HeightHow == "Sine Method (Larjavaara & Muller-Landau 2013) with Trimble LaserAce 1000 Rangefinder"| HeightHow == "Telemeter")
ModelsresultGuiana <-modelHD(
D = SineTelemeter$DBH,
H = SineTelemeter$Hauteur,
useWeight = TRUE
)
kable(ModelsresultGuiana) #log2 has the lowest RSE (4.699261 with all data, 4.316927 without "Visual E/estimate", 3.059860 with just "Sine Method", 3.809011 with sine and telemeter)
HDmodelGuiana <- modelHD(
D = SineTelemeter$DBH,
H = SineTelemeter$Hauteur,
method = "log2", # Compute the H-D model with the lowest RSE.
useWeight = TRUE
)
TreeHeightGuianaEstimation <- retrieveH(
D = TreeheightGuyafordata$DBH,
model = HDmodelGuiana
)
#Put these tree heights estimations in the inventory
TreeheightGuyafordataEstim <- TreeheightGuyafordata %>%
mutate(HauteurEstimee = TreeHeightGuianaEstimation$H)
plot(TreeheightGuyafordataEstim$DBH, TreeheightGuyafordataEstim$HauteurEstimee)
# obs/estim regression:
TreeheightGuyafordataEstim %>%
ggplot() +
aes(x = Hauteur, y = HauteurEstimee) +
geom_point(shape = "circle", size = 0.5, colour = "#440154") +
geom_smooth(span = 0.75) +
theme_minimal()
# We are well between 10 and 35m height
library(BIOMASS)
library(knitr)
# Just with Nouragues data (measurement: laser rangefinder):
# Fit first different height-diameter models
ModelsresultNouraguesdata <- modelHD(
D = NouraguesHD$D,
H = NouraguesHD$H,
useWeight = TRUE # model weights will be (D^2)*H (weights are proportional to tree volume, so that larger trees have a stronger influence during the construction of the model)
)
# choose the best model
kable(ModelsresultNouraguesdata) #log2 has the lowest RSE (4.222718) (Root Squared Error)
HDmodelNouraguesdata <- modelHD(
D = NouraguesHD$D,
H = NouraguesHD$H,
method = "log2", # Compute the local H-D model with the lowest RSE.
useWeight = TRUE
)
# Applicate the model to retrieve tree heights ---- HERE
TreeHeightEstimation <- retrieveH(
D = Paracou6_2016$DBH,
model = HDmodelNouraguesdata
)
## Put these tree heights estimations in the inventory
Paracou6_2016TreeHeight <- Paracou6_2016 %>%
mutate(TreeHeight = TreeHeightEstimation$H)
plot(Paracou6_2016TreeHeight$DBH, Paracou6_2016TreeHeight$TreeHeight)#amazing !! It's a log2 form
# Different guyanese data:
TreeheightGuyafordata <- read_delim("D:/VSC Kourou/DATA/TreeheightGuyafordata.csv",
";", escape_double = FALSE, col_types = cols(Circ = col_number(),
CircCorr = col_number(), Hauteur = col_number()),
locale = locale(encoding = "Latin1"),
trim_ws = TRUE) %>%
# filter(!(HeightHow == "Visual Estimate" | HeightHow == "Visual estimate"))
filter(HeightHow == "Sine Method (Larjavaara & Muller-Landau 2013) with Trimble LaserAce 1000 Rangefinder")
if (!("DBH" %in% names(TreeheightGuyafordata))) {add_column(TreeheightGuyafordata, DBH = NA) #if DBH (cm) doesn't exist create it
TreeheightGuyafordata$DBH = TreeheightGuyafordata$CircCorr/pi} # and compute it
unique(TreeheightGuyafordata$HeightHow)
ModelsresultGuiana <- modelHD(
D = TreeheightGuyafordata$DBH,
H = TreeheightGuyafordata$Hauteur,
useWeight = TRUE
)
kable(ModelsresultGuiana) #log2 has the lowest RSE (4.749107 with all data, 4.316927 without "Visual E/estimate", 3.059860 with just "Sine Method")
#Don't forget, you have the model form, and its parameters, and there are different according to the data.
HDmodelGuiana <- modelHD(
D = TreeheightGuyafordata$DBH,
H = TreeheightGuyafordata$Hauteur,
method = "log2", # Compute the local H-D model with the lowest RSE.
useWeight = TRUE
)
TreeHeightGuianaEstimation <- retrieveH(
D = Paracou6_2016$DBH,
model = HDmodelGuiana
)
## Put these tree heights estimations in the inventory
Paracou6_2016TreeHeight <- Paracou6_2016 %>%
mutate(TreeHeight = TreeHeightGuianaEstimation$H)
plot(Paracou6_2016TreeHeight$DBH, Paracou6_2016TreeHeight$TreeHeight)#amazing !! It's a log2 form
# Question : ya t'il des arbres mesurés par plsrs méthodes la même année ?
TreeheightMethods <- TreeheightGuyafordata %>%
select(Forest, Plot, SubPlot, TreeFieldNum, Genus, Species, CensusYear) %>%
unique()
# TreeHarvestableVolume = volume exploitable à partir du tarif de cubage (= a + b*DBH2) que représente chaque arbre (a et b dépendent de la localisation)
ForestZoneVolumeParametersTable <- data.frame (Forest = c("Acarouany","BAFOG","Kaw","Laussat","Montagne Plomb","Montagne Tortue","Nouragues","Organabo","Paracou","Régina St Georges","Risquetout","Tibourou"), # Guyafor forests
Zone = c("East", "West", "Central", "FrenchGuiana")) %>% # Volume formules zones
mutate(Zone = ifelse(Forest == "Acarouany"| Forest =="BAFOG" |Forest =="Laussat", "West", NA)) %>% # West zone
mutate(Zone = ifelse(Forest == "Montagne Plomb"| Forest =="Organabo"|Forest =="Paracou"|Forest =="Risquetout", "Central", Zone)) %>% # Central zone
mutate(Zone = ifelse(Forest == "Kaw"| Forest =="Montagne Tortue"|Forest =="Nouragues"|Forest =="Régina St Georges"|Forest =="Tibourou", "East", Zone)) %>% # East zone
left_join(VolumeParametersTable)
Paracou6_2016_volume <- Paracou6_2016 %>%
left_join(ForestZoneVolumeParametersTable) %>%
mutate(TreeHarvestableVolume = aCoef + bCoef * CircCorr^2) # compute the tree harvestable volume
# West = c("Acarouany", "BAFOG", "Laussat")
# Central = c("Montagne Plomb", "Organabo","Paracou", "Risquetout")
# East = c("Kaw", "Montagne Tortue","Nouragues", "Régina St Georges", "Tibourou")
# Guyane = le reste par défaut
# et possibilité à l'utilisateur de mettre ses propres tarifs
mutate(TrunkHeight = TreeHarvestableVolume/(π(DBH/2)^2)) # compute the trunk height (cylinderVolume= π(DBH/2)² x H)
mutate(CrownHeight = TreeHeight - TrunkHeight) # compute the crown height
spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale
filter(Taxo == "sp") %>%
select(-Family)
genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale
filter(Taxo == "gen") %>%
select(-Family, -ScientificName, -Species)
famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale
filter(Taxo == "fam") %>%
select(-Genus, -ScientificName, -Species)
Paracou6_2016 %>%
sample_n(100) %>% # short data
select(idTree, Family, Genus, Species) %>% #only useful columns
left_join(select(spParamCrownDiameter, -Taxo, -ScientificName, -beta), by = c("Genus","Species")) %>%
left_join(select(genParamCrownDiameter, -Taxo, -beta), by = c("Genus"), suffix = c(".species", ".genus")) %>%
left_join(select(famParamCrownDiameter, -Taxo, -beta) %>%
rename(alpha.family = alpha), by = c("Family")) %>%
mutate(alpha = alpha.species) %>% #mutate(alpha = alpha.species,beta = beta.species)
mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>%
mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>%
mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>%
select(-alpha.species, -alpha.genus, -alpha.family)
spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale
filter(Taxo == "sp") %>%
select(-Family)
genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale
filter(Taxo == "gen") %>%
select(-Family, -ScientificName, -Species)
famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale
filter(Taxo == "fam") %>%
select(-Genus, -ScientificName, -Species)
error <- rnorm(2000, mean = 0, sd = 0.0295966977) #Mélaine's model error distribution (hist(error))
Paracou6_2016_crown <- Paracou6_2016 %>%
unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>%
left_join(spParamCrownDiameter, by = c("ScientificName","Genus","Species")) %>% #Add parameters at species scale
left_join(genParamCrownDiameter, by = "Genus", suffix = c(".species", ".genus")) %>% #Add parameters at genus scale
left_join(famParamCrownDiameter, by = "Family") %>% #Add parameters at family scale
dplyr::rename(Taxo.family = Taxo, alpha.family = alpha, beta.family = beta) %>%
mutate(alpha = alpha.species, beta = beta.species, Taxo = Taxo.species) %>% # create the future & unique alpha,beta,taxo columns
mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>% #if species parameters are absent, take genus parameters
mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>% #if species&genus parameters are absent, take family parameters
mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-alpha.species, -alpha.genus, -alpha.family) %>% #remove obsolete columns
mutate(beta = ifelse(is.na(beta), beta.genus, beta)) %>% #if species parameters are absent, take genus parameters
mutate(beta = ifelse(is.na(beta), beta.family, beta)) %>% #if species&genus parameters are absent, take family parameters
mutate(beta = ifelse(is.na(beta), mean(beta, na.rm =T), beta)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-beta.species, -beta.genus, -beta.family) %>% #remove obsolete columns
mutate(Taxo = ifelse(is.na(Taxo), Taxo.genus, Taxo)) %>% #if species parameters are absent, take genus parameters
mutate(Taxo = ifelse(is.na(Taxo), Taxo.family, Taxo)) %>% #if species&genus parameters are absent, take family parameters
mutate(Taxo = ifelse(is.na(Taxo), "mean", Taxo)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-Taxo.species, -Taxo.genus, -Taxo.family) %>% #remove obsolete columns
mutate(CrownDiameter = exp((
(log(DBH)-alpha-error)/beta)
)/TreeHeight) # compute the crown diameter (CD) (ln(D) = 𝜶+ 𝜷 ln(H*CD) + 𝜺, with 𝜺~N(0,σ^2) and meanσ^2 = 0.0295966977. (Mélaine's allometries))
# Informations:
# How many individuals have find their Crown diameter parameters at species, genus, family or no scale:
sum(Paracou6_2016_crown$Taxo == "fam") #442 individuals with family scale parameters
sum(Paracou6_2016_crown$Taxo == "gen") #742 individuals with genus scale parameters
sum(Paracou6_2016_crown$Taxo == "sp") #2364 individuals with family scale parameters
sum(Paracou6_2016_crown$Taxo == "mean") #72 individuals no recognize at any scale, so with plot mean parameters
unique(Paracou6_2016_crown[(Paracou6_2016_crown$Taxo == "mean"),]$ScientificName) #which species are no recognized
length(unique(Paracou6_2016_crown$ScientificName)) #How many species are in the plot
# allometry parameters data preparation:
spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale
filter(Taxo == "sp") %>%
select(-Family)
genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale
filter(Taxo == "gen") %>%
select(-Family, -ScientificName, -Species)
famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale
filter(Taxo == "fam") %>%
select(-Genus, -ScientificName, -Species)
error <- rnorm(2000, mean = 0, sd = 0.0295966977) #Mélaine's model error distribution (hist(error))
# Variables computation:
InventoryAddtreedim <- InventoryClean %>%
# TreeHarvestableVolume (m^3)
left_join(ForestZoneVolumeParametersTable) %>%
mutate(TreeHarvestableVolume = aCoef + bCoef * (DBH/100)^2) %>% # compute the tree harvestable volume
# TrunkHeight & CrownHeight (m)
mutate(TrunkHeight = TreeHarvestableVolume/(pi*(((DBH/100)/2)^2))) %>% # compute the trunk height (cylinderVolume= π(DBH/2)² x H) DBH in cm, in m in the formule.
# mutate(CrownHeight = TreeHeight - TrunkHeight) %>% # compute the crown height TREEHEIGHT EN ATTENTE
# CrownDiameter (m)
# Add crown diameter parameters at the inventory:
unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>%
left_join(spParamCrownDiameter, by = c("ScientificName","Genus","Species")) %>% #Add parameters at species scale
left_join(genParamCrownDiameter, by = "Genus", suffix = c(".species", ".genus")) %>% #Add parameters at genus scale
left_join(famParamCrownDiameter, by = "Family") %>% #Add parameters at family scale
rename(Taxo.family = Taxo, alpha.family = alpha, beta.family = beta) %>%
mutate(alpha = alpha.species, beta = beta.species, Taxo = Taxo.species) %>% # create the future & unique alpha,beta,taxo columns
mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>% #if species parameters are absent, take genus parameters
mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>% #if species&genus parameters are absent, take family parameters
mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-alpha.species, -alpha.genus, -alpha.family) %>% #remove obsolete columns
mutate(beta = ifelse(is.na(beta), beta.genus, beta)) %>% #if species parameters are absent, take genus parameters
mutate(beta = ifelse(is.na(beta), beta.family, beta)) %>% #if species&genus parameters are absent, take family parameters
mutate(beta = ifelse(is.na(beta), mean(beta, na.rm =T), beta)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-beta.species, -beta.genus, -beta.family) %>% #remove obsolete columns
mutate(Taxo = ifelse(is.na(Taxo), Taxo.genus, Taxo)) %>% #if species parameters are absent, take genus parameters
mutate(Taxo = ifelse(is.na(Taxo), Taxo.family, Taxo)) %>% #if species&genus parameters are absent, take family parameters
mutate(Taxo = ifelse(is.na(Taxo), "mean", Taxo)) %>% #if species,genus&family parameters are absent, take parameters mean
select(-Taxo.species, -Taxo.genus, -Taxo.family) #remove obsolete columns
# mutate(CrownDiameter = exp((
# (log(DBH)-alpha-error)/beta)
# )/TreeHeight) # compute the crown diameter (CD) (ln(D) = 𝜶+ 𝜷 ln(H*CD) + 𝜺, with 𝜺~N(0,σ^2) and meanσ^2 = 0.0295966977. (Mélaine's allometries)) # Mettre finalement l'objet créé dans le zzz.R
# !!!!!!!TREEHEIGHT EN ATTENTE!!!!!!
scenariosparameters
pas de mode "manual",
si argument = NULL, prendre la valeur du scénario correspondante
if (is.null(winching)){
if (Type == "RIL1") winching = "0"
if (Type == "RIL2broken") winching = "0"
if (Type == "RIL2") winching = "1"
if (Type == "RIL3") winching = "2"
if (Type == "RIL3fuel") winching = "2"
if (Type == "RIL3fuelhollow") winching = "2"
}
if (is.null(directionalfelling)){
if (Type == "RIL1") directionalfelling = "0"
if (Type == "RIL2broken") directionalfelling = "1"
if (Type == "RIL2") directionalfelling = "1"
if (Type == "RIL3") directionalfelling = "2"
if (Type == "RIL3fuel") directionalfelling = "2"
if (Type == "RIL3fuelhollow") directionalfelling = "2"
}
if (is.null(objective)){
if (Type == "RIL1") objective = 25
if (Type == "RIL2broken") objective = 25
if (Type == "RIL2") objective = 25
if (Type == "RIL3") objective = 30
if (Type == "RIL3fuel") objective = 30
if (Type == "RIL3fuelhollow") objective = 30
}
if (is.null(diversification)){
if (Type == "RIL1") diversification = FALSE
if (Type == "RIL2broken") diversification = FALSE
if (Type == "RIL2") diversification = FALSE
if (Type == "RIL3") diversification = TRUE
if (Type == "RIL3fuel") diversification = TRUE
if (Type == "RIL3fuelhollow") diversification = TRUE
}
# Type SpatialDataType Winching DirectionalFelling Objective Diversification
# RIL1 SRTM 0 0 20-25 FALSE
# RIL2broken LIDAR 0 1 20-25 FALSE
# RIL2 LIDAR 1 1 20-25 FALSE
# RIL3 LIDAR 2 2 25-30 TRUE
# RIL3fuel LIDAR 2 2 25-30 TRUE
# RIL3fuelhollow LIDAR 2 2 25-30 TRUE
thomasgaquiere/Maria documentation built on Dec. 24, 2021, 1:20 a.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.
# si on avait une ancienne version de withr pour rstan remove.packages("withr") install.packages("withr")
library(EcoFoG) library(Maria) library(readr) library(dplyr) library(tibble) library(tidyr) library(schoolmath) library(ggplot2) library(reshape) library(sf) library(nngeo) library(raster) library(topoDistance) library(BIOMASS)
Paracou6_2016 <- Paracou6_2016 %>% dplyr::select(-TreeFieldNum,-Project, -idVern, -CommercialSp, -CensusDate) usethis::use_data(Paracou6_2016, overwrite = TRUE)
ParamCrownDiameterAllometry <- read_delim("D:/VSC Kourou/DATA/ParamCrownDiameterAllometry.csv", ";", escape_double = FALSE, col_types = cols(alpha = col_number(), beta = col_number()), locale = locale(decimal_mark = ","), trim_ws = TRUE) %>% separate(ScientificName, sep = "_", into = c("Genus", "Species"), remove = F) %>% mutate(Family = ifelse(Taxo == "fam", Genus, NA)) %>% mutate(Genus = ifelse(is.na(Family), Genus, NA)) usethis::use_data(ParamCrownDiameterAllometry)
SpeciesCriteria <- read_csv("D:/VSC Kourou/DATA/SpeciesCriteria.csv", col_types = cols(Commercial = col_factor(levels = c("0", "1", "2")), MinFD = col_number(), UpMinFD = col_number(), MaxFD = col_number()), locale = locale(encoding = "Latin1")) %>% separate(ScientificName, sep = " ", into = c("Genus", "Species"), remove = T) %>% mutate(Aggregative = ifelse(VernName == "Angélique" | VernName == "Gonfolo rose", TRUE, FALSE)) # unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) SpeciesCriteria$VernName <- tolower(SpeciesCriteria$VernName) # usethis::use_data(SpeciesCriteria, overwrite = TRUE)
ScenariosTable <- read_csv("D:/VSC Kourou/DATA/ScenariosTable.csv", col_types = cols(Winching = col_factor(levels = c("0", "1", "2")), DirectionalFelling = col_factor(levels = c("0", "1", "2"))), trim_ws = TRUE) # usethis::use_data(ScenariosTable,overwrite = TRUE)
Data loading
data(Paracou6_2016) #data inventory in the global environment data(ForestZoneVolumeParametersTable) # The volume parameters data in the global environment # load(system.file("extdata", "BrokenParacou6_2016.rda", package = "Maria")) # the test inventory in the global environment data(ParamCrownDiameterAllometry) data(SpeciesCriteria)
lapply(Paracou6_2016, class) #give the class of all the data variables
inventorycheckformat
INPUT: inventory = Inventaires botaniques Guyafor/ONF (data frame) OUTPUT: inventory = Inventaires botaniques type “Guyafor”(data frame)
vérification du type d’inventaire : Test présence des variables. Test class des variables. Si les tests précédents ont tous reçu =TRUE, STOP de la fonction inventorycheckformat Sinon mettre en forme “inventaire Guyafor” (nom des var, leurs classes, unités (mettre des limites avec range),...) re-vérification du type d’inventaire : Test présence des variables Test class des variables Si les tests précédent ont tous reçu =TRUE, STOP de la fonction inventorycheckformat Sinon Message : “Votre format de données renseignés à l’argument “inventory”, n’est pas compatible. Voir la vignette du package pour y trouver le format requis.”
'pouet' %in% names(Paracou6_2016) if (!is.null(Paracou6_2016$pouet)) stop("pouet variable should be an integer.")
# Variables presence check if(!("Plot" %in% names(Paracou6_2016))) stop("Plot variable is not found") if(!("CensusYear" %in% names(Paracou6_2016))) stop("CensusYear variable is not found") if(!("idTree" %in% names(Paracou6_2016))) stop("idTree variable is not found") if(!("Family" %in% names(Paracou6_2016))) stop("Family variable is not found") if(!("Genus" %in% names(Paracou6_2016))) stop("Genus variable is not found") if(!("Species" %in% names(Paracou6_2016))) stop("Species variable is not found") if(!("CircCorr" %in% names(Paracou6_2016))) stop("CircCorr variable is not found") if(!("CodeAlive" %in% names(Paracou6_2016))) stop("CodeAlive variable is not found") if(!("CommercialSp" %in% names(Paracou6_2016))) stop("CommercialSp variable is not found") if(!("UTMZone" %in% names(Paracou6_2016))) stop("UTMZone variable is not found") if(!('Lat' %in% names(Paracou6_2016))) stop("Lat variable is not found") if(!('Lon' %in% names(Paracou6_2016))) stop("Lon variable is not found") if(!('VernName' %in% names(Paracou6_2016))) stop("VernName variable is not found") if(!('Xfield' %in% names(Paracou6_2016))) stop("Xfield variable is not found") if(!('Yfield' %in% names(Paracou6_2016))) stop("Yfield variable is not found") if(!('Xutm' %in% names(Paracou6_2016))) stop("Xutm variable is not found") if(!('Yutm' %in% names(Paracou6_2016))) stop("Yutm variable is not found")
# Variables class check if (!inherits(Paracou6_2016$idTree, "integer") && !is.null(idTree)) #majuscule ? stop("idTree variable should be an integer.") if (!inherits(Paracou6_2016$Plot, "character") && !is.null(Plot)) stop("Plot variable should be a character.") if (!inherits(Paracou6_2016$Xfield, "numeric") && !is.null(Xfield)) stop("Xfield variable should be numeric.") if (!inherits(Paracou6_2016$Yfield, "numeric") && !is.null(Yfield)) stop("Yfield variable should be numeric.") if (!inherits(Paracou6_2016$Xutm, "numeric") && !is.null(Xutm)) stop("Xutm variable should be numeric.") if (!inherits(Paracou6_2016$Yutm, "numeric") && !is.null(Yutm)) stop("Yutm variable should be numeric.") if (!inherits(Paracou6_2016$UTMZone, "integer") && !is.null(UTMZone)) stop("UTMZone variable should be an integer.") if (!inherits(Paracou6_2016$Lat, "numeric") && !is.null(Lat)) stop("Lat variable should be numeric.") if (!inherits(Paracou6_2016$Lon, "numeric") && !is.null(Lon)) stop("Lon variable should be numeric.") if (!inherits(Paracou6_2016$Family, "character") && !is.null(Family)) stop("Family variable should be a character.") if (!inherits(Paracou6_2016$Genus, "character") && !is.null(Genus)) stop("Genus variable should be a character.") if (!inherits(Paracou6_2016$Species, "character") && !is.null(Species)) stop("Species variable should be a character.") if (!inherits(Paracou6_2016$VernName, "character") && !is.null(VernName)) stop("VernName variable should be a character.") if (!inherits(Paracou6_2016$CommercialSp, "logical") && !is.null(CommercialSp)) stop("CommercialSp variable should be logical.") if (!inherits(Paracou6_2016$CensusYear, "integer") && !is.null(CensusYear)) stop("CensusYear variable should be an integer.") if (!inherits(Paracou6_2016$CodeAlive, "logical") && !is.null(CodeAlive)) stop("CodeAlive variable should be logical.") if (!inherits(Paracou6_2016$CircCorr, "numeric") && !is.null(CircCorr)) #DBH stop("CircCorr variable should be numeric.") # # if (!inherits(Paracou6_2016$Forest, "character")) # stop("Forest variable should be a character.") #pas besoin # if (!inherits(Paracou6_2016$idVern, "integer") && !is.null(idVern)) # stop("idVern variable should be an integer.") #pas besoin # # # # if (!inherits(Paracou6_2016$PlotArea, "numeric") && !is.null(PlotArea)) # stop("PlotArea variable should be numeric.") #pas besoin # # if (!inherits(Paracou6_2016$SubPlot, "integer") && !is.null(SubPlot)) # stop("SubPlot variable should be an integer.") #pas besoin # # if (!inherits(Paracou6_2016$TreeFieldNum, "numeric") && !is.null(TreeFieldNum)) # stop("TreeFieldNum variable should be numeric.") #pas besoin # # if (!inherits(Paracou6_2016$Project, "character") && !is.null(Project)) # stop("Project variable should be a character.") #pas besoin # # if (!inherits(Paracou6_2016$Protocole, "character") && !is.null(Protocole)) # stop("Protocole variable should be a character.") #pas besoin
if(!("Plot" %in% names(Paracou6_2016))) stop("Plot variable is not found") if(!("CensusYear" %in% names(Paracou6_2016))) stop("CensusYear variable is not found") if(!("idTree" %in% names(Paracou6_2016))) stop("idTree variable is not found") if(!("Family" %in% names(Paracou6_2016))) stop("Family variable is not found") if(!("Genus" %in% names(Paracou6_2016))) stop("Genus variable is not found") if(!("Species" %in% names(Paracou6_2016))) stop("Species variable is not found") if(!("CircCorr" %in% names(Paracou6_2016))) stop("CircCorr variable is not found") if(!("CodeAlive" %in% names(Paracou6_2016))) stop("CodeAlive variable is not found") if(!("CommercialSp" %in% names(Paracou6_2016))) stop("CommercialSp variable is not found") if(!("UTMZone" %in% names(Paracou6_2016))) stop("UTMZone variable is not found") if(!('Lat' %in% names(Paracou6_2016))) stop("Lat variable is not found") if(!('Lon' %in% names(Paracou6_2016))) stop("Lon variable is not found") if(!('VernName' %in% names(Paracou6_2016))) stop("VernName variable is not found") if(!('Xfield' %in% names(Paracou6_2016))) stop("Xfield variable is not found") if(!('Yfield' %in% names(Paracou6_2016))) stop("Yfield variable is not found") if(!('Xutm' %in% names(Paracou6_2016))) stop("Xutm variable is not found") if(!('Yutm' %in% names(Paracou6_2016))) stop("Yutm variable is not found") if (!inherits(Paracou6_2016$idTree, "integer.")) #majuscule ? stop("idTree variable should be an integer.") if (!inherits(Paracou6_2016$Plot, "character")) stop("Plot variable should be a character.") if (!inherits(Paracou6_2016$Xfield, "numeric")) stop("Xfield variable should be numeric.") if (!inherits(Paracou6_2016$Yfield, "numeric")) stop("Yfield variable should be numeric.") if (!inherits(Paracou6_2016$Xutm, "numeric")) stop("Xutm variable should be numeric.") if (!inherits(Paracou6_2016$Yutm, "numeric")) stop("Yutm variable should be numeric.") if (!inherits(Paracou6_2016$UTMZone, "integer")) stop("UTMZone variable should be an integer.") if (!inherits(Paracou6_2016$Lat, "numeric")) stop("Lat variable should be numeric.") if (!inherits(Paracou6_2016$Lon, "numeric")) stop("Lon variable should be numeric.") if (!inherits(Paracou6_2016$Family, "character")) stop("Family variable should be a character.") if (!inherits(Paracou6_2016$Genus, "character")) stop("Genus variable should be a character.") if (!inherits(Paracou6_2016$Species, "character")) stop("Species variable should be a character.") if (!inherits(Paracou6_2016$VernName, "character")) stop("VernName variable should be a character.") if (!inherits(Paracou6_2016$CommercialSp, "logical")) stop("CommercialSp variable should be logical.") if (!inherits(Paracou6_2016$CensusYear, "integer")) stop("CensusYear variable should be an integer.") if (!inherits(Paracou6_2016$CodeAlive, "logical")) stop("CodeAlive variable should be logical.") if (!inherits(Paracou6_2016$CircCorr, "numeric")) #DBH stop("CircCorr variable should be numeric.")
cleaninventory
INPUT: inventory = Inventaires botaniques type “Guyafor” (Paracou6_2016)(data frame) OUTPUT: inventory = Inventaires botaniques exploitable (Paracou6_2016Clean) (data frame)
Ne prendre que les arbres vivants : lignes ayant la valeur “TRUE” à la colonne “CodeAlive” de inventory. Check si les valeurs de ID sont uniques. Si FALSE message d’erreur “Les identifiants des arbres ne sont pas uniques” Enlever les arbres hors de la parcelle s’il y en a : coordonnées de l’arbre non comprises dans les coordonnées de la parcelle. Ne prendre que les lignes pour lesquelles “DBH”> 10.
duplicated(Paracou6_2016) unique(Paracou6_2016) #same rows number #duplicated() gives "FALSE" when the row is unique # any() : at least one of the values true? any(duplicated(Paracou6_2016)) # at least one row is duplicated? all(Paracou6_2016$Plot==Paracou6_2016$Plot[1]) #all the Plot values are equal?
if (!("DBH" %in% names(Paracou6_2016))) {add_column(Paracou6_2016, DBH = NA) #if DBH (cm) doesn't exist create it Paracou6_2016$DBH = Paracou6_2016$CircCorr/pi} # and compute it InventoryClean <- Paracou6_2016 %>% filter(CodeAlive == "TRUE") %>% #only alive trees filter(DBH >= 10) # DBH >= 10 if (any(duplicated(Paracou6_2016)==TRUE)) stop ("Tree identifiers (idTree) are not unique.") # stop function if the tree identifiers (idTree) are not unique if (all(Paracou6_2016$Plot==Paracou6_2016$Plot[1])==FALSE) #all the Plot values are equal? stop ("Input inventory concerns different plots (Plot). Our function simulates logging at the plot level.") # Remove out-plot trees (A faire)
addtreedim
Fonction d’ajout de variables pour l’INPUT “inventory” : INPUT: inventory = InventoryClean (data frame) OUTPUT: inventory = Inventory_addtreedim (data frame)
calculer pour tous les arbres : + TreeBiomass : biomasse de tous les arbres (BIOMASS package) + TreeHeight : hauteur de l’arbre (H) (BIOMASS package) + TreeHarvestableVolume : volume exploitable à partir du tarif de cubage (= a + bDBH2) que représente chaque arbre (a et b dépendent de la localisation) + TrunkHeight : hauteur de fût (TreeHarvestableVolume = π(DBH/2)² x TrunkHeight) + CrownDiameter : diamètre de couronne (CD) (ln(D) = 𝜶+ 𝜷 ln(HCD) + 𝜺 (allométries de Mélaine) + CrownHeight : TreeHeight - TrunkHeight
library(BIOMASS) library(knitr) TreeheightGuyafordata <- read_delim("D:/VSC Kourou/DATA/TreeheightGuyafordata.csv", ";", escape_double = FALSE, col_types = cols(Circ = col_number(), CircCorr = col_number(), Hauteur = col_number()), locale = locale(encoding = "Latin1"), trim_ws = TRUE) %>% mutate(HeightHow = recode(HeightHow, 'Visual estimate' = 'Visual Estimate')) %>% # filter(!(HeightHow == "Visual Estimate" | HeightHow == "Visual estimate")) %>% unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>% mutate(DBH = ifelse(is.na(CircCorr), Circ/pi, CircCorr/pi)) # add DBH column any(is.na(TreeheightGuyafordata$DBH)) # all the DBH have been computed? SineTelemeter <- TreeheightGuyafordata %>% #data to fit model filter(HeightHow == "Sine Method (Larjavaara & Muller-Landau 2013) with Trimble LaserAce 1000 Rangefinder"| HeightHow == "Telemeter") ModelsresultGuiana <-modelHD( D = SineTelemeter$DBH, H = SineTelemeter$Hauteur, useWeight = TRUE ) kable(ModelsresultGuiana) #log2 has the lowest RSE (4.699261 with all data, 4.316927 without "Visual E/estimate", 3.059860 with just "Sine Method", 3.809011 with sine and telemeter) HDmodelGuiana <- modelHD( D = SineTelemeter$DBH, H = SineTelemeter$Hauteur, method = "log2", # Compute the H-D model with the lowest RSE. useWeight = TRUE ) TreeHeightGuianaEstimation <- retrieveH( D = TreeheightGuyafordata$DBH, model = HDmodelGuiana ) #Put these tree heights estimations in the inventory TreeheightGuyafordataEstim <- TreeheightGuyafordata %>% mutate(HauteurEstimee = TreeHeightGuianaEstimation$H) plot(TreeheightGuyafordataEstim$DBH, TreeheightGuyafordataEstim$HauteurEstimee) # obs/estim regression: TreeheightGuyafordataEstim %>% ggplot() + aes(x = Hauteur, y = HauteurEstimee) + geom_point(shape = "circle", size = 0.5, colour = "#440154") + geom_smooth(span = 0.75) + theme_minimal() # We are well between 10 and 35m height
library(BIOMASS) library(knitr) # Just with Nouragues data (measurement: laser rangefinder): # Fit first different height-diameter models ModelsresultNouraguesdata <- modelHD( D = NouraguesHD$D, H = NouraguesHD$H, useWeight = TRUE # model weights will be (D^2)*H (weights are proportional to tree volume, so that larger trees have a stronger influence during the construction of the model) ) # choose the best model kable(ModelsresultNouraguesdata) #log2 has the lowest RSE (4.222718) (Root Squared Error) HDmodelNouraguesdata <- modelHD( D = NouraguesHD$D, H = NouraguesHD$H, method = "log2", # Compute the local H-D model with the lowest RSE. useWeight = TRUE ) # Applicate the model to retrieve tree heights ---- HERE TreeHeightEstimation <- retrieveH( D = Paracou6_2016$DBH, model = HDmodelNouraguesdata ) ## Put these tree heights estimations in the inventory Paracou6_2016TreeHeight <- Paracou6_2016 %>% mutate(TreeHeight = TreeHeightEstimation$H) plot(Paracou6_2016TreeHeight$DBH, Paracou6_2016TreeHeight$TreeHeight)#amazing !! It's a log2 form # Different guyanese data: TreeheightGuyafordata <- read_delim("D:/VSC Kourou/DATA/TreeheightGuyafordata.csv", ";", escape_double = FALSE, col_types = cols(Circ = col_number(), CircCorr = col_number(), Hauteur = col_number()), locale = locale(encoding = "Latin1"), trim_ws = TRUE) %>% # filter(!(HeightHow == "Visual Estimate" | HeightHow == "Visual estimate")) filter(HeightHow == "Sine Method (Larjavaara & Muller-Landau 2013) with Trimble LaserAce 1000 Rangefinder") if (!("DBH" %in% names(TreeheightGuyafordata))) {add_column(TreeheightGuyafordata, DBH = NA) #if DBH (cm) doesn't exist create it TreeheightGuyafordata$DBH = TreeheightGuyafordata$CircCorr/pi} # and compute it unique(TreeheightGuyafordata$HeightHow) ModelsresultGuiana <- modelHD( D = TreeheightGuyafordata$DBH, H = TreeheightGuyafordata$Hauteur, useWeight = TRUE ) kable(ModelsresultGuiana) #log2 has the lowest RSE (4.749107 with all data, 4.316927 without "Visual E/estimate", 3.059860 with just "Sine Method") #Don't forget, you have the model form, and its parameters, and there are different according to the data. HDmodelGuiana <- modelHD( D = TreeheightGuyafordata$DBH, H = TreeheightGuyafordata$Hauteur, method = "log2", # Compute the local H-D model with the lowest RSE. useWeight = TRUE ) TreeHeightGuianaEstimation <- retrieveH( D = Paracou6_2016$DBH, model = HDmodelGuiana ) ## Put these tree heights estimations in the inventory Paracou6_2016TreeHeight <- Paracou6_2016 %>% mutate(TreeHeight = TreeHeightGuianaEstimation$H) plot(Paracou6_2016TreeHeight$DBH, Paracou6_2016TreeHeight$TreeHeight)#amazing !! It's a log2 form
# Question : ya t'il des arbres mesurés par plsrs méthodes la même année ? TreeheightMethods <- TreeheightGuyafordata %>% select(Forest, Plot, SubPlot, TreeFieldNum, Genus, Species, CensusYear) %>% unique()
# TreeHarvestableVolume = volume exploitable à partir du tarif de cubage (= a + b*DBH2) que représente chaque arbre (a et b dépendent de la localisation) ForestZoneVolumeParametersTable <- data.frame (Forest = c("Acarouany","BAFOG","Kaw","Laussat","Montagne Plomb","Montagne Tortue","Nouragues","Organabo","Paracou","Régina St Georges","Risquetout","Tibourou"), # Guyafor forests Zone = c("East", "West", "Central", "FrenchGuiana")) %>% # Volume formules zones mutate(Zone = ifelse(Forest == "Acarouany"| Forest =="BAFOG" |Forest =="Laussat", "West", NA)) %>% # West zone mutate(Zone = ifelse(Forest == "Montagne Plomb"| Forest =="Organabo"|Forest =="Paracou"|Forest =="Risquetout", "Central", Zone)) %>% # Central zone mutate(Zone = ifelse(Forest == "Kaw"| Forest =="Montagne Tortue"|Forest =="Nouragues"|Forest =="Régina St Georges"|Forest =="Tibourou", "East", Zone)) %>% # East zone left_join(VolumeParametersTable) Paracou6_2016_volume <- Paracou6_2016 %>% left_join(ForestZoneVolumeParametersTable) %>% mutate(TreeHarvestableVolume = aCoef + bCoef * CircCorr^2) # compute the tree harvestable volume # West = c("Acarouany", "BAFOG", "Laussat") # Central = c("Montagne Plomb", "Organabo","Paracou", "Risquetout") # East = c("Kaw", "Montagne Tortue","Nouragues", "Régina St Georges", "Tibourou") # Guyane = le reste par défaut # et possibilité à l'utilisateur de mettre ses propres tarifs
mutate(TrunkHeight = TreeHarvestableVolume/(π(DBH/2)^2)) # compute the trunk height (cylinderVolume= π(DBH/2)² x H) mutate(CrownHeight = TreeHeight - TrunkHeight) # compute the crown height
spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale filter(Taxo == "sp") %>% select(-Family) genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale filter(Taxo == "gen") %>% select(-Family, -ScientificName, -Species) famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale filter(Taxo == "fam") %>% select(-Genus, -ScientificName, -Species) Paracou6_2016 %>% sample_n(100) %>% # short data select(idTree, Family, Genus, Species) %>% #only useful columns left_join(select(spParamCrownDiameter, -Taxo, -ScientificName, -beta), by = c("Genus","Species")) %>% left_join(select(genParamCrownDiameter, -Taxo, -beta), by = c("Genus"), suffix = c(".species", ".genus")) %>% left_join(select(famParamCrownDiameter, -Taxo, -beta) %>% rename(alpha.family = alpha), by = c("Family")) %>% mutate(alpha = alpha.species) %>% #mutate(alpha = alpha.species,beta = beta.species) mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>% mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>% mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>% select(-alpha.species, -alpha.genus, -alpha.family)
spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale filter(Taxo == "sp") %>% select(-Family) genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale filter(Taxo == "gen") %>% select(-Family, -ScientificName, -Species) famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale filter(Taxo == "fam") %>% select(-Genus, -ScientificName, -Species) error <- rnorm(2000, mean = 0, sd = 0.0295966977) #Mélaine's model error distribution (hist(error)) Paracou6_2016_crown <- Paracou6_2016 %>% unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>% left_join(spParamCrownDiameter, by = c("ScientificName","Genus","Species")) %>% #Add parameters at species scale left_join(genParamCrownDiameter, by = "Genus", suffix = c(".species", ".genus")) %>% #Add parameters at genus scale left_join(famParamCrownDiameter, by = "Family") %>% #Add parameters at family scale dplyr::rename(Taxo.family = Taxo, alpha.family = alpha, beta.family = beta) %>% mutate(alpha = alpha.species, beta = beta.species, Taxo = Taxo.species) %>% # create the future & unique alpha,beta,taxo columns mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>% #if species parameters are absent, take genus parameters mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>% #if species&genus parameters are absent, take family parameters mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>% #if species,genus&family parameters are absent, take parameters mean select(-alpha.species, -alpha.genus, -alpha.family) %>% #remove obsolete columns mutate(beta = ifelse(is.na(beta), beta.genus, beta)) %>% #if species parameters are absent, take genus parameters mutate(beta = ifelse(is.na(beta), beta.family, beta)) %>% #if species&genus parameters are absent, take family parameters mutate(beta = ifelse(is.na(beta), mean(beta, na.rm =T), beta)) %>% #if species,genus&family parameters are absent, take parameters mean select(-beta.species, -beta.genus, -beta.family) %>% #remove obsolete columns mutate(Taxo = ifelse(is.na(Taxo), Taxo.genus, Taxo)) %>% #if species parameters are absent, take genus parameters mutate(Taxo = ifelse(is.na(Taxo), Taxo.family, Taxo)) %>% #if species&genus parameters are absent, take family parameters mutate(Taxo = ifelse(is.na(Taxo), "mean", Taxo)) %>% #if species,genus&family parameters are absent, take parameters mean select(-Taxo.species, -Taxo.genus, -Taxo.family) %>% #remove obsolete columns mutate(CrownDiameter = exp(( (log(DBH)-alpha-error)/beta) )/TreeHeight) # compute the crown diameter (CD) (ln(D) = 𝜶+ 𝜷 ln(H*CD) + 𝜺, with 𝜺~N(0,σ^2) and meanσ^2 = 0.0295966977. (Mélaine's allometries)) # Informations: # How many individuals have find their Crown diameter parameters at species, genus, family or no scale: sum(Paracou6_2016_crown$Taxo == "fam") #442 individuals with family scale parameters sum(Paracou6_2016_crown$Taxo == "gen") #742 individuals with genus scale parameters sum(Paracou6_2016_crown$Taxo == "sp") #2364 individuals with family scale parameters sum(Paracou6_2016_crown$Taxo == "mean") #72 individuals no recognize at any scale, so with plot mean parameters unique(Paracou6_2016_crown[(Paracou6_2016_crown$Taxo == "mean"),]$ScientificName) #which species are no recognized length(unique(Paracou6_2016_crown$ScientificName)) #How many species are in the plot
# allometry parameters data preparation: spParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at species scale filter(Taxo == "sp") %>% select(-Family) genParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at genus scale filter(Taxo == "gen") %>% select(-Family, -ScientificName, -Species) famParamCrownDiameter <- ParamCrownDiameterAllometry %>% #parameters at family scale filter(Taxo == "fam") %>% select(-Genus, -ScientificName, -Species) error <- rnorm(2000, mean = 0, sd = 0.0295966977) #Mélaine's model error distribution (hist(error)) # Variables computation: InventoryAddtreedim <- InventoryClean %>% # TreeHarvestableVolume (m^3) left_join(ForestZoneVolumeParametersTable) %>% mutate(TreeHarvestableVolume = aCoef + bCoef * (DBH/100)^2) %>% # compute the tree harvestable volume # TrunkHeight & CrownHeight (m) mutate(TrunkHeight = TreeHarvestableVolume/(pi*(((DBH/100)/2)^2))) %>% # compute the trunk height (cylinderVolume= π(DBH/2)² x H) DBH in cm, in m in the formule. # mutate(CrownHeight = TreeHeight - TrunkHeight) %>% # compute the crown height TREEHEIGHT EN ATTENTE # CrownDiameter (m) # Add crown diameter parameters at the inventory: unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>% left_join(spParamCrownDiameter, by = c("ScientificName","Genus","Species")) %>% #Add parameters at species scale left_join(genParamCrownDiameter, by = "Genus", suffix = c(".species", ".genus")) %>% #Add parameters at genus scale left_join(famParamCrownDiameter, by = "Family") %>% #Add parameters at family scale rename(Taxo.family = Taxo, alpha.family = alpha, beta.family = beta) %>% mutate(alpha = alpha.species, beta = beta.species, Taxo = Taxo.species) %>% # create the future & unique alpha,beta,taxo columns mutate(alpha = ifelse(is.na(alpha), alpha.genus, alpha)) %>% #if species parameters are absent, take genus parameters mutate(alpha = ifelse(is.na(alpha), alpha.family, alpha)) %>% #if species&genus parameters are absent, take family parameters mutate(alpha = ifelse(is.na(alpha), mean(alpha, na.rm =T), alpha)) %>% #if species,genus&family parameters are absent, take parameters mean select(-alpha.species, -alpha.genus, -alpha.family) %>% #remove obsolete columns mutate(beta = ifelse(is.na(beta), beta.genus, beta)) %>% #if species parameters are absent, take genus parameters mutate(beta = ifelse(is.na(beta), beta.family, beta)) %>% #if species&genus parameters are absent, take family parameters mutate(beta = ifelse(is.na(beta), mean(beta, na.rm =T), beta)) %>% #if species,genus&family parameters are absent, take parameters mean select(-beta.species, -beta.genus, -beta.family) %>% #remove obsolete columns mutate(Taxo = ifelse(is.na(Taxo), Taxo.genus, Taxo)) %>% #if species parameters are absent, take genus parameters mutate(Taxo = ifelse(is.na(Taxo), Taxo.family, Taxo)) %>% #if species&genus parameters are absent, take family parameters mutate(Taxo = ifelse(is.na(Taxo), "mean", Taxo)) %>% #if species,genus&family parameters are absent, take parameters mean select(-Taxo.species, -Taxo.genus, -Taxo.family) #remove obsolete columns # mutate(CrownDiameter = exp(( # (log(DBH)-alpha-error)/beta) # )/TreeHeight) # compute the crown diameter (CD) (ln(D) = 𝜶+ 𝜷 ln(H*CD) + 𝜺, with 𝜺~N(0,σ^2) and meanσ^2 = 0.0295966977. (Mélaine's allometries)) # Mettre finalement l'objet créé dans le zzz.R # !!!!!!!TREEHEIGHT EN ATTENTE!!!!!!
scenariosparameters
pas de mode "manual", si argument = NULL, prendre la valeur du scénario correspondante
if (is.null(winching)){ if (Type == "RIL1") winching = "0" if (Type == "RIL2broken") winching = "0" if (Type == "RIL2") winching = "1" if (Type == "RIL3") winching = "2" if (Type == "RIL3fuel") winching = "2" if (Type == "RIL3fuelhollow") winching = "2" } if (is.null(directionalfelling)){ if (Type == "RIL1") directionalfelling = "0" if (Type == "RIL2broken") directionalfelling = "1" if (Type == "RIL2") directionalfelling = "1" if (Type == "RIL3") directionalfelling = "2" if (Type == "RIL3fuel") directionalfelling = "2" if (Type == "RIL3fuelhollow") directionalfelling = "2" } if (is.null(objective)){ if (Type == "RIL1") objective = 25 if (Type == "RIL2broken") objective = 25 if (Type == "RIL2") objective = 25 if (Type == "RIL3") objective = 30 if (Type == "RIL3fuel") objective = 30 if (Type == "RIL3fuelhollow") objective = 30 } if (is.null(diversification)){ if (Type == "RIL1") diversification = FALSE if (Type == "RIL2broken") diversification = FALSE if (Type == "RIL2") diversification = FALSE if (Type == "RIL3") diversification = TRUE if (Type == "RIL3fuel") diversification = TRUE if (Type == "RIL3fuelhollow") diversification = TRUE } # Type SpatialDataType Winching DirectionalFelling Objective Diversification # RIL1 SRTM 0 0 20-25 FALSE # RIL2broken LIDAR 0 1 20-25 FALSE # RIL2 LIDAR 1 1 20-25 FALSE # RIL3 LIDAR 2 2 25-30 TRUE # RIL3fuel LIDAR 2 2 25-30 TRUE # RIL3fuelhollow LIDAR 2 2 25-30 TRUE
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.