Commit 725b315e authored by abc's avatar abc

corrected bugs in cloudTreeMetrics

parent 835bc8a7
......@@ -13,6 +13,7 @@ export(cimg2Raster)
export(circle2Raster)
export(cleanPredictionRaster)
export(cloudMetrics)
export(cloudTreeMetrics)
export(coregistration)
export(createDisk)
export(demFiltering)
......
......@@ -372,8 +372,8 @@ plotTreeInventory <-function(xy, height=NULL, diam=NULL, species=NULL, add=FALSE
#' raster::res(r) <- 1
#'
#' # xy positions
#' xy <- data.frame(c(10,20,31.25,15),
#' c(10,20,31.25,25))
#' xy <- data.frame(x=c(10,20,31.25,15),
#' y=c(10,20,31.25,25))
#' # compute mask
#' mask1 <- rasterXYMask(xy, c(5, 8, 5, 5), r)
#' mask2 <- rasterXYMask(xy, c(5, 8, 5, 5), r, binary=FALSE)
......@@ -390,7 +390,7 @@ plotTreeInventory <-function(xy, height=NULL, diam=NULL, species=NULL, add=FALSE
rasterXYMask <- function(xy, buff, r, binary=TRUE)
{
# convert vector to data.frame in case only one pair of coodinates is provided
if (length(xy) == 2) {xy <- data.frame(matrix(coord, 1, 2))}
if (is.null(dim(xy)) & length(xy)==2) {xy <- data.frame(matrix(xy, 1, 2))}
# compute squared buffers
buff2 <- buff^2
# compute XY coordinates of cell centers
......
......@@ -155,7 +155,7 @@ ABAmodelMetrics <- function(z, i, rn, c, hmin = 2, breaksH=NULL)
#'
stdTreeMetrics <- function(x, area.ha=NA)
{
data.frame(Tree.meanH=mean(x$h), Tree.sdH=stats::sd(x$h), Tree.giniH=reldist::gini(x$h), Tree.density=length(x$h)/area.ha, TreeInf10.density=sum(x$h<=10)/area.ha, TreeSup10.density=sum(x$h>10)/area.ha, TreeSup20.density=sum(x$h>20)/area.ha, TreeSup30.density=sum(x$h>30)/area.ha, Tree.meanCrownSurface=mean(x$s), Tree.CanopyCover=sum(x$sp)/area.ha/10000, Tree.meanCrownV=mean(x$v), Tree.meanCanopyH=sum(x$vp)/area.ha/10000, Tree.meanCrownH=sum(x$v)/sum(x$s))
data.frame(Tree.meanH=mean(x$h), Tree.sdH=stats::sd(x$h), Tree.giniH=ifelse(is.null(x), NA, reldist::gini(x$h)), Tree.density=length(x$h)/area.ha, TreeInf10.density=sum(x$h<=10)/area.ha, TreeSup10.density=sum(x$h>10)/area.ha, TreeSup20.density=sum(x$h>20)/area.ha, TreeSup30.density=sum(x$h>30)/area.ha, Tree.meanCrownSurface=mean(x$s), Tree.CanopyCover=sum(x$sp)/area.ha/10000, Tree.meanCrownV=mean(x$v), Tree.meanCanopyH=sum(x$vp)/area.ha/10000, Tree.meanCrownH=sum(x$v)/sum(x$s))
}
#
###############################
......@@ -219,24 +219,50 @@ points2terrainStats <- function(p, centre=NULL, r=NULL)
round(data.frame(altitude=altitude, azimut.gr=azimut, slope.gr=slope, adjR2.plane=summary(modlin)$adj.r.squared * 100),1)
} else { NULL}
}
################################
#' compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest). Call the treedetection function on each object and then arranges the results in a data.frame
#' Computation of tree metrics
#'
#' This function performs tree detection in each of LAS objects contained in a list and then computes summary statistics based on compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest).
#' - calls the tree detection and extraction functions on each object
#' - remove the detected trees located outside of the regions of interest defined by their centers and radius
#' - computes summary statistics of remaining detected trees.
#'
#' @param llasn list of LAS objects (from lidR package)
#' @param XY a dataframe or matrix with XY coordinates of plot centers
#' @param plot.radius a numeric plot radius in m
#' @param res CHM resolution for tree detection (m)
#' @param hmin min height for a tree (m)
#' @param h.points height for CHM thresholding (m)
#' @param plot.radius numeric. plot radius in same unit as coordinates of LAS objects
#' @param res numeric. resolution of canopy height model computed with \code{\link{points2DSM}} before tree segmentation
#' @param func a function to be applied to the attributes of extracted trees (return from internal call to treeExtraction function to compute plot level metrics
#' @param ... other parameters to be passed to \code{\link{treeSegmentation}}
#' @return a dataframe with tree metrics in columns corresponding to LAS objects of the list (lines)
#' @seealso \code{\link{treeSegmentation}}, \code{\link{treeExtraction}}, \code{\link{stdTreeMetrics}}
#' @examples
#' data(laschablais3)
#'
#' # extract three point clouds from LAS object
#' llas <- list()
#' llas[["A"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581680, 10)
#' llas[["B"]] <- lidR::lasclipCircle(laschablais3, 974390, 6581680, 10)
#' llas[["C"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581640, 10)
#' # normalize point clouds
#' llas <- lapply(llas, function(x) {lidR::lasnormalize(x, lidR::tin())})
#'
#' # compute tree metrics
#' cloudTreeMetrics(llas,
#' cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
#' 10, res=0.5)
#'
#' # compute metrics with user-defined function
#' # number of detected trees between 20 and 30 meters and their mean height
#' user.func <- function(x, area.ha=NA)
#' {
#' dummy <- x$h[which(x$h>10 & x$h<20)]
#' data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
#' }
#' cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
#' @export
#'
cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, hmin=5, h.points=60, func=stdTreeMetrics)
cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, func=stdTreeMetrics, ...)
{
# ADD example
# ALLOW ... parameters to be passed to treeSegmentation
plot.area.ha <- pi*plot.radius^2/10000
xy.list <- split(XY, seq(nrow(XY)))
# LOOP on list <- replace by lapply
......@@ -248,16 +274,18 @@ cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, hmin=5, h.points=6
# compute dsm
dummy <- points2DSM(x,res=res)
# replace NA, low and high values
dummy[is.na(dummy) | dummy<0 | dummy>h.points] <- 0
dummy[is.na(dummy) | dummy<0] <- 0
# tree detection
dummy <- treeSegmentation(dummy, hmin=hmin)
dummy <- treeSegmentation(dummy, ...)
# tree extraction
lsegms[[i]] <- treeExtraction(dummy$filled.dem, dummy$local.maxima, dummy$segments.id, rasterXYMask(coord, plot.radius, dummy$local.maxima, binary=T))
dummy <- treeExtraction(dummy$filled.dem, dummy$local.maxima, dummy$segments.id, rasterXYMask(coord, plot.radius, dummy$local.maxima, binary=T))
lsegms[[names(llasn)[i]]] <- dummy
}
# GERER les cas ou pas d'arbre detecte ??
# GERER les NA
# compute metrics
tree.metrics <- lapply(lsegms, FUN=function(x){func(x, area.ha=plot.area.ha)})
# bind data.frames
tree.metrics <- do.call(rbind,tree.metrics)
}
#
if (length(lsegms)>0)
{
tree.metrics <- lapply(lsegms, FUN=function(x){func(x, area.ha=plot.area.ha)})
tree.metrics <- do.call(rbind,tree.metrics)
} else {tree.metrics <- NULL}
return(tree.metrics)
}
\ No newline at end of file
......@@ -541,6 +541,11 @@ treeExtraction <- function(r.dem.nl, r.maxi, r.dem.w, r.mask=NULL)
#
# extract segment id, height and dom.radius
cells <- which(raster::values(r.maxi)>0)
# check if positive values are present
if (length(cells)==0)
{
segms <- NULL
} else {
coord <- raster::xyFromCell(r.maxi,cells)
segms <- data.frame(coord,id=raster::values(r.dem.w)[cells],h=raster::values(r.dem.nl)[cells],dom.radius=raster::values(r.maxi)[cells])
segms <- merge(segms, s, all.x = T)
......@@ -552,6 +557,7 @@ treeExtraction <- function(r.dem.nl, r.maxi, r.dem.w, r.mask=NULL)
}
sp::coordinates(segms) <- ~ x+y
segms@proj4string <- r.dem.nl@crs
}
segms
}
#
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/metrics.R
\name{cloudTreeMetrics}
\alias{cloudTreeMetrics}
\title{Computation of tree metrics}
\usage{
cloudTreeMetrics(llasn, XY, plot.radius, res = 0.5,
func = stdTreeMetrics, ...)
}
\arguments{
\item{llasn}{list of LAS objects (from lidR package)}
\item{XY}{a dataframe or matrix with XY coordinates of plot centers}
\item{plot.radius}{numeric. plot radius in same unit as coordinates of LAS objects}
\item{res}{numeric. resolution of canopy height model computed with \code{\link{points2DSM}} before tree segmentation}
\item{func}{a function to be applied to the attributes of extracted trees (return from internal call to treeExtraction function to compute plot level metrics}
\item{...}{other parameters to be passed to \code{\link{treeSegmentation}}}
}
\value{
a dataframe with tree metrics in columns corresponding to LAS objects of the list (lines)
}
\description{
This function performs tree detection in each of LAS objects contained in a list and then computes summary statistics based on compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest).
- calls the tree detection and extraction functions on each object
- remove the detected trees located outside of the regions of interest defined by their centers and radius
- computes summary statistics of remaining detected trees.
}
\examples{
data(laschablais3)
# extract three point clouds from LAS object
llas <- list()
llas[["A"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581680, 10)
llas[["B"]] <- lidR::lasclipCircle(laschablais3, 974390, 6581680, 10)
llas[["C"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581640, 10)
# normalize point clouds
llas <- lapply(llas, function(x) {lidR::lasnormalize(x, lidR::tin())})
# compute tree metrics
cloudTreeMetrics(llas,
cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
10, res=0.5)
# compute metrics with user-defined function
# number of detected trees between 20 and 30 meters and their mean height
user.func <- function(x, area.ha=NA)
{
dummy <- x$h[which(x$h>10 & x$h<20)]
data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
}
cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
}
\seealso{
\code{\link{treeSegmentation}}, \code{\link{treeExtraction}}, \code{\link{stdTreeMetrics}}
}
......@@ -28,8 +28,8 @@ raster::extent(r) <- c(0,40,0,40)
raster::res(r) <- 1
# xy positions
xy <- data.frame(c(10,20,31.25,15),
c(10,20,31.25,25))
xy <- data.frame(x=c(10,20,31.25,15),
y=c(10,20,31.25,25))
# compute mask
mask1 <- rasterXYMask(xy, c(5, 8, 5, 5), r)
mask2 <- rasterXYMask(xy, c(5, 8, 5, 5), r, binary=FALSE)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment