Nothing
#'Individual trees crown deliniation from LiDAR-derived Canopy Height Model (CHM)
#'
#'@description Delineate and compute ground-projected area of individual tree crowns detected from LiDAR-derived CHM
#'
#'@usage ForestCAS(chm, loc, maxcrown, exclusion)
#'
#'@param chm A LiDAR-derived Canopy Height Model (CHM) RasterLayer or SpatialGridDataFrame file.
#'@param loc A matrix or dataframe with three columns (tree xy coordinates and height).
#'@param maxcrown A single value of the maximum individual tree crown radius expected. Default 10.0 m.
#'@param exclusion A single value from 0 to 1 that represents the \out{\%} of pixel exclusion. E.g. a value of 0.5 will exclude all of the pixels for a single tree that has a height value of less than 50\out{\%} of the maximum height from the same tree. Default value is 0.3.
#'@return Returns a list that contains the individual tree canopy boundary polygons and the 4-column matrix with the tree xy coordinates, heights and ground-projected canopy area (with units of square meters).
#'@author Carlos Alberto Silva
#'@examples
#' \donttest{
#' # This takes > 5s
#' library(raster)
#'
#' # Import the LiDAR-derived CHM file
#' data(chm) # or set a CHM. e.g. chm<-raster("CHM_stand.asc")
#'
#' # Set the loc parameter
#' sCHM<-CHMsmoothing(chm, filter="mean", ws=5) # smoothing CHM
#' loc<-FindTreesCHM(sCHM, fws=5, minht=8) # or import a tree list
#'
#' # Set the maxcrown parameter
#' maxcrown=10.0
#'
#' # Set the exclusion parameter
#' exclusion=0.3 # 30
#'
#' # Compute individual tree detection canopy area
#' canopy<-ForestCAS(chm, loc, maxcrown, exclusion)
#'
#' #==================================================================================#
#' # Retrieving the boundary for individual tree detection and canopy area calculation
#' #==================================================================================#
#' boundaryTrees<-canopy[[1]]
#
#' # Plotting the individual tree canopy boundary over the CHM
#' plot(chm, main="LiDAR-derived CHM")
#'
#' # adding tree canopy boundary
#' plot(boundaryTrees, add = TRUE, border = 'red', bg = 'transparent')
#'
#' #============================================================================#
#' # Retrieving the list of individual trees detected for canopy area calculation
#' #============================================================================#
#' canopyList<-canopy[[2]] # list of ground-projected areas of individual tree canopies
#' summary(canopyList) # summary
#'
#' # Spatial location of the trees
#' library(sp)
#' XY<-SpatialPoints(canopyList[,1:2]) # Spatial points
#' plot(XY, col = "black", add = TRUE, pch = "*") # adding tree location to the plot
#' }
#'@importFrom spatstat.geom disc
#'@importFrom sp coordinates gridded Polygon Polygons SpatialPolygons over
#'@importFrom sp over SpatialPoints SpatialPolygonsDataFrame SpatialGridDataFrame SpatialPixelsDataFrame
#'@importFrom deldir tile.list deldir
#'@importFrom plyr ddply
#'@importFrom raster raster rasterToPolygons boundaries
#'@export
ForestCAS<-function(chm,loc,maxcrown=10,exclusion=0.3) {
if (class(chm)!="RasterLayer" & class(chm)!="SpatialGridDataFrame") {stop("The chm is invalid. It must to be a RasterLayer or SpatialGridDataFrame'")}
if (ncol(loc)!=3) {stop("The input loc is invalid. It must to be 3-column matrix or dataframe with the x,y coordinates and heights of the individual trees")}
if (class(maxcrown)!="numeric") {stop("The maxcrown parameter is invalid. It is not a numeric input")}
if (class(exclusion)!="numeric") {stop("The exclusion parameter is invalid. It is not a numeric input")}
if (exclusion >=1) {stop("The exclusion parameter is invalid. It must to be less than 1numeric input")}
if (class(chm)=="RasterLayer"){ chm<-methods::as(chm, "SpatialGridDataFrame")}
Hthreshold<-min(loc[,3])*exclusion
polys<-list()
width<-numeric()
DF2raster<-function(h.mH, i){
h.mHkl<-subset(h.mH,h.mH[,4]==levels(factor(h.mH[,4]))[i])
if (nrow(h.mHkl==1)==TRUE) {
h.mHkl<-rbind(h.mHkl,c(h.mHkl[,1],h.mHkl[,2]+0.005,h.mHkl[,3]+0.005,h.mHkl[,4]))}
spP <- cbind(h.mHkl[,2:3],h.mHkl[,1],h.mHkl[,4])
colnames(spP)<-c("x","y","z","id")
#coordinates(spP)<- c("x", "y")
#suppressWarnings(gridded(spP) <- TRUE)
m <- suppressWarnings(SpatialPixelsDataFrame(points=spP[c("x", "y")], data = spP))
rasterDF <- raster(m)
hhg<- boundaries(rasterDF, type='outer')
p <- rasterToPolygons(hhg, dissolve=TRUE)
sp.polys <- p[1,]
return(sp.polys)
}
for(i in 1:nrow(loc)) {
width[i]=maxcrown
discbuff<-disc(radius=width[i], centre=c(loc$x[i], loc$y[i]))
discpoly<-Polygon(rbind(cbind(discbuff$bdry[[1]]$x,
y=discbuff$bdry[[1]]$y), c(discbuff$bdry[[1]]$x[1],
y=discbuff$bdry[[1]]$y[1])))
polys<-c(polys, discpoly)
}
spolys<-list()
for(i in 1:length(polys)) {
spolybuff<-Polygons(list(polys[[i]]), ID=row.names(loc)[i])
spolys<-c(spolys, spolybuff)
}
polybuffs<-SpatialPolygons(spolys)
chmdf<-as.data.frame(chm,xy=TRUE)
Points.Ply<-over(SpatialPoints(chmdf[,2:3]),polybuffs)
Points.PlyD<-cbind(chmdf,Points.Ply)
Points.PlyD<-na.omit(Points.PlyD)
vor = deldir(loc[,1], loc[,2], z=loc[,3],suppressMsge=T)
tile = tile.list(vor)
polys = vector(mode='list', length=length(tile))
for (i in seq(along=polys)) {
pcrds = cbind(tile[[i]]$x, tile[[i]]$y)
pcrds = rbind(pcrds, pcrds[1,])
polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
veronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=loc[,1],
y=loc[,2], row.names=sapply(methods::slot(SP, 'polygons'),
function(x) methods::slot(x, 'ID'))))
chmV<-over(SpatialPoints(Points.PlyD[,2:3]), SP)
RpD<-cbind(Points.PlyD,chmV)
RpD.filter<-subset(RpD[,1:5],RpD[,1]>=Hthreshold)
RpD.filter<-cbind(RpD.filter[,1:3],RpD.filter[,5])
colnames(RpD.filter)<-c("z","x","y","g")
h.mH<-ddply(RpD.filter,"g", function (RpD.filter)
subset(RpD.filter,RpD.filter[,1]>= max(RpD.filter[,1])*exclusion))
for ( j in 1:nlevels(factor(h.mH[,4]))){
assign(paste0("SP.polys", j), DF2raster(h.mH,j))
cat (".");utils::flush.console()}
polygons <- methods::slot(get("SP.polys1"), "polygons")
for (i in 1:nlevels(factor(h.mH[,4]))) {
data.loc <- get(paste0("SP.polys",i))
polygons <- c(methods::slot(data.loc, "polygons"),polygons)
}
for (i in 1:length(polygons)) {
methods::slot(polygons[[i]], "ID") <- paste(i)
}
spatialPolygons <- SpatialPolygons(polygons)
spdf <- SpatialPolygonsDataFrame(spatialPolygons,
data.frame(Trees=1:length(polygons)))
options(scipen=10)
spdf<-spdf[spdf@data[-length(polygons),],]
areaList<-sapply(methods::slot(spdf, "polygons"), methods::slot, "area")
canopyTable<-cbind(loc,areaList)
colnames(canopyTable)<-c("x","y","z","ca")
result=list(spdf,canopyTable)
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.