clip_roi vs. filter keep_xy

[wiesehahn]

I have a large catalog just tried to read part of it via

# attempt1
las = clip_roi(catalog, bbox)
# attempt2
las = readLAS(catalog, filter = paste("-keep_xy", bbox$xmin, bbox$ymin, bbox$xmax, bbox$ymax))

I expected attempt 2 to be faster but instead I aborted it after 5 minutes while attempt 1 was ready within 7 seconds.

Why is that?

[Jean-Romain]

I think the behavior is legit and expected but you did not provide enough details to be 100% sure of my answer. Please show print(catalog), plot(catalog), plot(bbox, add = TRUE), las_check(catalog)

[wiesehahn]

> print(ctg)
class       : LAScatalog (v1.2 format 1)
extent      : 341000, 675000, 5682000, 5972000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 32N 
area        : 49468.59 km²
points      : 655.22 billion points
density     : 13.2 points/m²
density      : 9.4 pulses/m²
num. files  : 49588 

> las_check(ctg)

 Checking headers consistency
  - Checking file version consistency... ✓
  - Checking scale consistency... ✓
  - Checking offset consistency... ✓
  - Checking point type consistency... ✓
  - Checking VLR consistency... ✓
  - Checking CRS consistency... ✓
 Checking the headers
  - Checking scale factor validity... ✓
  - Checking Point Data Format ID validity... ✓
 Checking preprocessing already done 
  - Checking negative outliers...
    ⚠ 28609 file(s) with points below 0
  - Checking normalization... no
 Checking the geometry
  - Checking overlapping tiles... ✓
  - Checking point indexation... yes

I didn't plot it for now, I guess the the bbox could not be recognized at that scale anyway. Could still do it if needed.

[Jean-Romain]

It's ok. Give me 10 minutes to anwser

[Jean-Romain]

There is a lot to say here.

First of all you cannot expect the second query to be faster than the first one for the simple reason that internally, lidR already does what is the best and the most efficient. So in the best case you can only be as fast as lidR.

Ok, now the question is, why is it slow as hell?

clip() analyses the bounding box of the query and removes all the files that are not overlapped by the bounding box. In your case it corresponds to maybe one or two files only I guess. Let say two files. The clip() query is consequently performed in 2 files only. You performed the query in 50000 files!

But the collection is spatially indexed! It should be fast?!?

This is correct. But the spatial index aware way to make a rectangular query is -inside not -keep_xy. -keep_xy is a regular filter. It keeps the points inside a rectangle while reading the files. It means that your second query needs to read your 50000 files to extract the rectangle. In 49998 files it will retain 0 point but the files are read anyway. This is why it takes hours...

So if I'm using paste("-inside", bbox$xmin, bbox$ymin, bbox$xmax, bbox$ymax) it will work?

Probably better but not as fast as expected. readLAS() does not parse the filter argument at the R level. The string is sent to LASlib the library that read and write LAS and LAZ files. So lidR won't keep the 2 relevant files. Instead lidR will ask LASlib to read the 50000 files and perform the query in each. For 49998 files this query will likely be super fast and for the remaining two files that actually contain points of interest the query will be fast. But yet, you need to read the header of the 50000 files. This takes a fraction of second times 50000 which is no longer fast. You already know that because I guess readLAScatalog() took a while to load the 50000 files.

As a conclusion clip() performs a spatial index aware query on two files. You performed a non spatial index aware query on 50000 files.

How I'm supposed to know about -inside vs -keep_xy?

As far as I know you can't. This is not part of the documentation of LASlib the library that drives LAStools. I discovered many thing by reading the source code of LASlib and -inside is not part of the documented output of LASlib.

[wiesehahn]

Thanks for your super fast and super detailed answer!