r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
614 stars 134 forks source link

undocumented clip_roi behaviour with overlapping polygons #728

Closed mbedward closed 12 months ago

mbedward commented 1 year ago

When working with an 'sf' polygon layer where some features overlap, I hadn't realized that the clip_roi function only allocates a point to one of the features in an overlapping group. This isn't mentioned in the function docs and was the cause of unexpected values in downstream calculations (veg cover estimates) before I realized it was happening.

A simple example to reproduce the behaviour...

library(sf)
library(lidR)

# Load a LAS tile
# las <- lidR::readLAS(path)

# Create two overlapping polygons inside the tile bounds
bb <- st_bbox(las)
y <- mean(bb[c("ymin", "ymax")])
w <- bb["xmax"] - bb["xmin"]
x1 <- bb["xmin"] + 0.4*w
x2 <- bb["xmin"] + 0.6*w

pts <- st_sfc(st_point(c(x1, y)), st_point(c(x2, y)), crs = st_crs(las))
polys <- st_buffer(pts, dist = 0.2*w)

poly_las <- clip_roi(las, polys)

# Base R plot of a sample of the points
xyplot <- function(xlas, prop_points = 0.05) {
  xy <- cbind(xlas$X, xlas$Y)
  n <- nrow(xy)
  xy <- xy[sample(n, size=prop_points * n), ]
  plot(xy, asp=c(1,1), pch=16, cex=0.2)
}

xyplot(poly_las[[1]])  # full circle of points for first polygon
xyplot(poly_las[[2]])  # cresent shape (non-overlapping portion) for second polygon

A simple work-around is to run clip_roi via lapply...

poly_las <- lapply(seq_along(polys), function(i) clip_roi(las, polys[[i]] ))

Perhaps it would be good to mention this in the function docs?

Jean-Romain commented 1 year ago

This is not documented because it is completely unexpected. Thank you for the reproducible example.

Jean-Romain commented 12 months ago

Fixed

mbedward commented 12 months ago

Thanks - it's working perfectly for me now.