R: How do I loop through spatial points with a specific buffer? -
so problem quite difficult describe hope can make question clear possible.
use rlidar package
load .las
file r , afterwards convert spatialpointsdataframe using sp package
. spatialpointsdataframe quite dense.
now want define buffer of 0.5 meters , loop (iterate) him (the buffer) through points, choosing point highest z value within buffer, next point jump to.this should repeated until there isn't point within buffer higher z value current. all values (or perhaps x , y values) of "found" point should written list/dataframe , process should repeated until such highest points found. thats code got far:
>library(rlidar) >library(sp) >rlas<-readlas("test.las",short=false) >pointcloud<- data.frame(rlas) >coordinates(pointcloud) <- c("x", "y")
well googled extensively not find clues how proceed further...
dont know packages of help, guess perhaps spatstat
question go spatial point pattern analysis.
have ideas how archive in r? or not possible? (do perhaps have skip python make work?)
gladly appreciated.
if want set of points local maxima within 0.5m radius circle around each point, should work. gist of is:
- convert las points spatialpointsdataframe
- create buffered polygon set overlapping polygons
- loop through buffered polygons , find desired element within buffer -- in case, it's 1 maximum height.
code below:
library(rlidar) library(sp) library(rgeos) rlas <- readlas("test.las",short=false) pointcloud <- data.frame(rlas) coordinates(pointcloud) <- c("x", "y")
finish creating spatialpointsdataframe las source. i'm assuming field point height pointcloud$value
pointcloudspdf <- spatialpointsdataframe(data=pointcloud,xy)
use rgeos library intersection. it's important have byid=true
or polygons merged intersect
bufferedpoints <- gbuffer(pointcloudspdf,width=0.5,byid=true) # save our local maxima state (this updated) localmaxes <- rep(false,nrow(pointcloud)) i=0 (buff in 1:nrow(bufferedpoint@data)){ <- i+1 bufpolygons <- bufferedpoints@polygons[[i]] bufsppolygons <- spatialpolygons(list(bufpolygons)) bufsppolygondf <-patialpolygonsdataframe(bufsppolygons,bufferedpoints@data[i,]) ptsinbuffer <- which(!is.na(over(pointcloudspdf,sppolygondf))) # i'm assuming `value` field name containing point height localmax <- order(pointcloudspdf@data$value[ptsinbuffer],decreasing=true)[1] localmaxes[localmax] <- true } localmaxpointclouddf <- pointcloudspdf@data[localmaxes,]
now localmaxpointclouddf
should contain data original points if local maximum. warning -- isn't going super fast if have lot of points. if ends being concern may smarter pre-filtering points using smaller grid , extract
raster
package.
that this:
make cell size small enough each 0.5m buffer intersect @ least 4 raster cells -- err on smaller since comparing circles squares.
library(raster) numrows <- extent(pointcloudspdf)@ymax-extent(pointcloudspdf)@ymin/0.2 numcols <- extent(pointcloudspdf)@xmax-extent(pointcloudspdf)@xmin/0.2 emptyraster <- raster(nrow=numrows,ncol=numcols)
rasterize
create grid maximum value of given field within cell. because of square/circle mismatch starting point filter out obvious non-maxima. after have raster in local maxima represented cells. however, won't know cells maxima in 0.5m radius , don't know point in original feature layer came from.
r <- rasterize(pointcloudspdf,emptyraster,"value",fun="max")
extract
give raster values (i.e., highest value each cell) each point intersects. recall above local maxima in set, although values not 0.5m radius local maxima.
rastermaxes <- extract(r,pointcloudspdf)
to match original points raster maxes, subtract raster value @ each point point's value. if value 0, values same , have point potential maximum. note @ point merging points raster -- have throw of these out because "under" 0.5m radius higher local max though max in 0.2m x 0.2m cell.
potentialmaxima <- which(pointcloudspdf@data$value-rastermaxes==0)
next, subset original spatialpointsdataframe , we'll more exhaustive , accurate iteration on subset of points since should have thrown out bunch of points not have been maxima.
potentialmaximacoords <- coordinates(pointcloudspdf@coords[potentialmaxima,]) # using data.frame() constructor because example has 1 column potentialmaximadf <- data.frame(pointcloudspdf@data[potentialmaxima,]) potentialmaximaspdf <-spatialpointsdataframe(potentialmaximacoords,potentialmaximadf)
the rest of algorithm same buffering smaller dataset , iterating on it:
bufferedpoints <- gbuffer(potentialmaximaspdf, width=0.5, byid=true) # save our local maxima state (this updated) localmaxes <- rep(false, nrow(pointcloud)) i=0 (buff in 1:nrow(bufferedpoint@data)){ <- i+1 bufpolygons <- bufferedpoints@polygons[[i]] bufsppolygons <- spatialpolygons(list(bufpolygons)) bufsppolygondf <-patialpolygonsdataframe(bufsppolygons,bufferedpoints@data[i,]) ptsinbuffer <- which(!is.na(over(pointcloudspdf, sppolygondf))) localmax <- order(pointcloudspdf@data$value[ptsinbuffer], decreasing=true)[1] localmaxes[localmax] <- true } localmaxpointclouddf <- pointcloudspdf@data[localmaxes,]
Comments
Post a Comment