In thomasgaquiere/Maria: Simulate Different Forms of Logging on a Forest

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) %>% 
  unite(Genus, Species, col = "ScientificName", sep = "_", remove = F) %>% 
  mutate(DBH = ifelse(is.na(CircCorr), Circ/pi, CircCorr/pi)) %>% 
  mutate(HeightHow = recode(HeightHow, 'Visual estimate' = 'Visual Estimate'))

any(is.na(TreeheightGuyafordata$DBH)) # all the DBH have been computed?

# DBH range per method
TreeheightGuyafordata %>% 
  group_by(HeightHow) %>% #pour réaliser les opérations suivantes par groupes de modalités
  summarise(MinDBH = min(DBH)) %>% 
  arrange(MinDBH)

TreeheightGuyafordata %>% 
  group_by(HeightHow) %>% #pour réaliser les opérations suivantes par groupes de modalités
  summarise(MaxDBH = max(DBH)) %>% 
  arrange(desc(MaxDBH))
# Visual Estimate: 0.416986 - 256.23946
# Telemeter: 9.867606 - 159.15494
# Optical Telemeter Birnbaum: 7.591691 -131.46198
# Vertex Laser: 10.185916 - 95.73170
# Sine Method: 6.073353 - 88.80846


# measurements number per method
TreeheightGuyafordata %>% 
  count(HeightHow) %>% 
  arrange(desc(n))
# Visual Estimate   : 27986         
# Optical Telemeter Birnbaum: 4902  
# Telemeter: 1454           
# Sine Method:  998 
# Vertex Laser: 441

TreeheightGuyafordata %>% 
  group_by(HeightHow) %>% 
  n_distinct(ScientificName) %>% 
  count(ScientificName) %>% 
  arrange(desc(n))

TreeheightGuyafordata %>% 
  group_by(HeightHow) %>%
  n_distinct(TreeheightGuyafordata$ScientificName) %>% 
  count(ScientificName)

#Methods facet

# png(file = "D:/VSC Kourou/TreeHeights/MethodsFacet.png", width=1500, height=612)
TreeheightGuyafordata %>%
 ggplot() +
  aes(x = DBH, y = Hauteur) +
  geom_point(shape = "circle", size = 0.5, colour = "#440154") +
  geom_smooth(span = 0.75) + # using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
  theme_minimal() +
  facet_wrap(vars(HeightHow))

# dev.off()

# Colored multilignes
# png(file = "D:/VSC Kourou/TreeHeights/Colored multilignes.png", width=1500, height=612)

TreeheightGuyafordata %>%
  filter(HeightHow != "Visual Estimate") %>% #je l'enlève pcq ça cache tout
 ggplot() +
  aes(x = DBH, y = Hauteur, colour = HeightHow) +
  geom_line(size = 0.2) +
  scale_color_hue(direction = 1) +
  theme_minimal() +
  theme(legend.position = "bottom")

# dev.off()

# Methods Box-Plot
# png(file = "D:/VSC Kourou/TreeHeights/Methods Box-Plot.png", width=1500, height=612)

TreeheightGuyafordata %>%
ggplot(aes(x=DBH, y=Hauteur, fill=HeightHow)) +
  theme(legend.position = "bottom") +
  geom_boxplot()

# dev.off()

# DBH distribution per method
# png(file = "D:/VSC Kourou/TreeHeights/DBH distribution per method.png", width=1500, height=612)

TreeheightGuyafordata %>%
 ggplot() +
  aes(x = DBH) +
  geom_histogram(bins = 30L, fill = "#440154") +
  theme_minimal() +
  facet_wrap(vars(HeightHow), scales = "free") #or free_y, "fixed", or "free"

# dev.off()

# Ind/sp distribution per method
# png(file = "D:/VSC Kourou/TreeHeights/Ind-sp distribution per method.png", width=1500, height=612)

ggplot(TreeheightGuyafordata) +
  aes(x = ScientificName) +
  geom_bar(fill = "#440154") +
  theme_minimal() +
  facet_wrap(vars(HeightHow), scales = "free") #or "fixed"

# dev.off()

library(BIOMASS)
library(knitr)

ModelsresultGuiana <-modelHD(
  D = TreeheightGuyafordata$DBH,
  H = TreeheightGuyafordata$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")
#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 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)#amazing !! It's a log2 form

#facet obs/estim :
# png(file = "D:/VSC Kourou/TreeHeights/facet obs-estim.png", width=1500, height=612)

TreeheightGuyafordataEstim %>%
 ggplot() +
  aes(x = Hauteur, y = HauteurEstimee) +
  geom_point(shape = "circle", size = 0.5, colour = "#440154") +
  geom_smooth(span = 0.75) + # try with method = lm
  theme_minimal() +
  facet_wrap(vars(HeightHow))

# dev.off()

# Fit Height-Diameter model
library(BIOMASS)

SineTelemeter <- 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 == "Sine Method (Larjavaara & Muller-Landau 2013) with Trimble LaserAce 1000 Rangefinder"| HeightHow == "Telemeter") %>%
  mutate(DBH = ifelse(is.na(CircCorr), Circ/pi, CircCorr/pi)) # add DBH column


HDmodel <- modelHD(
  D = SineTelemeter$DBH,
  H = SineTelemeter$Hauteur,
  method = "log2", # Compute the H-D model with the lowest RSE.
  useWeight = TRUE
)

Dmodel[["coefficients"]]

# log 2: log(H) = a + b*log(D) + c*log(D)^2

# log(H) = 0.07359191 + 1.34241216*log(D) + -0.12282344*log(D)^2
# H = exp(0.07359191 + 1.34241216*log(D) + -0.12282344*log(D)^2)


# (Intercept)  0.07359191
# I(log(D)^1)  1.34241216
# I(log(D)^2) -0.12282344