R/PrepareTreeMap.R

In EmmanuelDuchateau/TLS: Provides Tools For Terrestrial LiDAR Data Analysis

#' Extract a slice of point cloud around DBH height for TreeMap
#'
#' Extract a slice centered on the DBH height (1.3m) with a thickness of 2* ThickDHP
#'
#' @param Mydata this dataset should contain at least 3 columns with X, Y, Z values.
#' The .las (.laz) format is better for faster processing but a data.frame (or data.table) is acceptable.
#' @param ThickDHP indicate the thickness of the slice around the DBH (1.3m height) (10cm by default, so the thick of the slice will be 20cm)
#' @param Hsoil Needed for the function Zsoil. Indicates the height of the points held above ground to create the DTM.
#' If the field is generally flat, the 4 m by default are sufficient, however for a fairly steep terrain,
#' it may be necessary to increase the value of Zsoil.
#'
#' @return Returns a data.table containing the points X, Y, Z, and Zr (created with Zsoil) located between (1.3 - ThickDHP) & (1.3 + ThickDHP)
#' @export
#'
DataInit<-function(Mydata, ThickDHP = 0.1, Hsoil=4){
  #Calculate real Z values#####
  tot2<-Zsoil(Mydata)
  tot2<-setDT(tot2@data)
  #####
  dat4<-tot2[Zr>=(1.3 - ThickDHP) & Zr<=(1.3 + ThickDHP)]
  #
  dat4<-dat4 %>%
    mutate(xa = round(X,2),
           ya = round(Y,2),
           za = round(Zr,2))
  dat4<-dat4 %>%
    group_by(xa,ya)%>%
    arrange(Z)%>%
    mutate(selec=seq(1,n(),1),
           alph=n())

  dat4<-dat4 %>%
    mutate(alph1=round(alph/2)+1,
           selec1=selec-alph1)

  dat4<-setDT(dat4)
  message("Calculate the real Z values and extract a slice of the points located around the DBH height (Zr = 1.3m) ... OK")

  return(dat4)
}




#' Find the best threshold (Alpha)
#'
#' Generated graphs for different threshold values to select the best value for cleaning the point cloud and make it usable for the TreeMap.
#' The selected final threshold value should be used to remove single points and branch points. Only isolated trunks should remain visible
#'
#' @param data must be the data generated by the function DataInit to work
#' @param standName Optional. If the tested plot has a name, it will be added to the name of the graphics (name = TestStd by default).
#' @param Threshold Indicates the threshold value to be tested (can be a vector of several values). A look at the graphs generated will allow to adjust the desired values and retest them.
#' @param saveGraphPath indicates the folder where the graphs will be generated.
#'
#' @return Returns graphs stored in the chosen folder
#' @export
#'
testAlphaThreshold<-function (data, standName=NULL, Threshold = c(3,5,8,10), saveGraphPath = "./TLSAnalysis"){

  if(is.null(standName)){standName = "TestStd"}
  dir.create(file.path(saveGraphPath), showWarnings = FALSE)
  dir.create(file.path(saveGraphPath, standName), showWarnings = FALSE)


  for(i in Threshold){
    nom<-paste(saveGraphPath,"/",standName, "/", i,"_", standName,"_ThresoldTest.jpg",sep="")
    jpeg(nom,width = 4000, height = 2500)
    suppressWarnings(print(
      ggplot(data[selec1==0 & alph>i])+
      geom_point(aes(xa,ya),size=0.8)+
      coord_fixed()+
      ggtitle(paste("LiDAR: Threshold = ", i, sep=""))+
      theme_bw(base_size=35)
    ))

    dev.off()
  }

  message("View the saved graphs in the selected folder to determine the best value of alpha (needed for the TreeMap)")
}