vrtstack writes GDAL multidimensional VRTs: small XML documents that let a collection of rasters be read as one array through any GDAL multidim reader, with no data copied. There are three ways in, for three situations.
The original job: you have 2D rasters (COGs, GeoTIFFs) that are
slices of one variable, and you want them as a single 3D array.
gdal_mdim_stack() builds the new axis from scratch.
# soil layers by depth
gdal_mdim_stack(
files = c("soil_0cm.tif", "soil_30cm.tif", "soil_100cm.tif"),
dim_name = "depth",
dim_values = c(0, 30, 100),
dim_type = "VERTICAL",
dim_direction = "DOWN",
dim_unit = "cm",
output = "soil.vrt"
)Evenly spaced values are detected and stored compactly as
<RegularlySpacedValues> rather than a long value
list; dimension values can be pulled from file names with a one-group
regex and parsed as dates; and http(s):// sources are
wrapped in /vsicurl/ automatically.
# a daily time series, dates parsed from the file names
gdal_mdim_stack(
files = list.files(pattern = "temp_\\d{8}\\.tif$"),
dim_name = "time",
dim_values = "(\\d{8})",
parse_format = "%Y%m%d",
dim_type = "TEMPORAL",
time_origin = "1970-01-01",
time_unit = "days",
output = "temperature.vrt"
)The convenience wrappers vrt_time_stack(),
vrt_depth_stack() and vrt_regular_stack()
cover the common cases with shorter calls.
Each source carries a <SourceSlab> with an
explicit count, so a read through the VRT opens only the
sources it covers rather than every file.
Some sources are already multidimensional (NetCDF, Zarr, HDF5) and may hold several arrays at once. For these, the fastest route is to take a small “first few files” VRT as a recipe and flesh it out to the full file set, re-writing only the concat dimension and the per-array source lists. Everything else - coordinate arrays, slab geometry, per-array scale/offset/unit/attributes - is carried through verbatim.
# `template = TRUE` auto-generates the recipe from the first n files via
# `gdal mdim mosaic`, then extends it to the whole set.
vrtstack(
files = all_files, # e.g. thousands of daily NetCDFs
concat = "(\\d{8})", parse_format = "%Y%m%d",
origin = "1978-01-01", unit = "days",
template = TRUE, concat_dim = "time",
output = "series.vrt"
)The concat dimension need not be the first axis, nor named
"time": the axis position is read from the recipe. A
ready-made recipe can be passed instead of TRUE (a path or
an XML string).
metadata_from chooses which source donates the
coordinate units and root-group provenance restored onto an
auto-generated recipe ("first", "last", or an
explicit path) - a repair for metadata that
gdal mdim mosaic drops. graft_coord_units()
and graft_group_attrs() expose that repair directly and are
no-ops once the metadata is already present.
When the grid is known from a catalog (STAC, CMR, or a URL template) rather than by opening a file, the emitter path builds the VRT from two plain structures and never touches the data: a sources table (one row per step) and a signature (the grid and array metadata for the collection).
src <- mdim_sources(
time = as.Date(c("2024-01-01", "2024-01-02")),
dsn = c("/vsicurl/.../20240101.tif", "/vsicurl/.../20240102.tif"),
source_ref = "1", kind = "band", # a COG band; use "array" for NetCDF
source_view = "[newaxis,...]" # rank a 2D band up to the 3D cube
)
sig <- mdim_signature(
collection = "demo", variable = "sst",
nx = 1440, ny = 720,
x0 = 0.125, dx = 0.25, y0 = 89.875, dy = -0.25,
srs = "EPSG:4326", dtype = "Int16",
nodata = -999, scale = 0.01, unit = "Celsius"
)
emit_mdim_vrt(src, sig, "demo.vrt")The contract is “identical grid or refuse”: the emitter trusts the
signature, so a collection whose assets are not on the same grid is a
refusal, never a resample. stac_sufficiency() reports
whether a parsed STAC item carries the fields the signature needs, and
verify_vrt() round-trips an emitted VRT through
gdal mdim info.
The pure-XML paths - gdal_mdim_stack(),
mdim_template_extend() and emit_mdim_vrt()
writing the document - work with any GDAL that gdalraster links. The
paths that drive GDAL’s unified command-line algorithms through
gdalraster::gdal_run() - verify_vrt(),
template = TRUE recipe generation, and the metadata grafts
- require GDAL >= 3.11.3.
The VRT opens through any GDAL multidim reader:
gdalmdiminfo and gdalmdimtranslate on the
command line, gdalraster in R, or an xarray GDAL backend in
Python. Selecting along the new dimension reads only the sources that
the selection covers.