Raster outputs not buffered in 0.13.0?

[samuelgrubinger]

Hi JR,
Thanks for all the work you have put in to keep lasR up and running. It is an amazing tool!

The following code was producing a seamless DTM last week. After the update to 0.13.0, it now produces raster outputs that appear not to be buffered, or with gaps between tiles. The point clouds look perfect. Any idea what could be happening?

Thanks!

    pipeline = reader(filter = "Classification == 2") +
            # create a point cloud of only the lowest ground point in each 50 cm grid cell
            filter_with_grid(1, operator = "min", filter = keep_ground()) +

            # generate a triangulation of these points
            (tri = triangulate(max_edge = 120, filter = "Classification == 2",
                             #use_low = TRUE, 
                              use_attribute = "Z")) +
            
            #rasterize the triangulation at a 1-m resolution 
           (dtm =  rasterize(1, tri, ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM.tif"))) +
            #(dtm = load_raster(paste0(root_dir, "/output/raster/DTM.tif"))) +

            # apply a laplachian filter to replace spikes and pits with median values
            (spikefree1 = pit_fill(dtm,
                        lap_size = 5L,
                      thr_lap = 0.04,
                      thr_spk = -0.04,
                      med_size = 5L,
                      dil_radius = 1L,
                                   ofile = "")) +

            # perform this filter again at a coarser scale to capture larger spikes
            pit_fill(spikefree1,
                       lap_size = 8L,
                      thr_lap = 0.3,
                      thr_spk = -0.12,
                      med_size = 8L,
                      dil_radius = 1L,
                     ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM_spikefree.tif"))
               

       
    lasR::exec(pipeline, 
         on = paste0(root_dir, "/output/_1_classified_points"),
         ncores = 16,
        with = list(buffer = 120),
         progress = FALSE)
    

[Jean-Romain]

It happens that I changed thousands of lines of code between 0.12.0 and 0.13.0 to re-engineer lasR and make it ready for being much less specific to LASlib and the LAS format and ready to support more generic point cloud formats. As mentioned in the change-log I'm expecting to find a few bugs in a close future. So here you found one.

Can you try to narrow down the problem before I start investigating. Do you have the problem with only one pit-fill iteration? Do you have the problem without any pit-fill stages (is DTM.tiff ok)? Can you reproduce with:

f <- system.file("extdata", "bcts/", package = "lasR")

[samuelgrubinger]

It occurs with all steps that produce raster outputs. I first noticed it with rasterize().

[Jean-Romain]

I need a reproducible example. All the unit tests are passing and the following is working just fine.

library(lasR)
f <- system.file("extdata", "bcts/", package = "lasR")
pipeline = rasterize(1)
chm = exec(pipeline, on = f)
terra::plot(chm)

[samuelgrubinger]

I tried it on these two test tiles and got the same result.

    pipeline = reader(filter = "Classification == 2") +
            # create a point cloud of only the lowest ground point in each 50 cm grid cell
            filter_with_grid(1, operator = "min", filter = keep_ground()) +

            # generate a triangulation of these points
            (tri = triangulate(max_edge = 120, filter = "Classification == 2",
                             #use_low = TRUE, 
                              use_attribute = "Z")) +
            
            #rasterize the triangulation at a 1-m resolution 
           (dtm =  rasterize(1, tri, ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM.tif")))
               

230526_161628_216_ground.zip

[Jean-Romain]

If I apply your pipeline on my data it works but I confirm I can reproduce on your data.

library(lasR)
f <- "/home/jr/Téléchargements/230526_161628_216_ground/"
pipeline = reader(filter = "Classification == 2")
tri = triangulate(max_edge = 120,  use_attribute = "Z")
dtm =  rasterize(1, tri)
pipeline = pipeline+tri+dtm

ans = exec(pipeline, on = f, verbose = T)
terra::plot(ans)

[Jean-Romain]

This was a very very silly bug. I wrote min_y instead of min_x somewhere in the code. The tests missed the bug because my test files are y-aligned and the coordinates are in a lucky range. Thank you for reporting. This is fixed.

[Jean-Romain]

By the way in

pipeline = reader(filter = "Classification == 2") +
            # create a point cloud of only the lowest ground point in each 50 cm grid cell
            filter_with_grid(1, operator = "min", filter = keep_ground()) +

            # generate a triangulation of these points
            (tri = triangulate(max_edge = 120, filter = "Classification == 2",
                             #use_low = TRUE, 
                              use_attribute = "Z")) +
            
            #rasterize the triangulation at a 1-m resolution 
           (dtm =  rasterize(1, tri, ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM.tif")))

"Classification == 2" is redundant. You already read only ground points. There is no other points in memory.

pipeline = reader(filter = "Classification == 2") +
            # create a point cloud of only the lowest ground point in each 50 cm grid cell
            filter_with_grid(1, operator = "min") +

            # generate a triangulation of these points
            (tri = triangulate(max_edge = 120,
                             #use_low = TRUE, 
                              use_attribute = "Z")) +
            
            #rasterize the triangulation at a 1-m resolution 
           (dtm =  rasterize(1, tri, ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM.tif")))

Last but not least, I did not remember that I added filter_with_grid. The new version of sampling_pixel does the same job and more. So filter_with_grid will likely be deprecated in future versions.

[samuelgrubinger]

Thanks so much for fixing this so quickly.
I reinstalled lasR. The kernel now crashes if I try to run a classification.

pipeline = reader(filter = "-thin_with_voxel .02
                            -keep_z_above 20
                            -drop_abs_scan_angle_above 15
                            -keep_return 1 2 3 4 5 6 7
                            -drop_number_of_returns 15
                            -drop_number_of_returns 14
                            -drop_number_of_returns 13
                            -drop_number_of_returns 12
                            -drop_number_of_returns 11
                            -drop_number_of_returns 10
                            -drop_number_of_returns 9
                            -drop_number_of_returns 8
                            -drop_withheld"
                                         ) +


        # classify points as noise if they have fewer than 200 points in their surrounding 10x10 m voxels,
        # for big clouds of above- or belowground noise
        classify_with_ivf(res = 10, n = 200L, class = 18L)`

The reader stage runs but the classification causes a crash.

[Jean-Romain]

Is it reproducible with the same files?

[samuelgrubinger]

Yep. It crashes instantaneously if I run on those files.

[Jean-Romain]

This is a random failure at read time because of an invalid read. I'll investigate. See #105

[samuelgrubinger]

Thanks JR. Any workaround for now to avoid reading the extra bytes? I tried running on CTG = readLAScatalog(path_als, select = "xyzacrnw") but it did not run.

[Jean-Romain]

readLAScatalog is from lidR. And the error seems occur when opening the file and it is not from my code. So even lidR should be affected I guess.

[samuelgrubinger]

I ran the following lidR code to produce clean input files without extra bytes:

lidR.progress = FALSE

path_als = paste0(root_dir, "/input")

# path to the output of this stage
data_path = paste0(root_dir, "/input_clean")

# define a function to drop extra bytes
drop_extra_bytes <- function(chunk) # user defined function
{
  las <- readLAS(chunk, select = "xyzarnw")  
  if (is.empty(las)) return(NULL)  
  return(las)
}

if (file.exists(data_path)) {

    CTG_clean = readLAScatalog(data_path)
    
} else {

    # create the directory if it doesn't exist
    dir.create(data_path, recursive = TRUE)

CTG = readLAScatalog(path_als)
#opt_select(CTG) = "xyzarnw"
opt_chunk_size(CTG) = 600
opt_output_files(CTG) = paste0(root_dir, "/input_clean/{XCENTER}_{YCENTER}")

CTG_clean = suppressWarnings(catalog_apply(CTG, drop_extra_bytes))
    }

It still crashes instantaneously when I try to run anything in lasR on the clean inputs.

[Jean-Romain]

The discussion about this bug is here #105