Cropping fails when attempting to crop more than 65535 layers at once

[ErikKusch]

Hiya,

trying to handle multi-layer SpatRasters and I am running into an issue when cropping the data. This error occurs as soon as the layer count exceeds 65535.

Minimal Working Example

Below is a minimal working example adapted from the terra::crop() documentation:

library(terra)

## data to crop
f <- system.file("ex/elev.tif", package="terra")
rr <- rast(f)
r <- c(rep(rr, 65535+1)) # repeat base layer one too many times for cropping

## shape to crop with
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)

## cropping that works
cm <- crop(r[[-1]], v, mask = TRUE, touches = TRUE) # reduce layer count by one

## cropping that fails
cm <- crop(r, v, mask = TRUE, touches = TRUE)
#Error: [crop] failed writing GTiff file
#In addition: Warning message:
#/private/var/folders/c9/5tlmp10s517_f5qmn120ctjc0000gp/T/RtmpUJljN2/spat_9a79465d898_39545_2.tif: Attempt to #create 94x88x65536 TIFF file, but bands must be lesser or equal to 65535. (GDAL error 1) 

Proposed Solution

How about adding a check to the cropping (and probably also masking) function and carry out cropping on maximally stacked SpatRasters and fusing them re-stacking them after cropping has been done? Like so:

library(terra)
library(sf)

## data to crop
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
r <- c(rep(r, ceiling(65535*1.7))) # a lot more layers than crop can handle at once

## shape to crop with
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)

## identifying layer indices for splitting into SpatRasters of maximum 65535 layers
Indices <- ceiling((1:nlyr(r))/65535)

## splitting SpatRaster into list of SpatRasters
r_ls <- terra::split(x = r, f = Indices)

## Cropping list of SpatRasters
cm_ls <- lapply(r_ls, terra::crop, y = v, mask = TRUE, touches = TRUE)

## Fusing list of cropped SpatRasters back together
cm <- do.call(c, cm_ls)

When checking that this cropping worked as intended, I find no issue:

all.equal(crop(r[[1]], v, mask = TRUE, touches = TRUE), cm[[1]]) [1] TRUE

Maybe this could even be combined with the addition of a cores argument akin to the implementation in the tapp function.

Operating System & Package Versions

I am on MacOS Sonoma (v14.5).

packageVersion("terra") [1] ‘1.7.78’

terra::gdal(lib="all") gdal proj geos "3.5.3" "9.1.0" "3.11.0"

[rhijmans]

I will look into using a different file format for temporary files, that can have more than than 65535 layers. Your solution can work, but I need something that automatically works (in the C++ code) for all raster methods.

[rhijmans]

I "fixed" this in the sense that you now get an error message when terra attempts to write a (temp) file with more than 65,535 layers.

Error: [writeRaster] cannot write more than 65535 layers

I should eventually implement a solution that can handle this, but that may take a while.