R/TreeMap.R

#' Automatic creation of a TreeMap
#'
#' Detect the coordinates XY and the Diameter at breast height (DBH) of all trees on the plot.
#' The point cloud had to be prepared and sliced with the function 'DataInit' and a value of threshold had to be found (use testAlphaThreshold())
#'
#' @param Mydata 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 alpha Alpha is the threshold value needed to clean the point cloud data
#' @param minPts Minimum distance between 2 points to belong to the same cluster
#' @param eps Minimal spacing between 2 points to create 2 clusters
#' @param minDbh Minimum DBH value saved (if the value is greater than 0, trees with a lower DBH will be removed from the TreeMap)
#' @param savePlot Saves a map of the position and DBH of all trees (This map will be used during the manual step to possibly detect adjustment errors or missing trees)
#' @param TreeValue Used to add particular trees to the map
#' @param savePath indicates the folder where the map and the data will be saved
#'
#'
#' @return Generate a .txt file containing the position and the DBH of all the trees as well as a map (Treedetection.jpg) positioning them on the plot
#' @export
#'
#' @importFrom dbscan dbscan
#' @importFrom RANN nn2
#' @importFrom rgl plot3d
#' @importFrom rgl aspect3d
TreeMap_Auto<-function(Mydata, standName=NULL, alpha, minPts = 8, eps= 0.1,minDbh=0,savePlot=T,TreeValue=NULL, savePath = "./TLSAnalysis/"){

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

  # minDbh in mm

  # # for debbugging
  # Mydata = Mydata; alpha=5; minPts = 8; eps= 0.1; minDbh = 0

#cluster
  vox_placette<-Mydata[selec1 %in% c(0,-5,5) & alph>alpha]
  keep_cols = c("X", "Y")
  vox_placette<- vox_placette[, ..keep_cols] #only keep the columns X Y
  db_pl1=dbscan::dbscan(vox_placette,minPts = minPts,eps=eps)


  #Tree map#
  tree_map_R=NULL
  for (i in   unique(db_pl1$cluster)[unique(db_pl1$cluster)>0]) {

    test=vox_placette[db_pl1$cluster==i,]
    # test 5 methods of circle fitting to calculate the DBH
    temp<-list(CircleFitByKasa(as.matrix(test)),CircleFitByPratt(as.matrix(test)),CircleFitByLandau(as.matrix(test)),
               CircleFitBySpath(as.matrix(test)),CircleFitByTaubin(as.matrix(test)))
    Criterion<-NULL
    for(j in 1:5){
      pts_circle=calculateCircle(temp[[j]][1],temp[[j]][2],temp[[j]][3])
      dist_Circle_pts=apply(test,1,function(x){
        min(sqrt((x[1]-pts_circle[,1])^2+(x[2]-pts_circle[,2])^2))})
      essai=sum(dist_Circle_pts)
      Criterion<-c(Criterion,essai)
    }
    #the criterion variable detect the fitting accuracy... The lower, the better
    k<-which(Criterion==min(Criterion))
    crit<-unique(round(Criterion[k]/(temp[[k]][,3]*nrow(test)),4))

    # if(fit_circle[,3]*2000>0&fit_circle[,3]*2000<10000) {
    pts_circle=calculateCircle(temp[[k]][1],temp[[k]][2],temp[[k]][3])
    #points(pts_circle,col="red",type="l",lwd=1)
    #text(temp[[k]][1],temp[[k]][2],paste("a.",i,sep=""),cex=1,pos=3)
    tree_map_R=rbind(tree_map_R,c(id=paste("a.",i,sep=""),X=temp[[k]][1],Y=temp[[k]][2],DBH=temp[[k]][3]*2000,Crit=crit,Num=i))
    #}
    # print(i)
  }
  tree_map_R=data.table(tree_map_R)



  ### 2nd round. An other fitting is done by selecting the center of the tree detected on the previous loops
  ### Only the closest points to the center are kept to limit the noise around trees
  tree_map_R1=NULL
  tree_map_R$Num<-as.numeric(paste(tree_map_R$Num))
  TreeDetect<-unique(tree_map_R$Num)
  for(k1 in TreeDetect){
    test=vox_placette[db_pl1$cluster==k1]
    Treepos<-data.table(xc=as.numeric(tree_map_R[Num==k1]$X),
                        yc=as.numeric(tree_map_R[Num==k1]$Y))

    n1<-dim(test)[1]
    test$dist<-c(nn2(test, query = Treepos, k = n1)$nn.dists)

    test$group<-1
    try(test[test$ya<=yc]$group<-2,silent=T)

    test<-test%>%
      mutate(xa = round(X,3),
             ya = round(Y,3))

    test<-test %>%
      arrange(group,xa,dist)%>%
      mutate(selec=seq(1,n(),1))

    test<-setDT(test)
    val<-mean(test[selec==1]$dist)+sd(test[selec==1]$dist)/4
    test$selec2<-0
    test[test$dist<=val]$selec2<-1
    test1<-test[selec==1 & selec2==1,1:2]
    #test1<-test[test$selec==1,1:2]

    if(length(test1$X)<=3)
    {test1<-test[,1:2]}

    temp<-list(CircleFitByKasa(as.matrix(test1)),CircleFitByPratt(as.matrix(test1)),CircleFitByLandau(as.matrix(test1)),
               CircleFitBySpath(as.matrix(test1)),CircleFitByTaubin(as.matrix(test1)))
    Criterion<-NULL
    for(j in 1:5){
      pts_circle=calculateCircle(temp[[j]][1],temp[[j]][2],temp[[j]][3])
      dist_Circle_pts=apply(test1,1,function(x){
        min(sqrt((x[1]-pts_circle[,1])^2+(x[2]-pts_circle[,2])^2))})
      essai=sum(dist_Circle_pts)
      Criterion<-c(Criterion,essai)
    }
    k<-min(which(Criterion==min(Criterion)) )
    crit<-unique(round(Criterion[k]/(temp[[k]][,3]*nrow(test1)),4))
    tree_map_R1=rbind(tree_map_R1,c(id=paste("a.",k1,sep=""),X=temp[[k]][1],Y=temp[[k]][2],DBH=temp[[k]][3]*2000,Crit=crit,Num=k1))

  }
  tree_map_R1=data.table(tree_map_R1)

  #select the best between tree_map_R and Tree_map_R1
  tree_map_R$Num<-as.numeric(paste(tree_map_R$Num))
  tree_map_R$DBH<-as.numeric(paste(tree_map_R$DBH))
  tree_map_R$X<-as.numeric(paste(tree_map_R$X))
  tree_map_R$Y<-as.numeric(paste(tree_map_R$Y))
  tree_map_R$Crit<-as.numeric(paste(tree_map_R$Crit))

  tree_map_R1$Num<-as.numeric(paste(tree_map_R1$Num))
  tree_map_R1$DBH<-as.numeric(paste(tree_map_R1$DBH))
  tree_map_R1$X<-as.numeric(paste(tree_map_R1$X))
  tree_map_R1$Y<-as.numeric(paste(tree_map_R1$Y))
  tree_map_R1$Crit<-as.numeric(paste(tree_map_R1$Crit))

  tree_map_R2=NULL
  for(k2 in unique(tree_map_R$Num)){
    temp1<-rbind(tree_map_R[tree_map_R$Num==k2,],tree_map_R1[tree_map_R1$Num==k2,])
    temp1<- temp1 %>%
      arrange(Crit) %>%
      mutate(sel = seq(1,n(),1))
    test1<-temp1[temp1$sel==1,]
    tree_map_R2=rbind(tree_map_R2,test1)
  }
  tree_map_R2=data.table(tree_map_R2)


  ###
  #some trees may overlap. in this case, only one of the two
  tree_map_R2$Num<-as.numeric(paste(tree_map_R2$Num))
  tree_map_R2$X<-as.numeric(paste(tree_map_R2$X))
  tree_map_R2$Y<-as.numeric(paste(tree_map_R2$Y))
  tree_map_R2$DBH<-as.numeric(paste(tree_map_R2$DBH))
  tree_map_R2$recoup<-NA;tree_map_R2$dist<-NA
  tree_map_R2$recoup<-as.numeric(tree_map_R2$recoup)
  for(t1 in tree_map_R2$Num){
    Treepos<-data.table(xt1=tree_map_R2[Num==t1]$X,
                        yt1=tree_map_R2[Num==t1]$Y)

    dbht1<-tree_map_R2[Num==t1]$DBH/1000

    n1<-dim(tree_map_R2)[1]
    tree_map_R2$dist<-c(nn2(tree_map_R2[,2:3], query = Treepos, k = n1)$nn.dists)
    tree_map_R2$dist<-tree_map_R2$dist - dbht1
    tree_map_R2[Num==t1]$dist<-9999

    if(length(tree_map_R2[!is.na(tree_map_R2$dist) & tree_map_R2$dist<0 ,]$recoup)>0)
    {tree_map_R2[!is.na(tree_map_R2$dist) & tree_map_R2$dist<0 ,]$recoup<-t1}
  }
  #we select duplicates
  tree_map1<-tree_map_R2[!is.na(recoup)]
  if(dim(tree_map1)[1]>0)
  {

    for(t3 in 1:length(tree_map1$recoup)){
      temp<-sort(c(as.numeric(paste(tree_map1$Num[t3])),as.numeric(tree_map1$recoup[t3])))
      tree_map1$merg1[t3]<-paste(temp[1],"_",temp[2],sep="")
    }
    tree_map1  <- tree_map1 %>%
      group_by(merg1) %>%
      arrange(-DBH) %>%
      mutate(selec=seq(1,n(),1))
    tree_map1<-tree_map1[tree_map1$selec==1,]

    treesuppr<-as.numeric(tree_map1$recoup)
    tree_map_R2$selec1<-NA
    tree_map_R2$selec1<-as.numeric(tree_map_R2$selec1)
    for(i in treesuppr){
      tree_map_R2[Num==i]$selec1<-1
    }
    tree_map_R2<-tree_map_R2[is.na(tree_map_R2$selec1),1:6]
    #we put a high criterion to the remaining tree to check is in the manual selection
    for(i in unique(tree_map1$Num)){
      if(length(tree_map_R2[tree_map_R2$Num==i,]$Crit)>0)
      {tree_map_R2[tree_map_R2$Num==i,]$Crit<-10}
    }

  }
  # tree_map_R2$recoup<-NULL
  # tree_map_R2$dist<-NULL
  # tree_map_R2$sel<-NULL
  # filter trees with minimum dbh
  tree_map_R2<-tree_map_R2[DBH > minDbh]
  #tree_map_R2=data.frame(tree_map_R2)
  nom3<-paste(savePath,standName, "/",standName,"_tree_map.txt",sep="")
  write.table(tree_map_R2,nom3,row.names=F,dec=".")

  #plot TreeMap
  if (savePlot==TRUE){
    nom<-paste(savePath,standName, "/",standName,"_Tree_detection.jpg",sep="")
    jpeg(nom,width = 10000, height = 8000)

    plot(Mydata[selec1 %in% c(0,-5,5) & alph>(alpha/2)],cex=1.6,pch=20,asp=1, cex.axis=5, panel.first = grid())

    if(!is.null(TreeValue)){
      points(TreeValue$Xmodif, TreeValue$Ymodif, pch=21, bg="blue",cex=5,lwd=3)
      text(TreeValue$Xmodif, TreeValue$Ymodif, TreeValue$nar,cex=3,pos=3, col="blue")
    }

    for(k2 in unique(tree_map_R2$Num)){
      pts_circle=calculateCircle(tree_map_R2[tree_map_R2$Num==k2,]$X,tree_map_R2[tree_map_R2$Num==k2,]$Y,tree_map_R2[tree_map_R2$Num==k2,]$DBH/2000)
      points(pts_circle,col="red",type="l",lwd=4)
      text(tree_map_R2[tree_map_R2$Num==k2,]$X,tree_map_R2[tree_map_R2$Num==k2,]$Y,paste(tree_map_R2[tree_map_R2$Num==k2,]$id),cex=3,pos=3, col="green")
    }
    dev.off()
  }

  message("TreeMap automatic detection ... OK")
}




