Nicola Puletti 17 dicembre 2018 @ Arezzo

Studio su estrazioni TLS

# General scope Assessment of poplar profiles sensitivity to competition anisotropy and to plantation density. Three plots have been established in a poplar plantation of 10 years in Viadana (Mantova). Poplar rows are alternated with rows of other valuable (but less fast growing) species, one every two poplar rows. Distances from stem to stem within the plantation rows vary between the plots: 4 m, 4.5 m and 5 m. Distance between the the plantation rows is 9 m. Each plot includes 12 stems from two facing rows. Based on TLS-data, stem profiles from two orthogonal vertical crosssections have been determined: one along the rows direction and the other on the between-rows direction. Profiles have been sampled systematically adopting different segment lengths (Segm_lenght) estimating stem crosssection width (diam_btw_rows and diam_wti_rows) at different heigths along the stem (Sect_height)

[MEMO: credo sia opportuno ricordarsi di guardare Rubio-Cuadrado Á, Bravo-Oviedo A, Mutke S, Del Río M (2018). Climate effects on growth differ according to height and diameter along the stem in Pinus pinaster Ait. iForest 11: 237-242. – doi: 10.3832/ifor2318-011 [online 2018-03-12]]

Setup and fetch the data

library(rgl); library(hexbin); library(magrittr); library(reshape2); library(tibble)
library(geometry); library(rgeos); library(sp); library(ggplot2); library(readxl)

# (subFolds <- list.files('H:/TLS Plantations/AAA_StemAnalysis/stem_ruotati'))
SCDfolder <- "StemCloudData"
(subFolds <- list.files(SCDfolder))

## [1] "V400" "V450" "V500"

Hvox <- .10

started.at=proc.time()
out <- NULL
# for(k in 5:7) {
#  sF <- subFolds[k]
for(sF in subFolds) {
  myf <- paste(paste(SCDfolder, sF, sep='/'), list.files(paste(SCDfolder, sF, sep='/')), sep='/')
  
  for( i in 1:length(myf) ) {
    
    inData <- read.table(myf[i], sep = '\t')
    
    colnames(inData) <- c('X','Y','Z')
    
    inData$Xn <- inData$X-min(inData$X)
    inData$Yn <- inData$Y-min(inData$Y)
    inData$Zn <- inData$Z-min(inData$Z)
    
    Htop <- max(inData$Zn)
    
    bT <- seq(0,Htop,by=Hvox); aT <- (bT+Hvox)
    
    dx <- dy <- NULL
    for(j in 1:length(bT)) {
      
      sTj <- inData[inData$Zn>=bT[j] & inData$Zn<aT[j], c('Xn', 'Yn', 'Zn')]
    
      xytab <- data.frame(sTj$Xn, sTj$Yn)
      POL.j <- readWKT( sprintf("MULTIPOINT (%s )", paste(sprintf("(%f %f)", xytab[,1], xytab[,2]), collapse = ", ")) )
      convGeo <- gConvexHull(POL.j)
      mx <- my <- max(c(max(inData$Xn),max(inData$Yn)))
      
      dx <- c(dx, convGeo@bbox[1,2]-convGeo@bbox[1,1])
      dy <- c(dy, convGeo@bbox[2,2]-convGeo@bbox[2,1])
    }
    #tesi = substr(out$stem_id,1,4), stem = substr(out$stem_id,6,8
    out <- rbind(out, data.frame(tesi = sF, treid = substr(myf[i],21,22), length_toppo=Hvox, slice = bT+Hvox/2, diam_btw_rows=dx, diam_wti_rows=dy))
  }
}
cat("Finished in", data.table::timetaken(started.at),"\n")

## Finished in 24.9sec

Diameters <- out %>% 
  melt(id.vars=c('tesi', 'treid', 'slice', 'length_toppo')) %>%
  rename(Treat = 1, TreeId = 2, direction = 5, diam = 6)
#  colnames()[c(1,2,5,6)] <- c('Treat','TreeId','direction', 'diam')

Diameters %>%
  ggplot(aes(slice, diam)) + 
  ggtitle(paste('slice height:', Hvox*100, 'cm')) +
  # scale_y_continuous(limits = c(0.25, 0.5)) +
  scale_x_continuous(limits = c(0, 4.5)) +
  xlab('cross-section height (m)') +
  geom_line(aes(color=direction)) +
  facet_grid(TreeId ~ Treat, scales = 'free')

## Warning: Removed 22 rows containing missing values (geom_path).

right = function(text, num_char) {
  substr(text, nchar(text) - (num_char-1), nchar(text))
}

