#'@title
#' clip plots or polygons form lidar
#'
#'@description
#' clip plots or polygons form lidar - requires that one first make a project using lasR_project which
#' includes both lidar and dtms
#'
#'@details
#' There are three ways to provide the input shapes to this function - either as a
#' list of plot ids, coordinates and plot diameters (4 column data.frame), as a list of vertices
#' organized by plot_id (3 column data.frame), or as a polygon shapefile (path to shapefile) inwich case
#' the plot id field must be specified "id_field_plots".
#'
#'\cr
#'Revision History
#' \tabular{ll}{
#'1.0 \tab date and revisions.. \cr
#'}
#'
#'@author
#'
#'Jacob Strunk <Jstrunk@@fs.fed.us>
#'
#'@param idxyd dataframe with point id, x and y coordinates, and circle diameter
#'@param idxy (optional) dataframe with polygon vertices (x,y) organized by point id
#'@param plot_polys (optional) shapefile path, or sp object with plot boundaries
#'@param id_field_plots (optional) in the event that a polygon is provided
#'@param lasR_project_polys spatial polygons dataframe (intersection between dtm and las tiles)
#'@param lasR_project path to spatial polygons dataframe (intersection between dtm and las tiles)
#'@param plot_tile_intersect path to spatial polygons dataframe representing intersection between plots and lasR_project
#'@param dir_out where to send clipped plots
#'@param height T/F subtract the ground
#'@param do_plot T/F make a plot of the intersections
#'@param return T/F return intersectiosn
#'@param n_core number of cores to use in clipping
#'@param dir_dtm in case path to dtms has changed from lasR_project
#'@param dir_las in case path to las has changed from lasR_project
#'@param skip_existing skip csv files that have been processed already - in case processing was interupted
#'@return
#' <Delete and Replace>
#'
#'@examples
#' <Delete and Replace>
#'
#'@import rgdal plyr rgeos raster lidR data.table parallel
#'
#'@export
#
#'@seealso \code{\link{another_function}}\cr \code{\link{yet_another_function}}\cr
#Desired upgrades to this function:
#
# x = function(x){}
#copy function arguments and use this code to format arguments
##writeClipboard(paste(gsub("^[[:space:]]*[,]*","#'@param ",gsub("=.*"," ",readClipboard())),collapse="\n"))
clip_plots=function(
idxyd=NA #id,x,y,diameter
,idxy=NA #id,xy coordinates of vertices for each polygon
,plot_polys=NA #shapefile path, or sp object for plots
,id_field_plots="plot"
,lasR_project=NA
,lasR_project_polys=NA
,plot_tile_intersect=NA
,dir_out=NA
,height=F
,do_plot=F
,return=F
,n_core=6
,dir_dtm=NA #in case drive paths are wrong (External drives...)
,dir_las=NA #in case drive paths are wrong (External drives...)
,skip_existing=T
,fix_dsm_bug = F # laz files for hood canal DSMs are screwed up
){
if(interactive()){
require(rgdal)
require(plyr)
require(rgdal)
require(rgeos)
require(raster)
require(lidR)
require(data.table)
require(parallel)
}
if(!file.exists(dir_out)) try(dir.create(dir_out, recursive=T),silent=T)
dir_skip=file.path(dir_out,"skip")
if(!file.exists(dir_skip)) try(dir.create(dir_skip, recursive=T),silent=T)
if(is.na(plot_tile_intersect)){
#load lasR_project
if(!is.na(lasR_project) & is.na(lasR_project_polys[1])){
if(!inherits(lasR_project,"SpatialPolygonsDataFrame")){
proj=read.csv(lasR_project,stringsAsFactors =F)
proj_polys0=bbox2polys(proj[,c("tile_id","mnx","mxx","mny","mxy")])
row.names(proj)=proj[,"tile_id"]
proj_polys=SpatialPolygonsDataFrame(proj_polys0,proj)
}
if(inherits(lasR_project,"SpatialPolygonsDataFrame")) proj_polys=lasR_project
}
if(!is.na(lasR_project_polys[1])){
if(!inherits(lasR_project_polys,"SpatialPolygonsDataFrame")) proj_polys=readOGR(lasR_project_polys,stringsAsFactors=F)#readOGR(dirname(lasR_project_polys),basename(lasR_project_polys))
if(inherits(lasR_project_polys,"SpatialPolygonsDataFrame")) proj_polys=lasR_project_polys
}
print("load lasR_project");print(Sys.time())
#create sp objects for plots
plot_polys_in=NULL
if(!is.na(plot_polys[1])){
if(inherits(plot_polys,"SpatialPolygonsDataFrame")) plot_polys_in=plot_polys
if(inherits(plot_polys,"character")) plot_polys_in=readOGR(plot_polys,stringsAsFactors=F)#readOGR(dirname(plot_polys),gsub("[.]shp","",basename(plot_polys)))
if(is.null(plot_polys)) stop("Class of plot_polys must be character or SpatialPolygonsDataFrame")
}
if(!is.na(unlist(idxy)[1]) & !inherits(plot_polys_in,"SpatialPolygonsDataFrame")){
plot_polys_in=points2polys(idxy)
}
if(!is.na(unlist(idxyd)[1]) & !inherits(plot_polys_in,"SpatialPolygonsDataFrame")){
seq1=c(seq(-pi,pi,.1),pi+.1)
circle1=data.frame(sin(seq1),cos(seq1))
spl_idxyd=split(idxyd,1:nrow(idxyd))
idxy=rbind.fill(lapply(spl_idxyd,function(x,circle1)data.frame(id=x[,1],x=circle1[,1]*x[,4]+x[,2],y=circle1[,2]*x[,4]+x[,3]),circle1))
row.names(idxyd)=idxyd[,1]
plot_polys_in0=points2polys(idxy)
plot_polys_in=SpatialPolygonsDataFrame(plot_polys_in0,data=idxyd)
}
#fix row names
row.names(plot_polys_in)=as.character(plot_polys_in@data[,id_field_plots])
print("Get / create Plot Polys");print(Sys.time())
#clean up self intersections
proj_polys_b=gBuffer(proj_polys, byid=TRUE, width=0)
plot_polys_b=gBuffer(plot_polys_in, byid=TRUE, width=0)
#clip data to match extents
ext_tile=as(extent(as.vector(t(bbox(proj_polys_b)))), "SpatialPolygons")
ext_plot=as(extent(as.vector(t(bbox(plot_polys_b)))), "SpatialPolygons")
plot_polys_ext=gIntersection(plot_polys_b,ext_tile, byid=T,drop_lower_td=T)
proj_polys_ext=gIntersection(proj_polys_b,ext_plot, byid=T,drop_lower_td=T)
#patch data back onto clipped polygons
row.names(proj_polys_ext)=gsub(" 1","",row.names(proj_polys_ext))
row.names(plot_polys_ext)=gsub(" 1","",row.names(plot_polys_ext))
keep_proj=row.names(proj_polys_b@data) %in% row.names(proj_polys_ext)
keep_plot=row.names(plot_polys_b@data) %in% row.names(plot_polys_ext)
proj_polys_b1=proj_polys_b[keep_proj,]
plot_polys_b1=plot_polys_b[keep_plot,]
plot_polys_spdf=SpatialPolygonsDataFrame(plot_polys_ext,plot_polys_b1@data)
proj_polys_spdf=SpatialPolygonsDataFrame(proj_polys_ext,proj_polys_b1@data)
print("clip, recombine plots with data");print(Sys.time())
#intersect plots with tiles
test = gIntersection(plot_polys_spdf,proj_polys_spdf, byid=T)
proj_plot_x=gIntersects(plot_polys_spdf,proj_polys_spdf, byid=T,returnDense=F)
proj_plot_x1=proj_plot_x[sapply(proj_plot_x,function(x)length(x)>0)]
print("intersect plots and tiles");print(Sys.time())
#parse intersections and compile data
fn_parse=function(x,y,plots,tiles){
data.frame(plots[as.character(x),,drop=F],tiles[y,],row.names=NULL)
}
plots_tiles=rbind.fill(mapply(fn_parse
,x=names(proj_plot_x1)
,y=proj_plot_x1
,SIMPLIFY = F
, MoreArgs =
list(
plots=plot_polys_spdf@data
,tiles=proj_polys_spdf@data
)
)
)
#merge duplicate records
dup_id=.dup2(as.character(plots_tiles[,id_field_plots]))
dups=plots_tiles[dup_id,]
no_dups_df=plots_tiles[!dup_id,]
spl_dups=split(dups,as.character(dups[,id_field_plots]))
dups_df=.fn_merge(spl_dups)
plots_tiles_unq=rbind(no_dups_df,dups_df)
plots_tiles_unq[,id_field_plots]=as.character(plots_tiles_unq[,id_field_plots])
row.names(plots_tiles_unq)=plots_tiles_unq[,id_field_plots]
print("merge duplicates");print(Sys.time())
#add records to geometry
good_polys=names(plot_polys_ext) %in% as.character(plots_tiles_unq[,id_field_plots])
plot_polys_merge=SpatialPolygonsDataFrame(plot_polys_ext[good_polys,],plots_tiles_unq)
}else{
plot_polys_merge=readOGR(plot_tile_intersect,stringsAsFactors=F)
}
print("fix paths");print(Sys.time())
#fix drive paths in lasR_project
if(!is.na(dir_dtm)) plot_polys_merge@data[,"dtm_file"]=unlist(lapply(plot_polys_merge@data[,"dtm_file"],function(...,dir_dtm)paste(file.path(dir_dtm,basename(strsplit(...,",")[[1]])),collapse=","),dir_dtm=dir_dtm))
if(!is.na(dir_las)) plot_polys_merge@data[,"las_file"]=unlist(lapply(plot_polys_merge@data[,"las_file"],function(...,dir_las)paste(file.path(dir_las,basename(strsplit(...,",")[[1]])),collapse=","),dir_las=dir_las))
if(do_plot ){
if(! is.na(plot_tile_intersect)) plot(proj_polys_spdf)
plot(plot_polys_merge,border="red",add=T,lwd=10)
}
print("skip existing");print(Sys.time())
#skip existing files
if(skip_existing){
files_out_dir=unlist(c(list.files(dir_out,pattern="[.]las"),list.files(dir_skip,pattern="[.]las")))
out_nms=paste("plot_",plot_polys_merge@data[,id_field_plots],".las",sep="")
plot_polys_merge=plot_polys_merge[!out_nms %in% files_out_dir,]
}
print("write shapefile");print(Sys.time())
#write shapefile of intersections
dir_overlap=file.path(dir_out,"plot_tile_overlap")
if(!dir.exists(dir_overlap)) dir.create(dir_overlap)
if(!file.exists(paste(dir_overlap,"plot_merge_tiles.shp",sep="\\"))) writeOGR(plot_polys_merge,dir_overlap,"plot_merge_tiles", driver="ESRI Shapefile")
print("split polygons");print(Sys.time())
#clip points
spl_plots=sp::split(plot_polys_merge,1:nrow(plot_polys_merge))
print("initiate clipping");print(Sys.time())
if(n_core>1){
clus=makeCluster(n_core)
clusterEvalQ(clus,library(lasR))
res=parLapply(clus,spl_plots,.try_clip_plots,dir_out = dir_out,height=height,id=id_field_plots,fix_dsm_bug=fix_dsm_bug)
stopCluster(clus)
}
if(n_core<2){
lapply(spl_plots,.try_clip_plots,dir_out = dir_out,height=height,id=id_field_plots,fix_dsm_bug=fix_dsm_bug)
}
#.try_clip_plots(spl_plots[[2]],dir_out = dir_out)
print("clip plots");print(Sys.time())
print("write outputs");print(Sys.time())
if(return) return(plot_polys_merge)
}
.dup2=function(x,type=c("all","first","last"),value=F,invert=F){
if(type[1]=="all") index = duplicated(x) | duplicated(x, fromLast=TRUE)
if(type[1]=="first") index = duplicated(x)
if(type[1]=="last") index = duplicated(x, fromLast=TRUE)
if(invert) index = !index
if(value) return(x[index])
if(!value) return(index)
}
.fn_merge=function(x){
if(interactive()){ require(data.table)}
if(class(x)=="list") x_in=data.frame(rbindlist(lapply(x,.fn_merge)))
else{
x_in=x[1,]
x_in[,"tile_id"]=paste(as.character(x[,"tile_id"]),collapse=",")
x_in[,"las_file"]=paste(unique(strsplit(paste(x[,"las_file"],collapse=","),",")[[1]]),collapse=",")
x_in[,"dtm_file"]=paste(unique(strsplit(paste(x[,"dtm_file"],collapse=","),",")[[1]]),collapse=",")
if("x.1" %in% names(x_in)) x_in[,"x.1"]=paste(x[,"x.1"],collapse=",")
if("y.1" %in% names(x_in)) x_in[,"y.1"]=paste(x[,"y.1"],collapse=",")
x_in[,"mnx"]=min(x[,"mnx"])
x_in[,"mxx"]=max(x[,"mxx"])
x_in[,"mny"]=min(x[,"mny"])
x_in[,"mxy"]=max(x[,"mxy"])
}
return(x_in)
}
.try_clip_plots=function(...){
#require(lasR)
.clip_plots=function(x,id,dir_out,return=F,height=T, fix_dsm_bug = F){
if(interactive()){
require(lidR)
#require(lasR)
require(plyr)
}
print()
las_files_in = grep("[.]la.$",as.character(unlist(strsplit(x@data[,"las_file"],",")[1])),value=T)
dtm_files_in = grep("[.].{,4}$",unlist(strsplit(x@data[,"dtm_file"],",")[1]),value=T)
las_in=readLAS(files = las_files_in)
if(fix_dsm_bug) las_in@header@PHB['Header Size'] = 235
dtm_in=read_dtm(dtm_files_in)
dtm_poly=try(crop(dtm_in,gBuffer(x,width=20)))
if(class(dtm_poly)=="try-error"){warning("plot and dem do not intersect, plot: ",x@data[,1]);return()}
if(length(x@polygons[[1]]@Polygons)>1){
#clip to plot bbox
x_ext=as(extent(x),"SpatialPolygons")
las_clip_in=lasclip(las_in, x_ext@polygons[[1]]@Polygons[[1]], inside = TRUE)
#clip individual sub polygons
las_poly_list=lapply(x@polygons[[1]]@Polygons,function(x,las){lasclip(las, x@polygons[[1]]@Polygons[[1]], inside = TRUE)}, las_clip_in)
#merge las from subplots back together
points_merge=do.call(rbind,lapply(las_poly_list,function(x)x@data))
las_poly=LAS(points_merge,las_poly_list[[1]]@header)
}else{
#simple clip for contiguous polygon
#poly_coords=x@polygons[[1]]@Polygons[[1]]@coords
#las_poly=lasclip(las_in, "polygon", poly_coords , inside = TRUE)
las_poly=lasclip(las_in, x@polygons[[1]]@Polygons[[1]] , inside = TRUE)
}
#write empty file if there are no points in las_poly
is_skip = class(las_poly)!="LAS"
if(!is_skip) is_skip = length(las_poly$X) == 0
if(is_skip){
skip_file=file.path(dir_out,"skip",paste(id,"_",x@data[1,id],".las",sep=""))
file.create(skip_file)
return()
}
if(height) las_hts = lasnormalize(las_poly, dtm_poly)
if(!height) las_hts = las_poly
#write to file
if(!dir.exists(dir_out)) dir.create(dir_out)
if(class(las_hts)=="LAS"){
out_file_i=file.path(dir_out,paste(id,"_",x@data[1,id],".laz",sep=""))
err=try(writeLAS(las_hts,out_file_i))
}else{
skip_file=file.path(dir_out,"skip",paste(id,"_",x@data[1,id],".las",sep=""))
file.create(skip_file)
}
}
try(.clip_plots(...))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.