Large rasters cause integer overflow in raster_cell_from_xy()

[bfrank5]

Hello,

I noticed an issue attempting to use a large DTM (~100,000 ha, 0.3 m resolution) during the normalization of points. The error I receive is just a warning indicating that the normalization algorithm switches during the operation to some default. Regardless of this specific context, I narrowed down the issue to a problem in raster_cell_from_xy() which must convert (in my case) extremely large cell indices into integers, which exceeds the capacity of as.integer, forcing the indices to return as NA values. Here is a reproducible example:

library(terra)

extent <- ext(c(0, 1e05, 0, 1e05))
raster <- rast(extent, resolution = 1, vals = 0)

# Create an arbitrary point
x <- 50
y <- 50

# Get the cell index
lidR:::raster_cell_from_xy(raster, x, y)

Producing the output:

[1] NA
Warning message:
In lidR:::raster_cell_from_xy(raster, x, y) :
  NAs introduced by coercion to integer range

I believe the particular line in question is located here:

lidR/R/utils_raster.R

Line 75 in 09a878b

cell <- as.integer(row * ncol + col + 1)

[Jean-Romain]

Can you clip your dtm to the extent of your point cloud? At some point R has on 32 bits signed integer and 64 bits double. This limits how many integers we can handle. Do you know how terra deals with that ?

[bfrank5]

Hey JR,

I think the clipping is a workaround, but I was having trouble with parallel processing in catalog_apply using terra::crop. I run into a std::bad_alloc error, but I don't have a lot of specific information yet (nor anything reproducible). I could potentially tile the dtm into separate files and dynamically fetch the file path in catalog_apply, but trying a few other ideas first.

(Edit: I resolved the std::bad_alloc error, it turned out I was clipping with the wrong bounding box)

I am playing with simply removing the as.integer() call right now using a fork of the repository, and it appears to work fine, obviously you will want a more rigorous solution. I am still checking a larger dataset now.

As far as I can tell, at least for my application using normalize_heights, the next step is here:

lidR/R/utils_raster.R

Lines 94 to 121 in 09a878b

	raster_value_from_cells = function(raster, cells, layer = 1)
	{
	if (inherits(raster, "Raster"))
	{
	if (raster_nlayer(raster) > 1)
	return(raster::extract(raster, cells, layer = layer, nl = 1)[,1])
	else
	return(raster::extract(raster, cells, layer = layer, nl = 1))
	}
	if (is(raster, "stars_proxy"))
	stop("stars_proxy not yet supported in 'raster_value_from_cells()'", call. = FALSE) # nocov
	if (is(raster, "stars"))
	{
	dims <- dim(raster)
	if (is.na(dims[3]))
	return(raster[[1]][cells])
	return(raster[,,,layer][[1]][cells])
	}
	if (is(raster, "SpatRaster"))
	return(terra::extract(raster, cells)[[layer]])
	raster_error() # nocov
	}

where cells is the vector if indices attained previously. Drilling down, the terra::crop implementation ultimately goes here:

https://github.com/rspatial/terra/blob/e5458c22fc97a51bb73523ac8c971e975cda0afa/R/extract_single.R#L138-L155

I think the operable part here is rowColFromCell, which you can see the typing for here in the header file here:

https://github.com/rspatial/terra/blob/e5458c22fc97a51bb73523ac8c971e975cda0afa/src/spatRaster.h#L426

and appears to take a double. I can't speak to the other methods (e.g., raster::extract).

I wonder if it is sufficient to check if it is a whole number? And even then, terra::crop rounds the indices initially to the nearest whole number.

Up to you! I mostly wanted to stick this here in case anyone else runs into the same problem.