#stems_basic_measurements <-   readxl::read_xlsx("C:/Users/Workstation01/OneDrive - CREA/AIT2018 short comm
stems_basic_measurements <- 
  readxl::read_xlsx("Pioppi_Viadana_CaratteriDendrometriciBase.xlsx") 
# %>%
  # dplyr::mutate(TreeId = paste0("s", right(as.character(1000+TreeId),2))) %>%
  # dplyr::mutate_at(c("Tesi", "TreeId"), as.factor)

# Diameters0 <- data.frame(tesi = substr(out$stem_id,1,4), stem = substr(out$stem_id,6,8), out[,2:5])
# Diameters0 <- Diameters0 %>% dplyr::rename(treid=stem)
# expData3 <- Diameters0[,c("tesi","length_toppo","treid","slice","diam_btw_rows","diam_wti_rows")]
# SOSTANZIALMENTE  out -> Diameters0 -> expData3 -> tff, non mi sembra che i passaggi intermedi servano!

tff <- out %>%
  # riodina e rinomina colonne
  select(
         Treat = tesi,               # distance between rows: Treatment
         Segm_length = length_toppo, # length of the segments
         TreeId = treid,
         Sect_height = slice,        # height along the stem of diameter measurements (cross Sections)
         diam_btw_rows, 
         diam_wti_rows 
         ) %>%    
# ATTENZIONE! Non ho capito questo 'dbh'!!
# expData3$dbh <- apply(expData3[, c('diam_btw_rows','diam_wti_rows')],1,mean)
# a me sembra un 'avg_diam' a quota 'slice', che non è necessariamente 1.3!!
  mutate(avg_diam = (diam_btw_rows + diam_wti_rows)/2,
         delta_d_cm = 100 * (diam_btw_rows - diam_wti_rows)) %>% 
  # Trees in different treatments are different subjects, better using distinct Ids
  # mutate(TreeId = paste0(case_when(tesi == "V400" ~ "A",
  #                                  tesi == "V450" ~ "B",
  #                                  tesi == "V500" ~ "C"),
  #                        right(as.character(treid),2))) %>%
#  data.frame() %>%
  filter(Sect_height <= 5.5) %>% 
  inner_join(mutate(stems_basic_measurements,
                              Treat = as.factor(Tesi),
                              TreeId = as.factor(str_sub(stems_basic_measurements$TreeId +100, 2,3)),
                              dbh = DBH))

## Joining, by = c("Treat", "TreeId")

             # ,
             # by = c( "Treat" = "Treat", "treid" = "treid"))
head(tff)

library(tidyverse)
library(magrittr)
library(ggplot2)
library(knitr)
library(broom)

set.seed(432); tff %>%
  select(Treat, TreeId, dbh) %>%
  unique() %>%
  ggplot(aes(y = dbh, x = Treat)) + geom_boxplot(colour = "darkgrey") +
  geom_jitter(width = .05, height = 0) + 
  ggtitle("Effect of plantation density on brest height diametr (dbh)") +
  ylab("dbh [cm]")

tff %>%
  select(Treat, TreeId, dbh = DBH) %>%
  unique() %>%
  lm(dbh ~ Treat, data = .) %>%
  summary() %>% tidy() %>%
  kable(digits = 3, format.args = list(zero.print = F),
        caption = "Plantation density effect on 'dbh'")

TLSdbh_df <- tff %>% 
  dplyr::filter(abs(1.3 - Sect_height)<.05) %>% 
  select(Treat, TreeId, dbh, starts_with('diam')) %>%
  gather(direction, TLSdbh, starts_with('diam')) %>%
  mutate(TLSdbh = 100 * TLSdbh)

TLSdbh_df %>%
  lm(TLSdbh ~ Treat + direction, data = .) %>%
  summary() %>% tidy() %>%
  kable(digits =c(0, 2, 3, 3, 5), 
        caption = "Plantation density and measurement directions effects on TLS estimated dbh")

To be careful about!

TLSdbh_df %>% 
  spread(direction, TLSdbh) %>%
  mutate( dbh = jitter(dbh, factor = 1)) %>%
  ggplot(aes(x = dbh, color = Treat)) +
  geom_abline(slope=1, color = "grey") +
  geom_errorbar(aes(ymin=diam_wti_rows, ymax=diam_btw_rows), 
                width = .2, linetype =2, size = .5) +
  geom_line(aes(y = I((diam_wti_rows + diam_btw_rows)/2))) +
  geom_point(aes(y = I((diam_wti_rows + diam_btw_rows)/2)), shape = 1) +
  scale_x_continuous(breaks = seq(0, 100, 1), minor_breaks = NULL) +
  scale_y_continuous(breaks = seq(0, 100, 1), minor_breaks = NULL) +
  ggtitle("Field measured dbh and corresponding TLS estimates",
          subtitle = "between- and within-rows estiamates by plantation density")

