How to get a circular subset of a las-dataset with specific area using lidR::clip

To get a circular subset of a las-dataset with specific area, I would like to use lidR::clip_circular(). To do so, I first calculate the central point of my las-dataset, then I define the radius I want to use, to get a subset of exactly 500m^2 from the centroid of my las-dataset. The operation works and does not throw any error, the result however is not correct, see base::print() at the end of my short code snippet.

I have tried to use lidR::clip_roi() as well, providing a polygon representing my region of interest, but got the same, incorrect result. Now I do not have any clue how to go on. I could imagine, that it is about the crs I am using (EPSG:25832) or because the area is circular and not rectangular...

las_ctpt <- sf::st_coordinates(sf::st_centroid(sf::st_as_sfc(sf::st_bbox(las_norm, crs = crs_epsg25832), crs = crs_epsg25832), crs = crs_epsg25832))
# get the centroid of the (normalized) las-dataset

buff_radius <- base::sqrt(500/pi)
# calculate the radius to get a circular subset with an area of exactly 500m^2 of the las-dataset

las_subset <- lidR::clip_circle(las = las_norm, radius = buff_radius, xcenter = las_ctpt[, 1], ycenter = las_ctpt[, 2])
# subset the (normalized) las-file

base::print(las_subset)
class        : LAS (v1.4 format 6)
memory       : 3.7 Mb 
extent       : 78????.?, 78????.?, 528????, 528???? (xmin, xmax, ymin, ymax)
coord. ref.  : ETRS89 / UTM zone 32N 
area         : 604 m²
points       : 38.9 thousand points
density      : 64.45 points/m²
density      : 47.87 pulses/m²

las_subset_rect <- ((base::sqrt(500/pi))*2)^2
# this should be a rectangular box around the subset I want to have with an area of: (radius*2)^2

base::print(las_subset_rect)
[1] 636.6198

(Please note, that I have blackened out the exact location of my study area because of privacy issues. However, I am using EPSG:25832 which is a metric CRS used in Germany.)

CodePudding user response：

I can reproduce your issue with example data

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = lidR::readLAS(LASfile)
las_ctpt <- sf::st_coordinates(sf::st_centroid(sf::st_as_sfc(sf::st_bbox(las))))
buff_radius <- base::sqrt(500/pi)
las_subset <- lidR::clip_circle(las = las, radius = buff_radius, xcenter = las_ctpt[, 1], ycenter = las_ctpt[, 2])
las_subset
#> class        : LAS (v1.2 format 1)
#> memory       : 51.5 Kb 
#> extent       : 684867.3, 684892.2, 5017878, 5017903 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NAD83 / UTM zone 17N 
#> area         : 576 m²
#> points       : 863  points
#> density      : 1.5 points/m²
#> density      : 0.96 pulses/m²
#> sf::st_area(sf::st_convex_hull(las_subset))

However it is only an artifact of area estimation inaccuracies. The mathematical area of a point cloud is 0 m². All other measurements are approximations with various strategies. lidR computes the occupancy grid using a resolution of 2 meters. This is a change made in version 4.0 requested by users to mimic LAStools. It is more accurate for large dataset unevenly sampled but overestimates the area for small datasets. If you compute it with a convex hull (what lidR did before version 4.0) you get

sf::st_area(sf::st_convex_hull(las_subset))
#> 484.1633 [m^2]

Which is closer to the expected value