LAScatalog function (catalog_map) using filter_poi(), "Cannot write a file with 0 points"

[SteinarVatne]

I'm trying to identify and remove points found by the segment_shapes() function (lines and planes) on a LAScatalog, but I get the error message above. opt_output_files is set, and I've included the "if (is.empty(las)) return(NULL), but my current reasoning is that files are stored/overwritten for every step of the calculation:

las<-segment_shapes(las, shp_plane(k = 25), "planar")
las<-segment_shapes(las, shp_line(th = 8, k = 15), "linear")
las<-filter_poi(las, planar==FALSE, linear==FALSE)

If the result of filter_poi is a chunk with 0 points, then the last line will cause the error? Is there some way around this, is my reasoning wrong, or is there some other way to do this? Is there some way to process the intermediate steps in memory and only write the final chunk to disk?

[Jean-Romain]

Hard to help without a minimal reproducible example but I guess you could add if (is.empty(las)) return(NULL) after filter_poi()

[SteinarVatne]

Too bad! ;-) If files were saved only on the return from the function, the code should work...which would honestly make more sense. But it doesn't work.

Could this be different with a catalog_apply function, which has a separate handling of LAS and LAScluster objects?

[Jean-Romain]

The files are save on the return of the function. catalog_map is nothing else than catalog_apply with stuff internally done for the user. Again, without a minimal reproducible example I can't help. Here is my attempt at making a reproducible example myself. It works as expected

library(lasR)
test = function(las)
{
  las = filter_poi(las, X > 500000)
  if (is.empty(las)) return(NULL)
  return(las)
}

f1 <- system.file("extdata", "Megaplot.laz", package="lidR")
f2 <- system.file("extdata", "MixedConifer.laz", package="lidR")
ctg = readLAScatalog(c(f1,f2))

opt_output_files(ctg) <- paste0(tempdir(), "/*")
u = catalog_map(ctg, test)
u
plot(u)

[SteinarVatne]

data.zip
I scraped together some test-data, attached.

To get inside the size requirements for attachments I had to cut some tiles, this results in some warnings about different files that are not part of the problem, just ignore those.

The issue is affected by tilesize and -placement, but should appear with this script:

library(lidR)
library(terra)
library(sf)
library(future)
plan(multisession)

myfilt<-function(las){
    las <- segment_shapes(las, shp_plane(k = 25), "planar")
    las <- segment_shapes(las, shp_line(th = 8, k = 15), "linear")
    las<-filter_poi(las, planar==FALSE, linear==FALSE)
    if (is.empty(las)) return(NULL)
    return(las)
}

tmpfolder<-"D:/data/tmp/"

las = readLAScatalog("D:/data", filter="-keep_class 1 2")

opt_chunk_buffer(las) <- 20
opt_chunk_size(las) <- 800
opt_laz_compression(las) <- TRUE

opt_output_files(las)<- paste0(tmpfolder,"{XLEFT}_{YBOTTOM}_edit")
las<-catalog_map(las,myfilt)

[Jean-Romain]

You found a bug. It is a limit case. In chunk 2 when exiting your myfilt function it remains a single point. This point is in the buffer and is removed automatically by catalog_map that handles the buffer for you. So it remains 0 points inside the tile outside your function. You can handle that using the good old catalog_apply function but I fixed it in 4.1.1 available on github. Install the latest version or use:

myfilt<-function(las)
{
  if (is(las, "LAS"))
  {
    las <- segment_shapes(las, shp_plane(k = 25), "planar")
    las <- segment_shapes(las, shp_line(th = 8, k = 15), "linear")
    las<-filter_poi(las, planar==FALSE, linear==FALSE)
    return(las)
  }
  else if (is(las, "LAScluster"))
  {
    las = readLAS(las)
    if (lidR::is.empty(las)) return(NULL)
    las = myfilt(las)
    if (lidR::is.empty(las)) return(NULL)
    las = filter_poi(las, buffer == 0)
    if (lidR::is.empty(las)) return(NULL)
    return(las)
  }
  else if (is(las, "LAScatalog"))
  {
    options <- list(need_buffer = TRUE, drop_null = TRUE, need_output_file = TRUE)
    output  <- catalog_sapply(las, myfilt, .options = options)
    return(output)
  }
}

[SteinarVatne]

Thank you, that works! For posterity, I got an error from the "..." part in catalog_sapply, perhaps something you forgot to remove?

[Jean-Romain]

Yes, remove it