Profile’s difference analysed by Hierarchical Linear Modeling (HLM)

Tree level modelling has been developed based on crossesctions estimated dividing the stem in segments of 0.1 meters.

oll <- 0.1 # optimal Log Length for profiles discrimination
# Segm_length ottimale: compromesso tra min lag correlation e conservazione del dettaglio della forma
stff <- tff  %>%
  filter(Segm_length == oll)

library(lme4)
library(HLMdiag)

mc <- function(frml, main = paste0(as.expression(frml)), ...) {
# evaluate the model (specified as formula) and produce the residulas control plot
  fm <- stff %>%  lmer(data = ., frml, ... )
  p <- qplot(x = stff$Sect_height, y = residuals(fm), 
        geom = c("point", "smooth")) + 
    xlab("cross-section height [m]") +
    ylab("level-1 residuals") +
    ggtitle(main)
  print(p)
return(fm)
}

# prova CON Sect_height, al cubo (con RML=T, non converge!)
fm13 <- mc(delta_d_cm ~ Treat + (1 | TreeId) 
           + Sect_height *(1 | TreeId) 
           + I(Sect_height^2)* (1 | TreeId)
           + I(Sect_height^3)* (1 | TreeId)
           , REML = F
           , main = "delta_d ~ Treat + (poly(Sect_height, 3)| TreeId)")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

stff %>% ggplot() +
  geom_point(aes(x = Sect_height, y = delta_d_cm)) +
  geom_smooth(aes(x = Sect_height, y = delta_d_cm)) +
  geom_line(aes(x = Sect_height, y = fitted(fm13)), color = "yellow", size = 1) + 
  facet_grid(rows = vars(TreeId), cols = vars(Treat), scales = "free") +
  ggtitle("Nested models fitting profile differences by tree")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

resid2_fm13 <- HLMresid(object = fm13, level = "TreeId")
names(resid2_fm13) <- paste0("pw",0:3,"coeff")
# head(resid2_fm13)

stff %>% 
  select(Treat, TreeId) %>% 
  unique() %>% 
  inner_join(mutate(resid2_fm13, TreeId = row.names(resid2_fm13))) %>%
  gather(param, value,  3:6) %>%
  ggplot(aes(x = Treat, y = value)) +
  geom_boxplot(color = "darkgrey") +
  geom_jitter(width = .15, height = 0) + 
  geom_hline(yintercept = 0) +
  facet_wrap(~param, scales = "free") +
  ggtitle("Level 2 (Treat = plantation distance) residuals")

## Joining, by = "TreeId"

## Warning: Column `TreeId` joining factor and character vector, coercing into
## character vector

Treat_anova <- stff %>%  
  lmer(data = ., delta_d_cm ~ 
             1 * (1 | TreeId) 
           + Sect_height *(1 | TreeId) 
           + I(Sect_height^2)* (1 | TreeId)
           + I(Sect_height^3)* (1 | TreeId)
           , REML = F) %>%
  anova(fm13)

Treat_anova %>%
  mutate_all(funs(ifelse(is.na(.), 0, .))) %>%
  kable(digits = c(0, 2, 2, 2, 2, 4, 0, 4), format.args = list(zero.print = F), 
    caption = cat(paste("Contribution of plantation distance (Treat) to profiles differences variance explaination", paste(attributes(Treat_anova)$heading[-1], collapse = "\n"), collapse = "\n")))

## Contribution of plantation distance (Treat) to profiles differences variance explaination Models:
## .: delta_d_cm ~ 1 * (1 | TreeId) + Sect_height * (1 | TreeId) + 
## .:     I(Sect_height^2) * (1 | TreeId) + I(Sect_height^3) * (1 | 
## .:     TreeId)
## fm13: delta_d_cm ~ Treat + (1 | TreeId) + Sect_height * (1 | TreeId) + 
## fm13:     I(Sect_height^2) * (1 | TreeId) + I(Sect_height^3) * (1 | 
## fm13:     TreeId)

fm13 %>%
  augment() %>% 
  ggplot(aes(x = Treat, y = .fitted, colour = Treat)) + 
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", shape = 13, size = 4) +
  ylab("BetweenRowsDiam - WithinRowsDiameter [cm] Model estimates") +
  ggtitle("In the direction with less competition stem cross-sections are 1.6-1.8 cm wider but,\n even accounting for tree-level effects, the variability of the difference between orthogonal profiles\n still completely hides the effect of plantation density (means = crossed circles)")