#**************************************************************************** # # Copyright (C) 2017-2020 Dr. Jan Hackenberg, free software developer # All rights reserved. # # Contact : https://gitlab.com/simpleForest # # Developers : Dr. Jan Hackenberg # # This file is developed for SimpleForest plugin Version 1 for Computree. # # SimpleForest plugin is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # SimpleForest plugin is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with SimpleForest plugin. If not, see . # # SimpleForest is an extended version of the SimpleTree platform. # # *****************************************************************************/ # If you use this script, please cite: # SimpleTree —An Efficient Open Source Tool to Build Tree Models from TLS Clouds # I work on a more updated peer review paper and as soon I have a dedicated # SimpleForest publication please rather cite the updated one. Not done yet though # *****************************************************************************/ library(nls2) setwd('/home/drsnuggles/clouds/evaluation/sphereGV/qsm/detailed/') file_name = 'Q12.csv' #################################################################################################### #total volume and total length per branch order ################### #################################################################################################### branch_order_histogram <- function(data){ maxOrder = max(data$branchOrder) order = data.frame(order=as.numeric(), sum_length=as.numeric(),sum_volume=as.numeric()) for(i in 1:maxOrder){ order = rbind(order,c(0,0,0)) } for(i in 1:maxOrder){ data_order = data[which(data$branchOrder==i),] order[i,1] = i; order[i,2] = sum(data_order$length); order[i,3] = sum(data_order$volume); } names(order) = c("order","length", "volume") order } median_radius<-function(segment) { median(segment$radius) } #################################################################################################### #creates a matrix with an entry for every segment ################### #the average segments area is one row entry, the sum of its child segments areas the other########## #################################################################################################### area_parent_children <- function (data, binWidth){ number_segments <- max(data$segmentID) area = data.frame(area_parent=as.numeric(),area_child=as.numeric()) for(i in 1:number_segments) { area = rbind(area,c(0,0)) } for(i in 1:number_segments) { data_parent = data[which(data$segmentID==i),] radius = median_radius(data_parent) area[i,1] = radius*radius*3.1416349*10000 data_children = data[which(data$parentSegmentID==i),] segmentID_children = unique(data_children$segmentID) if(!is.null(segmentID_children)) for(j in segmentID_children) { data_child = data[which(data$segmentID==j),] radius = median_radius(data_child) area_child = radius*radius*3.1416349*10000 area[i,2] = area[i,2] + area_child } } names(area) = c("parent_area","child_area") area <- area[which(area$child_area> 0),] area } #################################################################################################### #counts for every stem height bin the total volume of branches growing out of ################### #this segment ######################## #################################################################################################### branch_histogram <- function (data, binWidth){ min <- min(data$startZ) data$startZ <- data$startZ - min data$endZ <- data$endZ - min max_height <- max(data$endZ) branches = data.frame(height_low=as.numeric(), height_high=as.numeric(),sum_vol_branches=as.numeric()) height = 0.0; #################################################################################################### #generate an empty frame with height bins of the stem to store the branch volume ################### #per bin ######################## #################################################################################################### while(height < max_height) { branches = rbind(branches, c(height, height + binWidth, 0 )) height <- height + binWidth; } #################################################################################################### # we get the number of branches ########################################## #################################################################################################### number_branches = max(data$branchID) branchID = unique(data$branchID) branchID = branchID[which(branchID>0)] #################################################################################################### # every stem cylinders volume and radius is added to the ########################################## # height bin containing the cylinders middle height ########################################## # a counter is added for each bin counting the sum of added cylinders ############################## #################################################################################################### for(i in branchID) { data_branch = data[which(data$branchID==i),] height_branch = data_branch$startZ[1] index <- ceiling(height_branch/binWidth) branches[index,3] <-branches[index,3] + sum(data_branch$volume) } names(branches) = c("branch_min_height", "branch_max_height", "branch_volume") branches } #################################################################################################### #Diameter (and Volume) distribution of the stem along the stem height ################### #################################################################################################### stemCurve <- function (data, binWidth){ min <- min(data$startZ) data$startZ <- data$startZ - min data$endZ <- data$endZ - min max_height <- max(data$endZ) stem_curve = data.frame(height_low=as.numeric(), height_high=as.numeric(),sum_radius=as.numeric(), sum_cylinders=as.numeric(), radius=as.numeric(), diameter=as.numeric(), volume=as.numeric()) height = 0.0; #################################################################################################### #generate an empty frame, sum_radius and sum_cylinders are helper variables ######################## #the empty frame contains one row for every height bin ######################## #################################################################################################### while(height < max_height) { stem_curve = rbind(stem_curve, c(height, height + binWidth, 0, 0,0,0,0)) height <- height + binWidth; } #################################################################################################### # we create a helper frame containing only stem cylinders ########################################## #################################################################################################### data.stem = data[which(data$branchID==0),] #################################################################################################### # every stem cylinders volume and radius is added to the ########################################## # height bin containing the cylinders middle height ########################################## # a counter is added for each bin counting the sum of added cylinders ############################## #################################################################################################### for(i in 1:length(data.stem[,1])) { height_cylinder <- (data.stem$startZ[i] + data.stem$endZ[i])/2 index <- ceiling(height_cylinder/binWidth) stem_curve[index,3] <- stem_curve[index,3] + data.stem$radius[i]; stem_curve[index,4] <- stem_curve[index,4] + 1; stem_curve[index,7] <- stem_curve[index,7] + data.stem$volume[i]; } #################################################################################################### # we remove empty bins #### #################################################################################################### stem_curve <- stem_curve[which(stem_curve[,4]!=0),] #################################################################################################### # mean radius is computed for each bin, radius [m] is converted to diameter [cm] #### #################################################################################################### for(i in 1:length(stem_curve[,1])) { stem_curve[i,5] <- stem_curve[i,3]/stem_curve[i,4]; stem_curve[i,6] <- 200 * stem_curve[i,5]; } stem_curve <- data.frame(stem_curve[,1],stem_curve[,2],stem_curve[,6],stem_curve[,7]) names(stem_curve) = c("min_height","max_height", "diameter", "volume") stem_curve } #################################################################################################### #Allometric model from one tree ################### #################################################################################################### allom <- function (data){ # Only use unmodified radii data data.sub <- data[which(as.character(data$FittingType) == " SPHEREFOLLOWING"),] x<-data$radius y<-data$growthVolume lm <- lm(y.log~x.log) a <- coef(lm)[1] b <- coef(lm)[2] expA <- exp(a) data2 <- data.frame(x,y) nls.vol <- nls(y ~ b1*x**b2,data = data2,start = list(b1 = expA,b2 = b), control = list(maxiter=500,warnOnly = TRUE)) nls.vol } data <- read.csv(file_name, sep=",") binWidth = 0.5; order = branch_order_histogram(data) plot(order$order, order$length, xaxt="n", xlab = "Branch Order", ylab = "Accumulated Branch Length [m]", main = "Branch Order Histogram Length", pch = 20) axis(1, at = order$order, las=2) plot(order$order, order$volume, xaxt="n", xlab = "Branch Order", ylab = "Accumulated Branch Volume [m^3]", main = "Branch Order Histogram Volume", pch = 20) axis(1, at = order$order, las=2) branches = branch_histogram(data,binWidth) max(branches$branch_volume) plot(branches$branch_min_height, branches$branch_volume, xlab = "Branch Height [m]", ylab = "Accumulated Branch Volume [m^3]", main = "Branch Histogram", pch = 20) stem = stemCurve(data,binWidth) plot(stem$min_height, stem$diameter, xlab = "Height [m]", ylab = "Diameter [cm]", main = "Stem Curve", pch = 20) area = area_parent_children(data) plot(area$parent_area, area$child_area, xlab = "Cross sectional area before branch junction", ylab = "Summed Cross sectional areas after branch junction", main = "Pipe model theory", pch = 20) lm.area <- lm(area$child_area ~ area$parent_area) abline(lm.area, lwd = 2) max<- max(max(area$parent_area),max(area$child_area)) lines(c(0,max),c(0,max), lwd = 2, lty = 2) summary(lm.area) leg.txt = c("linear model", "1:1 line") legend("left", leg.txt ,lwd = c(2) ,lty = c(1,2) ,pt.cex =1 ) ## Be Careful here and watch for: ## * Its radius and not diameter ## * We use meter here allomModel <- allom(data) plot(data$radius, data$growthVolume, xlab = "Radius [m]", ylab = "Volume [m^3]", main = "Tree specific allometry", sub = "Please watch units", pch = 20) r2nls <- 1-(deviance(allomModel)/sum((data$growthVolume-mean(data$growthVolume))^2)) curve( (coefficients(allomModel)[1] *x** coefficients(allomModel)[2]),0,max(x),add=TRUE,lwd=2) r2 <- paste("r-squared :", format(r2nls, digits = 2)) a <- paste("coefficient a :", format(coefficients(allomModel)[1], digits = 2)) b <- paste("coefficient b :", format(coefficients(allomModel)[2], digits = 2)) leg.txt = c(r2, a, b) legend("topleft", leg.txt ,col = c("white") ,pt.cex =1 )