rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
538 stars 89 forks source link

Random segfaults using terra::focal on big rasters #1563

Open vidlb opened 3 months ago

vidlb commented 3 months ago

Hi there,

I'm trying to apply a focal function (sum) to a big raster input / dtm. I sometimes face segmentation faults, but it doesn't seems to be a problem with the data since it may complete successfully with another attempt.

This is quite annoying since the process takes several hours and it seems it often occurs near the end. The final output can reach from 12 to 18GB (compressed), with a size up to 100 000 x 100 000 pixels. Here is a screenshot of the error : rcpp_error

And here is the code.

smooth <- function(inpath, outpath) {
  ######## Mean smooth DTM in a 25m window ##########
  dtm <- terra::rast(inpath)
  # Prepare output
  nodata <- -99999
  # Create kernel
  radius <- 12.5
  kernel <- terra::focalMat(dtm, radius, "circle")
  # Apply focal function and write output
  copt <- c("COMPRESS=ZSTD", "ZLEVEL=1", "PREDICTOR=3", "TILED=YES", "BIGTIFF=IF_SAFER")
  opts <- list(gdal = copt, NAflag = nodata)
  terra::focal(dtm, kernel, "sum", filename = outpath, overwrite = TRUE, wopt = opts)
}

I'm running the script in a docker container based on the rocker/geospatial:4.4 image, terra version is 1.7.78. Any ideas ? Is my raster just too big ? It's a VRT pointing to a lot of files. Or should I use a different setting for na.policy than the default one ? The error message indeed contains NA related functions calls. Yet I don't understand why it doesn't always occur.

vidlb commented 3 months ago

FYI I also tried to use focalCpp without success, I tried to adapt the weighted sum function from the docs (just needed a sum), then tried to literally copy paste the example from the docs but it always failed (verbose compiled but then failed to run with a message like "test failed").

mdsumner commented 3 months ago

How big is kernel? Are you possibly asking for a radius in metres in a longlat context?

vidlb commented 3 months ago

The kernel radius is 12.5, and my CRS is in meters. Resolution is 1m so this is a 25*25px window.