#' Semi-automatic TreeMap corrections
#'
#'From the TreeMap created with the function 'TreeMap_Auto', all trees having a too high value
#'of the variable 'critR' will be considered as potentially problematic and will be displayed to be manually corrected.
#'To correct the displayed tree, please click on the center of the tree, then click on Esc.
#'It is possible to click up to 5 times if necessary (i.e. 5 trees need to be positionned on the same image).
#'In this case, click on Esc only at the end of the multiple selection.
#'
#'
#' @param Mydata 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).
#' WARNING! If you work with multiple plots and do not specify a unique name for each new plot, the default name will be used each time and the TreeMap of the previous plot will be overwritten.
#' @param alpha Alpha is the threshold value needed to clean the point cloud data
#' @param critR Indicates the threshold value of the variable 'critR' for which a tree will be displayed and will require a correction
#' @param Dist_view Indicates the distance displayed around the target tree (Dist_view x DBH)
#' @param add3Dhelp if TRUE, a 3d plot of the tree will be displyed to help correction
#' @param savePath indicates the folder where the map and the data will be saved
#' @param TreeValue Used to add particular trees to the map
#'
#' @return Generate a .txt file containing the position and the DBH of all the trees as well as a map (Treedetection.jpg) positioning them on the plot
#' @export
#'
TreeMap_semiManual<-function(Mydata, standName=NULL, alpha, critR = 0.1,
                             Dist_view = 2,add3Dhelp=F, savePath = "./TLSAnalysis/",TreeValue=NULL){

  message("To correct the displayed tree, please click on the center of the tree, then click on Esc. it is possible to click up to 5 trees if necessary (in this case, click on Esc only at the end of the multiple selection)")

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

  nom<-paste(savePath,standName, "/",standName,"_tree_map.txt",sep="")
  tree_map_R2<-fread(nom)

  #map_R3 contains trees without problems, R4 those potentially problematic
  tree_map_R3<-tree_map_R2[tree_map_R2$Crit<critR,]
  tree_map_R4<-tree_map_R2[tree_map_R2$Crit>=critR ,]
  if(dim(tree_map_R4)[1] !=0){
    tree_map_R5<-NULL
    sel<-as.numeric(tree_map_R4$Num)
    for(i in sel){

      xmin2<-tree_map_R2[Num==i]$X - (tree_map_R2[Num==i]$DBH*Dist_view/1000)
      xmax2<-tree_map_R2[Num==i]$X + (tree_map_R2[Num==i]$DBH*Dist_view/1000)
      ymin2<-tree_map_R2[Num==i]$Y - (tree_map_R2[Num==i]$DBH*Dist_view/1000)
      ymax2<-tree_map_R2[Num==i]$Y + (tree_map_R2[Num==i]$DBH*Dist_view/1000)

      test1<-Mydata[X>=xmin2 & X<=xmax2 & Y>=ymin2 & Y<=ymax2  & alph>alpha/2 & selec1>-4 & selec1<4]
      #test0<-Mydata[Mydata$X>=xmin5 & Mydata$X<=xmax5 & Mydata$Y>=ymin5 & Mydata$Y<=ymax5  & Mydata$alph>alpha/2 ,]
      test<-Mydata[X>=xmin2 & X<=xmax2 & Y>=ymin2 & Y<=ymax2 & alph>alpha/2 ,]

      pts_circle=calculateCircle(tree_map_R2[Num==i]$X,tree_map_R2[Num==i]$Y,(tree_map_R2[Num==i]$DBH/2000))
      if(isTRUE(add3Dhelp)){rgl::plot3d(test$X, test$Y, test$Z ,size=3);aspect3d("iso")}

      plot(test$X, test$Y,asp=1,pch=20, col="grey",main=paste("Crit: ",tree_map_R4[Num==i]$Crit, "  Id: ",tree_map_R4[Num==i]$id, sep=""))
      points(test1$X, test1$Y,pch=20)
      #points(pts_circle[,1],pts_circle[,2] ,col="red",lwd=2,type="l")
      lines(pts_circle, col="red",lwd=3)
      for(comp in tree_map_R2[Num!=i]$Num){
        pts_circle=calculateCircle(tree_map_R2[Num==comp]$X,tree_map_R2[Num==comp]$Y,(tree_map_R2[Num==comp]$DBH/2000))
        lines(pts_circle, col="orange", lwd=3)
      }

      #we click manually on the center of the tree stem to correct with locator. 5 times maximum
      ps<-as.data.frame(locator(5))

      ####If 0 trees (only branches), clicks on "Echap"
      if(length(ps)==0)
      {tree_map_R5=rbind(tree_map_R5,c(id=0,X=0,Y=0,DBH=0,Crit=0,Num=i))}

      ###If at least one tree to correct...
      if(dim(ps)[1]>=1)
      {
        #
        nloc<-dim(ps)[1]
        for(p1 in 1:nloc){

          test1<-test1%>%
            mutate(dist=sqrt((X-ps[p1,1])^2+ (Y-ps[p1,2])^2))
          #we will take the average distance of the nearest 10% of the point where we clicked
          #to avoid some problems with points in the middle of the tree that reduces the diameter
          pct<-round(0.1*length(test1$X))
          val<-mean(sort(test1$dist)[1:pct])
          n1<-length(test1[test1$dist< (1.5*val),]$X)
          if(n1<4)
          {test2<-test1[test1$dist< (2*val),1:2]}
          if(n1>=4)
          {test2<-test1[test1$dist< (1.5*val),1:2]}

          temp<-list(CircleFitByKasa(as.matrix(test2)),CircleFitByPratt(as.matrix(test2)),CircleFitByLandau(as.matrix(test2)),
                     CircleFitBySpath(as.matrix(test2)),CircleFitByTaubin(as.matrix(test2)))
          Criterion<-NULL
          for(j in 1:5){
            pts_circle=calculateCircle(temp[[j]][1],temp[[j]][2],temp[[j]][3])
            dist_Circle_pts=apply(test2,1,function(x){
              min(sqrt((x[1]-pts_circle[,1])^2+(x[2]-pts_circle[,2])^2))})
            essai=sum(dist_Circle_pts)
            Criterion<-c(Criterion,essai)
          }
          k<-which(Criterion==min(Criterion))
          crit<-unique(round(Criterion[k]/(temp[[k]][,3]*nrow(test)),4))
          pts_circle=calculateCircle(temp[[k]][1],temp[[k]][2],temp[[k]][3])
          points(test2[,1:2],pch=20,col="blue")
          lines(pts_circle,col="green",type="l",lwd=2)
          tree_map_R5=rbind(tree_map_R5,c(id=paste("a.",i,"_",p1,sep=""),X=temp[[k]][1],Y=temp[[k]][2],DBH=temp[[k]][3]*2000,Crit=crit,Num=paste(i,"_",p1,sep="")))
        }

      }

    }

    tree_map_R5=data.table(tree_map_R5)
    tree_map_R5$Num<-as.factor(paste(tree_map_R5$Num))
    tree_map_R5$DBH<-as.numeric(paste(tree_map_R5$DBH))
    tree_map_R5$X<-as.numeric(paste(tree_map_R5$X))
    tree_map_R5$Y<-as.numeric(paste(tree_map_R5$Y))

    tree_map_R5<-tree_map_R5[tree_map_R5$X!=0,]
    #tree_map_R3$sel<-NULL
    tree_map<-rbind(tree_map_R3,tree_map_R5)
    tree_map<-tree_map
    nom3<-paste(savePath,standName, "/",standName,"_tree_map.txt",sep="")
    write.table(tree_map,nom3,row.names=F,dec=".")

    #tree_map<-read.table("tree_map.txt",header=T)



    vox_placette<-Mydata[selec1 %in% c(0,-4,4) & alph>alpha/2]
    keep_cols = c("X", "Y")
    vox_placette<- vox_placette[, ..keep_cols]

     nom<-paste(savePath,standName, "/",standName,"_Tree_detection.jpg",sep="")
    jpeg(nom,width = 10000, height = 8000)
    plot(Mydata[selec1 %in% c(0,-5,5) & alph>(alpha/2)],cex=1,pch=20,asp=1, panel.first = grid())

    if(!is.null(TreeValue)){
      points(TreeValue$Xmodif, TreeValue$Ymodif, pch=21, bg="blue",cex=5,lwd=3)
      text(TreeValue$Xmodif, TreeValue$Ymodif, TreeValue$nar,cex=3,pos=3, col="blue")
    }
    for(k2 in unique(tree_map$Num)){
      pts_circle=calculateCircle(tree_map[Num==k2]$X,tree_map[Num==k2]$Y,tree_map[Num==k2]$DBH/2000)
      points(pts_circle,col="red",type="l",lwd=4)
      text(tree_map[tree_map$Num==k2,]$X,tree_map[Num==k2]$Y,paste(tree_map[Num==k2]$id),cex=3,pos=3, col="green")
    }

    dev.off()
    message("All trees with a too high 'critR' parameter have been fixed.")
    message("However, visual validation of the tree map should be done to determine if manual corrections are required (missing or poorly adjusted trees).")
  }

  if(dim(tree_map_R4)[1] ==0) {
    message("According to the parameter 'critR' no trees need to be corrected.")
    message("However, visual validation of the tree map should be done to determine if manual corrections are required (missing or poorly adjusted trees).")
  }
}
EmmanuelDuchateau/TLS documentation built on May 14, 2019, 2:01 p.m.