How to homogenize points using lidR to reduce streaking artifacts

96 views Asked by At

I am trying to homogenize an aerial LAS catalog to achieve a consistent point density across a large study area. The dataset ranges from less than 4 points per meter to greater than 8 points per meter, with many areas of extreme point density in overlapping swathes. Certain raster metrics I calculate generate streaking artifacts with lines of low value pixels occurring where high point density meets lower point density (see screenshot of output 10m IQR). I hope to fix this by thinning points to a more even consistency by using the homogenize option in lidR::decimate_points(). While running this function, a came across another issue where summary() reports a successful decimation, but the output files appear the same as the input files.

I tried numerous density and resolution combinations to ensure I wasn't making simple USft to meter conversion errors. The R book mentions that attempting to achieve higher point densities than original will result in an identical file being created, but I have attempted some fairly modest densities/resolutions without luck.

To ensure that my input tiles weren't overlapping and creating duplicate points, I also tried adjusting the opt_chunk_size to both higher and lower values.

In summary:

  1. Can homogenize remove streaking at the edges of overlapping points?

Update: Homogenizing points didn't help much and I suspect errors in calibration/post-processing

  1. Why does R indicate successful homogenization, while the output files remain the same?

Update: I solved this problem by simply restarting PC and rerunning code with a different working directory.

10m interquartile range [IQR] with streaking at edges with point cloud overlaid and colored by point density visualized in Global Mapper

In R:

library(lidR) #lidR version 4.0.4
library(future)
#load normalized LAZ tiles
ctg <- readLAScatalog("path/", filter = "-drop_class 2") #drop ground classes from normalized tiles
opt_chunk_buffer(ctg)<- 200
opt_chunk_size(ctg)<- 0 #maintain chunk size 750m = 2460.63ft

summary(ctg) #input catalog
class       : LAScatalog (v1.1 and 1.2 format 1)
extent      : x1, x2, y1, y2 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83 / New York Central (ftUS) 
area        : 2275.8 kus-ft²
points      : 10.74 billion points
density     : 4.7 points/us-ft²
density     : 2.5 pulses/us-ft²
num. files  : 422 
proc. opt.  : buffer: 200 | chunk: 0
input opt.  : select: * | filter: -drop_class 2
output opt. : on disk | w2w guaranteed | merging enabled

plan(multisession, workers = 3L) 
opt_output_files(ctg)<- paste0("path/retile\_{XLEFT}\_{YBOTTOM}") # select output and naming for tiles 

ctg_thinned = decimate_points(ctg, homogenize(density = 1, res= 10.76, use_pulse = FALSE)) #1 pt/~1m^2 in feet

summary(ctg_thinned) #output catalog
class       : LAScatalog (v1.1 and 1.2 format 1)
extent      : x1, x2, y1, y2 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83 / New York Central (ftUS) 
area        : 2257.83 kus-ft²
points      : 1.69 billion points
density     : 0.7 points/us-ft²
density     : 0.5 pulses/us-ft²
num. files  : 421 
proc. opt.  : buffer: 5 | chunk: 0
input opt.  : select: * | filter: -drop_class 2
0

There are 0 answers