grid_canopy(p2r()) output size depends on subcircle tweak / breaks with res=raster

[Hannes-Ole]

When using grid_canopy(), the output is dependend on the subcircle tweak used in p2r(). This makes sense, since after "expanding" each point to a disk, the extent of the point cloud is larger - it can be argued that it makes sense, that the user gets as chm the full extent of the "extended" point cloud - although having a larger output extent than input extent may still be surprising to some.

However the algorithm breaks when using a reference raster to determine alignment/resolution, and I think that shouldn't happen (I'd rather get the "cropped" chm). Here is a code example, if it isn't reproducible on your data I can send you some:

# create reference raster
ref_raster <- raster(extent(tile_dtm), res = 1)
crs(ref_raster) = crs(tile_dtm)

# make sure the las cloud extent is the same as the reference raster
tile_las = clip_roi(tile_las, extent(ref_raster))

Calling the following causes the error:

>   r_chm_bottom = grid_canopy(tile_las, res=ref_raster, p2r(0.5))
Fehler in C_rasterize(las, layout, subcircle, 1L) : 
  C++ unexpected internal error in 'rasterize': point out of raster.

Here we see, what the algorithm would like to output, without being restricted by an extent:

>   r_chm = grid_canopy(tile_las, res=1, p2r(0.5))
>   extent(r_chm_bottom)
class      : Extent 
xmin       : 2692999 
xmax       : 2694001 
ymin       : 1124999 
ymax       : 1126001 

When using grid_metrics or setting p2r(subcircle = 0) , which is the default, it works without complaints. A side question: Is grid_metrics(~max(Z)) basically equivalent to grid_canopy(p2r())?

>   r_npoints_all = grid_metrics(tile_las, ~max(Z), res=ref_raster)
>   extent(r_npoints_all)
class      : Extent 
xmin       : 2693000 
xmax       : 2694000 
ymin       : 1125000 
ymax       : 1126000 

Also the pitfree algorithm apparently has a crop back to the original extent implemented, as grid_canopy(pitfree(subcircle > 0)) also works fine:

>   r_chm = grid_canopy(tile_las, res=ref_raster, pitfree(thresholds = c(0, 2, 5, 10, 15), max_edge = c(0, 1), subcircle = 0.5) )
>   extent(r_chm_bottom)
class      : Extent 
xmin       : 2693000 
xmax       : 2694000 
ymin       : 1125000 
ymax       : 1126000 

I suggest having a p2r optional parameter cropToInputExtent = TRUE (you may have a more concise name for it ;] ), that allows for extension by the subcircle tweak, but defaults to extent(input) = extent(output) [as generally may be expected and allows for usage of a reference raster grid_canopy(las, res=raster(extent(las), res=1), p2r(0.5)) ].

[Jean-Romain]

C++ unexpected internal error in 'rasterize': point out of raster.

This problem has been fixed in https://github.com/r-lidar/lidR/issues/483 and released in 3.2.2. Which version are you using?

t can be argued that it makes sense, that the user gets as chm the full extent of the "extended" point cloud - although having a larger output extent than input extent may still be surprising to some.

I do agree and this is something I was thinking about for a while. I removed this in v4. The extent of the point cloud is now fixed and is not affected by the subcircle tweak

grid_metrics(~max(Z)) basically equivalent to grid_canopy(p2r())

Yes, but ten times slower approximately.

I suggest having a p2r optional parameter cropToInputExtent = TRUE

I think that do not extending the raster is simpler, more consistent and more logic.

[Hannes-Ole]

You are right, I was working on a server where still lidR 3.2.1 was installed, so this is obsolete then. I also agree that not extending the raster is consistent and the optional extension isn't really necessary. Thank you!