Split and re-merge RasterLayer(s)

splitRaster divides up a raster into an arbitrary number of pieces (tiles). Split rasters can be recombined using do.call(merge, y) or mergeRaster(y), where y <- splitRaster(x).

mergeRaster(x, fun = NULL)

# S4 method for list
mergeRaster(x, fun = NULL)

splitRaster(
  r,
  nx = 1,
  ny = 1,
  buffer = c(0, 0),
  path = NA,
  cl,
  rType = "FLT4S",
  fExt = ".tif"
)

Arguments

x

A list of split raster tiles (i.e., from splitRaster).

fun

Function (e.g. mean, min, or max that accepts a na.rm argument. The default is mean.

r

The raster to be split.

nx

The number of tiles to make along the x-axis.

ny

The number of tiles to make along the y-axis.

buffer

Numeric vector of length 2 giving the size of the buffer along the x and y axes. If values greater than or equal to 1 are used, this is interpreted as the number of pixels (cells) to use as a buffer. Values between 0 and 1 are interpreted as proportions of the number of pixels in each tile (rounded up to an integer value). Default is c(0, 0), which means no buffer.

path

Character specifying the directory to which the split tiles will be saved. If missing, the function will write to memory.

cl

A cluster object. Optional. This would generally be created using parallel::makeCluster() or equivalent. This is an alternative way, instead of beginCluster(), to use parallelism for this function, allowing for more control over cluster use.

rType

Data type of the split rasters. Defaults to FLT4S.

fExt

file extension (e.g., ".grd" or ".tif") specifying the file format.

Value

mergeRaster returns a RasterLayer object.

splitRaster returns a list (length nx*ny) of cropped raster tiles.

Details

mergeRaster differs from merge in how overlapping tile regions are handled: merge retains the values of the first raster in the list. This has the consequence of retaining the values from the buffered region in the first tile in place of the values from the neighbouring tile. On the other hand, mergeRaster retains the values of the tile region, over the values in any buffered regions. This is useful for reducing edge effects when performing raster operations involving contagious processes.

This function is parallel-aware using the same mechanism as used in raster: NOTE: This may not work as expected as we transition away from raster. Specifically, if you start a cluster using raster::beginCluster(), then this function will automatically use that cluster. It is always a good idea to stop the cluster when finished, using raster::endCluster().

See also

terra::merge(), terra::mosaic()

do.call(), terra::merge().

Author

Yong Luo, Alex Chubaty, Tati Micheletti & Ian Eddy

Alex Chubaty and Yong Luo

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)
set.seed(1462)

## an example with dimensions: nrow = 77,  ncol = 101, nlayers = 3
b <- rast(system.file("ex/logo.tif", package = "terra"))
r <- b[[1]] # use first layer only
nx <- 3
ny <- 4

tmpdir <- dir.create(file.path(tempdir(), "splitRaster-example"), showWarnings = FALSE)

y0 <- splitRaster(r, nx, ny, path = file.path(tmpdir, "y0")) # no buffer

## buffer: 10 pixels along both axes
y1 <- splitRaster(r, nx, ny, c(10, 10), path = file.path(tmpdir, "y1"))

## buffer: half the width and length of each tile
y2 <- splitRaster(r, nx, ny, c(0.5, 0.5), path = file.path(tmpdir, "y2"))

## the original raster:
if (interactive()) plot(r) # may require a call to `dev()` if using RStudio

## the split raster:
layout(mat = matrix(seq_len(nx * ny), ncol = nx, nrow = ny))
plotOrder <- unlist(lapply(split(1:12, rep(1:nx, each = ny)), rev))

if (interactive()) {
  invisible(lapply(y0[plotOrder], terra::plot))
}

## parallel splitting
if (requireNamespace("raster", quietly = TRUE) &&
    requireNamespace("parallel")) {
  if (interactive()) {
    n <- pmin(parallel::detectCores(), 4) # use up to 4 cores
    raster::beginCluster(n, type = "PSOCK")
    y3 <- splitRaster(r, nx, ny, c(0.7, 0.7), path = file.path(tmpdir, "y3"))
    raster::endCluster()
    if (interactive()) {
      invisible(lapply(y3[plotOrder], terra::plot))
    }
  }
}

## can be recombined using `terra::merge`
m0 <- do.call(merge, y0)
all.equal(m0, r) ## TRUE
#> [1] TRUE

m1 <- do.call(merge, y1)
all.equal(m1, r) ## TRUE
#> [1] TRUE

m2 <- do.call(merge, y2)
all.equal(m2, r) ## TRUE
#> [1] TRUE

## or recombine using mergeRaster
n0 <- mergeRaster(y0)
all.equal(n0, r) ## TRUE
#> [1] TRUE

n1 <- mergeRaster(y1)
all.equal(n1, r) ## TRUE
#> [1] TRUE

n2 <- mergeRaster(y2)
all.equal(n2, r) ## TRUE
#> [1] TRUE

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)
unlink(tmpdir, recursive = TRUE)