Open source and S4 extensible framework for efficient spatial and pixel-based tiling operations on large datasets. {tilework} provides easy-to-use tile planners that enable memory-efficient processing of spatial data through parallelizable tile-based operations.
For another approach to spatially tiled computation, see: chopin
- Flexible Tiling: Support for spatial extent-based, pixel-exact, arbitrary point-centered, and adaptive variable-size tiling
- Memory Efficient: Process tilewise or batchwise without loading entire datasets into memory
- Stateful Iteration: Iterator patterns for streaming and batch processing
- Parallel Processing: Built-in support for parallel execution via the {future} framework
- Flexible Padding: Add padding around tiles to handle edge effects
- Metadata Support: Attach custom metadata to tiles for advanced workflows
- Terra Integration: Seamless integration with the {terra} package for spatial data handling
# Install from GitHub
devtools::install_github("drieslab/tilework")library(tilework)
library(terra)
# Load a raster
f <- system.file("ex/elev.tif", package = "terra")
r <- rast(f)
# Create a spatial tile iterator
tp <- spatialTilePlan(ext = ext(r), n = 16)
# Check tile layout
tp
dim(tp)
plot(tp)
# Extract specific tiles
tile_ext <- tp[5] # Get 5th tile
tile_grid <- tp[1, 2:3] # Get specific grid positions
# Apply a function across tiles
outdir <- tempdir()
tileApply(r, tiles = tp, FUN = function(x, .I) {
writeRaster(x, file.path(outdir, sprintf("tile_%03d.tif", .I)))
})# Create a pixel-based tile iterator
px <- pixelTilePlan(pxdims = c(500, 500), ncols = 100, nrows = 100)
# Check dimensions
dim(px) # Grid dimensions
length(px) # Total number of tiles
plot(px) # Visualize grid
# Apply processing with pixel tiles
tileApply(r, tiles = px, FUN = function(x) {
# Process each 100x100 pixel tile
mean(values(x), na.rm = TRUE)
})tilePlan
Virtual base class for all tile iterators with common functionality:
- Tile indexing with
[i]and[i,j]notation - Padding with
+and-operators - Metadata management with
$accessor - Plotting capabilities
spatialTilePlan
For spatial extent-based tiling:
- Define tiles using geographic coordinates
- Automatic grid layout optimization
- {terra}
SpatExtentintegration
tp <- spatialTilePlan(ext = c(0, 100, 0, 100), n = 9)
dim(tp) # Returns [3, 3] - actual grid layoutpixelTilePlan
For pixel-exact tiling:
- Define tiles using pixel coordinates
- Precise control over tile dimensions
- Ideal for image processing workflows
pti <- pixelTilePlan(pxdims = c(1000, 1000), ncols = 250, nrows = 250)pointTilePlan
For tiling centered on arbitrary (x, y) coordinates with uniform tile dimensions.
Tile placement is driven by supplied point locations rather than a grid.
Supports spatial (CRS units) or pixel coordinate modes via input / output
toggles; cross-mode conversion is resolved automatically at extraction time when
a raster is provided. See the
tile plans vignette
for full coverage of coordinate modes.
tp <- pointTilePlan("spatial",
coords = cbind(x = c(10, 50, 90), y = c(10, 50, 90)),
width = 20,
height = 20
)freeTilePlan
For explicit per-tile bounds with no required uniformity in size or spacing.
Tile bounds are the canonical representation — positions are not computed from
a formula. The primary use case is adaptive decomposition via quadtreePlan(),
where high-density regions use small tiles and sparse regions use large tiles.
# Manual bounds (e.g. from an external partitioning algorithm)
tp <- freeTilePlan()
tp$bounds <- rbind(
c(0, 50, 0, 50),
c(50, 100, 0, 50),
c(0, 50, 50, 100),
c(50, 100, 50, 100)
)
length(tp) # 4
tp[2] # SpatExtent for tile 2
plot(tp)quadtreePlan() builds a freeTilePlan automatically by iteratively
subdividing tiles whose FUN value exceeds a threshold, then merging
neighboring leaf tiles back together when their combined value stays ≤
threshold. The last FUN value per leaf is stored in $n_records.
# Points on disk (required for tileApply dispatch)
pts <- terra::vect(f, proxy = TRUE)
fp <- quadtreePlan(
x = pts,
threshold = 500L,
min_tile_size = 1
)
plot(fp)
fp$n_records # point count per leaf tiletileGroup
Organize tiles groups. These are created on top of tilePlan
classes. These serve as lazy selection(s) of particular tiles of the
underlying tilePlan.
# Create groups that should be processed separately or differently
tg <- tileGroup(tp, groups = list(
"quadrant1" = 1:4, # First 4 tiles
"quadrant2" = c(5, 6, 9, 10), # Specific tiles
"border" = list(c(1,3), c(1,3)) # Grid-based selection
))
# Set active group for easy access
tg$active <- "quadrant1"
length(tg) # Returns length of active group
tg[, 2] # Second tile from active grouptileIterator
Stateful iterator for streaming processing. Created on top of tilePlan or
tileGroup (with $active set). Use with tileApply() for distribution of
batches across parallelized {future} workers. A setup_FUN argument initializes
per-worker state (e.g. loading a model) once before batch processing begins —
see the ML vignette
for a worked example.
# Create iterator for batch processing
iter <- tileIterator(tp, batch_size = 3)
# Check status
iter$has_next
iter$remaining
iter$progress
# Process in batches
while (iter$has_next) {
batch <- iter$next_batch()
cat("Processing", length(batch), "tiles\n")
}
# Reset for another pass
iter$reset()tileSelection
Lazy drop = FALSE selection wrapper — preserves a subset of tile indices without materialising bounds. Useful for selecting specific tiles to process without modifying the underlying plan.
# Select tiles of interest (e.g. from a spatial query)
sel <- tp[c(1, 3, 7), drop = FALSE]
length(sel) # 3
tileApply(r, tiles = sel, FUN = function(x) terra::global(x, "mean"))# Load a raster file
f <- system.file("ex/elev.tif", package="terra")
r <- terra::rast(f)
# Create tile plan matching raster
tp <- pixelTilePlan(pxdims = dim(r)[1:2], nrows = 100, ncols = 100)
# Extract tiles
tiles <- getTile(r, tp, i = 1:4) # Get first 4 tiles
tile_data <- getTile(r, tp, i = 3, j = 5) # Get specific grid positionApply functions across tiles with automatic parallellization.
# Process tiles in parallel
results <- tileApply(r, tiles = tp, FUN = function(x, .I) {
# x is the tile raster data
# .I is the tile number
# Example: calculate mean value per tile
terra::global(x, "mean", na.rm = TRUE)
})
# Save tiles to files
outdir <- tempdir()
tileApply(r, tiles = tp, FUN = function(x, .I, .R, .C) {
# .R and .C provide row/col indices
filename <- file.path(outdir, sprintf("tile_r%d_c%d.tif", .R, .C))
terra::writeRaster(x, filename)
})Depending on the class of tp, different parallelization schemes are used.
# Process different groups with different strategies
tileApply(r, tiles = tg,
parallel_strategy = "groups", # Parallelize across groups
FUN = function(x, .GROUP) {
# Process based on group
if (.GROUP == "border") {
# Special processing for border tiles
terra::focal(x, w = matrix(1/9, 3, 3))
} else {
# Standard processing
x
}
}
)Add padding around tiles to handle edge effects:
# Add 10-unit padding to all tiles
padded_ti <- tp + 10
# Remove 5-unit padding
reduced_ti <- tp - 5
# Preview padded tiles
plot(padded_ti, alpha = 0.3)Attach custom metadata to tiles:
# Add metadata
tp$processing_priority <- sample(1:3, length(tp), replace = TRUE)
tp$output_format <- "GTiff"
# Access metadata
tp$processing_priority
tp[[5]] # Metadata for tile 5
# Metadata is preserved during tile selection
tile_data <- tp[5]
attr(tile_data[[1]], "processing_priority")Split iterators. This pattern is what powers the tileApply() method for tileIterator
# Create base iterator
iter <- tileIterator(tp, batch_size = 5)
# Split across 4 workers
worker_iters <- iterSplit(iter, n = 4, distribute = TRUE)
# Each worker gets independent iterator with subset of tiles
sapply(worker_iters, function(x) x$remaining)Built-in support for parallel execution:
library(future)
plan(multisession, workers = 4)
# Parallel tile processing
results <- tileApply(r, tiles = tp,
FUN = function(x) {
# Your processing function
mean(values(x), na.rm = TRUE)
})Multiple ways to access tiles:
# Single tile by index
tile_5 <- tp[5]
# Multiple tiles
tiles_1_to_3 <- tp[1:3]
# Grid-based selection (row, column)
corner_tiles <- tp[1, c(1, ncol(tp))]
# All tiles
all_tiles <- tp[]- Memory Management: Use appropriate tile sizes to balance memory usage and processing efficiency
- Pad Planning: Consider padding requirements for spatial operations to avoid edge effects
- Parallel Strategy: Choose between parallelizing across groups vs. within groups based on your workflow — see the decision table in the orchestration vignette
- Metadata Usage: Leverage metadata for complex processing logic and file organization
- Iterator Patterns: Use stateful iterators for streaming large datasets that don't fit in memory
Processing Large Satellite Images
# Load large satellite image
large_raster <- rast("large_satellite_image.tif")
# Create efficient tiling scheme
tp <- spatialTilePlan(ext = ext(large_raster), n = 100)
# Add padding for edge effects
tp <- tp + 50 # 50-unit padding
# Process tiles in parallel
plan(multisession, workers = 8)
results <- tileApply(large_raster, tiles = tp,
FUN = function(x, .I) {
# Apply NDVI calculation
ndvi <- (x[[4]] - x[[3]]) / (x[[4]] + x[[3]])
# Save processed tile
writeRaster(ndvi,
sprintf("ndvi_tile_%03d.tif", .I))
# Return summary statistics
c(mean = mean(values(ndvi), na.rm = TRUE),
sd = sd(values(ndvi), na.rm = TRUE))
})# High-resolution image processing
image <- rast("high_res_image.tif")
# Create pixel-exact tiles
pti <- pixelTilePlan(pxdims = c(nrow(image), ncol(image)), ncols = 512, nrows = 512)
# Process each tile
texture_metrics <- tileApply(image, tiles = pti,
FUN = function(x) {
# Calculate texture metrics
vals <- values(x)
list(
contrast = var(vals, na.rm = TRUE),
homogeneity = 1 / (1 + var(vals, na.rm = TRUE))
)
})Dependencies
- terra: Spatial data handling and raster operations
- checkmate: Input validation
- future.apply: Parallel processing support
- Choosing and Creating a Tile Plan — when to use
spatialTilePlan,pixelTilePlan,pointTilePlan, orfreeTilePlan; adaptive quadtree decomposition; coordinate modes; padding - Selecting, Grouping, and Iterating Tiles —
tileSelection,tileGroupparallelization strategies,tileIteratorstreaming and per-worker setup - Patch-Based Feature Extraction for Machine Learning — end-to-end ML inference pipeline using
tileGroup,tileIterator, andsetup_FUN
Contributions are welcome! Please feel free to submit issues, feature requests, or pull requests.