Package 'vrtstack'

Emit a multidimensional VRT from a sources table and a signature

Description

Writes a GDAL multidimensional VRT for one variable with dimensions ⁠(time, Y, X)⁠. Y/X are ⁠<RegularlySpacedValues>⁠ from the signature; the concat coordinate is written as inline string values from the sources table. One ⁠<Source>⁠ is emitted per file: a multi-step file contributes a single source whose full extent lands at the DestSlab offset of its first step.

Usage

emit_mdim_vrt(sources, sig, path)

Arguments

sources	A sources table from mdim_sources().
sig	A signature from mdim_signature().
path	Output VRT path.

Value

path, invisibly.

See Also

verify_vrt()

---

Create a Multidimensional VRT from 2D Rasters

Description

Stack 2D rasters along a new dimension, producing a GDAL Multidimensional VRT. This enables treating a collection of 2D rasters (COGs, GeoTIFFs, etc.) as a single multidimensional array, readable by GDAL's multidim API and xarray backends.

Usage

gdal_mdim_stack(
  files,
  dim_name,
  dim_values = NULL,
  dim_start = NULL,
  dim_increment = NULL,
  dim_type = NULL,
  dim_direction = NULL,
  dim_unit = NULL,
  time_origin = NULL,
  time_unit = "days",
  parse_format = NULL,
  array_name = NULL,
  output = NULL
)

Arguments

files	Character vector of file paths (COGs, GeoTIFFs, /vsicurl/ URIs, etc.)
dim_name	Name for the new dimension (e.g., "time", "depth", "level")
dim_values	Dimension values. One of: Numeric vector of values (length must match files) A regex pattern with one capturing group to extract from filenames A function taking filename (basename) and returning a numeric value NULL if using dim_start + dim_increment
dim_start, dim_increment	For regular axes, specify start and increment instead of explicit values. Files are assumed to be in the correct order.
dim_type	Optional dimension type: "TEMPORAL", "VERTICAL", "PARAMETRIC"
dim_direction	Optional direction: "NORTH", "SOUTH", "EAST", "WEST", "UP", "DOWN", "FUTURE", "PAST"
dim_unit	Unit string for the dimension (e.g., "days since 1970-01-01", "cm")
time_origin	For temporal dimensions, origin date as "YYYY-MM-DD" or POSIXct. When set with parse_format, values are encoded as offsets from this origin.
time_unit	For temporal dimensions: "seconds", "minutes", "hours", "days"
parse_format	For date/time extraction, a strptime() format string (e.g., "%Y%m%d" for "20230115")
array_name	Name for the data array in the VRT. Default: derived from common filename prefix, or "data".
output	Output VRT path. If NULL, returns XML string without writing.

Details

The function generates a GDAL Multidimensional VRT that virtually stacks 2D rasters along a new dimension. Key features:

• Regular axes: When dimension values are evenly spaced (or specified via dim_start/dim_increment), the VRT uses ⁠<RegularlySpacedValues>⁠ for compact representation.

• Auto-detection: Even with explicit dim_values, regular spacing is detected automatically and converted to the compact form.

• Filename parsing: Extract dimension values from filenames using regex, with optional date/time parsing. Useful for time series like "temperature_20230101.tif".

• Cloud-native: Supports ⁠/vsicurl/⁠ URIs for remote files. HTTP(S) URLs are automatically prefixed with ⁠/vsicurl/⁠.

The resulting VRT can be read by:

• GDAL's multidim API: gdalmdiminfo, gdalmdimtranslate

• Python xarray backends: gdx, xgdal

• Any tool supporting GDAL multidimensional rasters

Value

Path to VRT file (invisibly), or XML string if output is NULL.

Examples

## Not run: 
# Stack 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_stack.vrt"
)

# Time series with regex extraction
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"
)

# Regular axis - monthly data
gdal_mdim_stack(
  files = sprintf("precip_%02d.tif", 1:12),
  dim_name = "month",
  dim_start = 1,
  dim_increment = 1,
  output = "monthly.vrt"
)

# Cloud-hosted COGs
gdal_mdim_stack(
  files = paste0(
    "https://storage.googleapis.com/solus100pub/claytotal_",
    c(0, 5, 15, 30, 60, 100, 150), "_cm_p.tif"
  ),
  dim_name = "depth",
  dim_values = c(0, 5, 15, 30, 60, 100, 150),
  dim_type = "VERTICAL",
  dim_unit = "cm",
  array_name = "claytotal",
  output = "solus_clay.vrt"
)

## End(Not run)

---

Inject <Unit> on coordinate arrays that lost it in gdal mdim mosaic

Description

A coordinate array is any <Array> with no <Source>. Units are read from one source via gdal mdim info (or supplied directly via units). Arrays that already carry a <Unit> are left untouched, so this is idempotent and a no-op once GDAL propagates the unit itself.

Usage

graft_coord_units(x, source = NULL, units = NULL)

Arguments

x	Mosaic recipe: xml2 document or XML string.
source	One source path to read units from (via gdal mdim info).
units	Optional named character vector (array name -> unit); bypasses source and needs no GDAL (handy for tests).

Value

Enriched x as an xml_document.

---

Inject invariant root-group provenance attributes dropped by mosaic

Description

Inject invariant root-group provenance attributes dropped by mosaic

Usage

graft_group_attrs(
  x,
  source = NULL,
  attrs = NULL,
  skip = c("Conventions", .MDIM_PER_FILE_GLOBALS)
)

Arguments

x	Mosaic recipe: xml2 document or XML string.
source	One source path to read globals from (via gdal mdim info).
attrs	Optional named character vector (name -> value); bypasses source (handy for tests). All treated as String attributes.
skip	Attribute names to exclude (default: Conventions + per-file).

Value

Enriched x as an xml_document.

---

Construct a grid signature for a multidimensional VRT

Description

One row (as a list) per collection: the grid geometry and array metadata the emitter needs, derived from catalog metadata or a single cached probe. The emitter trusts these values, so a collection with non-aligned assets is a signature-gate refusal, never a resample.

Usage

mdim_signature(
  collection,
  variable,
  nx,
  ny,
  x0,
  dx,
  y0,
  dy,
  srs,
  dtype,
  nodata = NA,
  scale = NA,
  offset = NA,
  unit = NA_character_,
  steps_per_file = 1L,
  chunk_nx = NA_integer_,
  chunk_ny = NA_integer_,
  provenance = "probe"
)

Arguments

collection, variable	Identifiers for the collection and data array.
nx, ny	Grid dimensions in X and Y.
x0, dx, y0, dy	Grid origin and spacing as pixel-centre coordinates in the array's own orientation (dy is typically negative for north-up grids).
srs	Spatial reference (e.g. "EPSG:4326" or a WKT string).
dtype	VRT data type of the array (e.g. "Int16", "Float32").
nodata, scale, offset, unit	Optional array metadata; NA to omit.
steps_per_file	Number of concat-axis steps stored in each source file.
chunk_nx, chunk_ny	Optional source chunk (block) shape, recorded for later byte-range harvesting.
provenance	Free text recording where the grid facts came from (catalog, probe, or a published grid).

Value

A named list describing one collection's grid and array metadata.

See Also

mdim_sources(), emit_mdim_vrt()

---

Construct a sources table for a multidimensional VRT

Description

One row per output step along the concat (time) axis. This is the interface contract produced by adapters and consumed by emit_mdim_vrt().

Usage

mdim_sources(
  time,
  dsn,
  source_ref,
  kind,
  step_in_file = 0L,
  dest_index = NULL,
  source_view = NA_character_
)

Arguments

time	Vector of per-step coordinate values (e.g. a Date vector).
dsn	Source data source names (file paths or ⁠/vsicurl/⁠ URIs).
source_ref	Array path (e.g. "/sst") when kind = "array", or band number (e.g. "1") when kind = "band".
kind	Source flavour: "array" (⁠<SourceArray>⁠, a multidimensional source) or "band" (⁠<SourceBand>⁠, a classic 2D source).
step_in_file	Zero-based step index within a multi-step file.
dest_index	Zero-based destination position along the concat axis; defaults to the rank order of time.
source_view	Optional GetView expression that brings each source to exactly the destination rank (e.g. "[newaxis,...]" to rank a 2D band up to a 3D cube, "[:,0,:,:]" to drop a singleton dimension). NA for none.

Value

A data frame ordered by dest_index, one row per output step.

See Also

mdim_signature(), emit_mdim_vrt()

---

Extend an mdim VRT recipe to a full file set

Description

Pure XML-tree rewrite (no GDAL needed). Sets the concat dimension size, rewrites the concat coordinate array's values, and regenerates the ⁠<Source>⁠ list of every data array (an ⁠<Array>⁠ that owns a ⁠<Source>⁠). Per-array SourceArray path and slab geometry are read from the array's first source, never assumed from the array name.

Usage

mdim_template_extend(xml, files, concat_values = NULL, concat_dim = "time")

Arguments

xml	Template VRT as an XML string (or anything xml2::read_xml() takes).
files	Full character vector of sources.
concat_values	Values for the concat coordinate array, or NULL to keep the template's values (only valid when length matches files).
concat_dim	Concat dimension name (default "time"); if not found, the dimension whose size equals the template's per-array source count is used.

Value

An xml_document (serialise with as.character() / xml2::write_xml()).

---

Report mdim-VRT metadata sufficiency of a STAC item

Description

The emitter's requirements expressed as a checklist over one already-parsed STAC item (as from jsonlite::fromJSON(..., simplifyVector = FALSE)). No HTTP is performed here; fetching items is an adapter's job.

Usage

stac_sufficiency(item, asset = NULL)

Arguments

item	A parsed STAC item (a list).
asset	Optional asset key; when supplied, asset-level fields take precedence over item-level properties.

Details

proj:epsg versus proj:code is a proj-extension version marker, not a gap: either satisfies the CRS requirement. datetime may be null for interval items, with start_datetime/end_datetime carrying the extent.

Value

A data frame with one row per checked field: field, found_at (one of "item-properties", "asset", "absent"), and a short value_preview.

---

Verify a multidimensional VRT with ⁠gdal mdim info⁠

Description

Opens path through gdalraster's ⁠mdim info⁠ algorithm as a round-trip check that the emitted VRT parses and can be interpreted by GDAL.

Usage

verify_vrt(path)

Arguments

path	Path to a VRT.

Value

The ⁠gdal mdim info⁠ output as a character string, invisibly.

---

Stack Rasters by Depth/Level into Multidimensional VRT

Description

A convenience wrapper for gdal_mdim_stack() for vertical profile data (soil layers, atmospheric levels, ocean depths, etc.)

Usage

vrt_depth_stack(
  files,
  values,
  direction = "DOWN",
  unit = NULL,
  dim_name = NULL,
  array_name = NULL,
  output = NULL
)

Arguments

files	Character vector of file paths
values	Numeric vector of depth/level values, or regex pattern to extract
direction	Vertical direction: "DOWN" (depth increases), "UP" (elevation). Default: "DOWN"
unit	Unit string (e.g., "m", "cm", "hPa")
dim_name	Dimension name. Default: "depth" if direction is DOWN, "level" if UP.
array_name	Name for the data array
output	Output VRT path

Value

Path to VRT file (invisibly), or XML string if output is NULL.

Examples

## Not run: 
# Soil layers at specific depths
vrt_depth_stack(
  files = c("soil_0cm.tif", "soil_30cm.tif", "soil_100cm.tif"),
  values = c(0, 30, 100),
  unit = "cm",
  output = "soil_profile.vrt"
)

# Extract depth from filename
vrt_depth_stack(
  files = list.files(pattern = "temp_\\d+m\\.tif$"),
  values = "(\\d+)m",  # extracts "50" from "temp_50m.tif"
  unit = "m",
  output = "ocean_profile.vrt"
)

## End(Not run)

---

Stack Rasters with Regular Spacing

Description

A convenience wrapper for gdal_mdim_stack() when files represent regularly spaced values (e.g., monthly data, hourly snapshots, uniform vertical levels).

Usage

vrt_regular_stack(
  files,
  dim_name,
  start,
  increment,
  unit = NULL,
  array_name = NULL,
  output = NULL
)

Arguments

files	Character vector of file paths, in order
dim_name	Name for the new dimension
start	Starting value
increment	Spacing between values
unit	Unit string
array_name	Name for the data array
output	Output VRT path

Value

Path to VRT file (invisibly), or XML string if output is NULL.

Examples

## Not run: 
# 12 monthly files
vrt_regular_stack(
  files = sprintf("precip_%02d.tif", 1:12),
  dim_name = "month",
  start = 1,
  increment = 1,
  output = "monthly_precip.vrt"
)

# Pressure levels at 50 hPa intervals
vrt_regular_stack(
  files = level_files,
  dim_name = "pressure",
  start = 1000,
  increment = -50,
  unit = "hPa",
  output = "pressure_levels.vrt"
)

## End(Not run)

---

Stack Time Series Rasters into Multidimensional VRT

Description

A convenience wrapper for gdal_mdim_stack() tailored for time series data. Extracts dates from filenames and encodes them as numeric offsets from an origin.

Usage

vrt_time_stack(
  files,
  pattern = "(\\d{8})",
  format = "%Y%m%d",
  origin = "1970-01-01",
  unit = "days",
  array_name = NULL,
  output = NULL
)

Arguments

files	Character vector of file paths
pattern	Regex pattern with one capturing group for the date portion. Default: "(\\d{8})" matches 8-digit dates like 20230115.
format	strptime format for parsing. Default: "%Y%m%d"
origin	Time origin as "YYYY-MM-DD" string. Default: "1970-01-01"
unit	Time unit: "days", "hours", "minutes", "seconds". Default: "days"
array_name	Name for the data array
output	Output VRT path

Value

Path to VRT file (invisibly), or XML string if output is NULL.

Examples

## Not run: 
# Daily temperature files: temp_20230101.tif, temp_20230102.tif, ...
vrt_time_stack(
  files = list.files(pattern = "temp_\\d{8}\\.tif$"),
  output = "temperature_timeseries.vrt"
)

# Custom pattern for ISO dates: data_2023-01-15.tif
vrt_time_stack(
  files = my_files,
  pattern = "(\\d{4}-\\d{2}-\\d{2})",
  format = "%Y-%m-%d",
  output = "timeseries.vrt"
)

## End(Not run)

---

Stack rasters into a multidimensional VRT (two paths, one entry point)

Description

vrtstack() builds a GDAL multidimensional VRT that concatenates a set of sources along one axis. Two construction paths share a single concat-axis resolver:

Usage

vrtstack(
  files,
  concat = NULL,
  template = NULL,
  concat_dim = "time",
  n_template = 10,
  parse_format = NULL,
  origin = NULL,
  unit = "days",
  metadata_from = "first",
  output = NULL,
  ...
)

Arguments

files	Character vector of source paths/URIs (full set).
concat	Concat-axis spec, resolved by resolve_concat(): a numeric/ character vector (length == files), a single regex with one capture group to pull from basenames, a function of basename, or NULL to use concat_start/concat_increment (from-scratch) or the template's own values (template-extend, when concat is NULL).
template	One of: NULL (from-scratch); a path to a .vrt; an XML string; or TRUE to auto-generate the recipe from head(files, n_template) via ⁠gdal mdim mosaic⁠.
concat_dim	Name of the concatenation dimension. Default "time".
n_template	Files used to build the auto-generated recipe when template = TRUE. Default 10.
parse_format, origin, unit	Passed to resolve_concat() for date parsing.
metadata_from	Which source donates the coordinate units and root-group provenance restored onto an auto-generated recipe (template = TRUE): one of "first" (default), "last", or an explicit source path. Ignored when template is a supplied VRT or XML string, which already carries its own metadata.
output	Output VRT path, or NULL to return the XML string.
...	Passed to the from-scratch builder (dim_type, dim_unit, array_name, etc.) and ignored on the template path.

Details

• from-scratch (template = NULL) - the original GeoTIFF/COG path. Each 2D source contributes one slab of a single new 3D array. Sources use the classic ⁠<SourceBand>⁠ + ⁠<SourceTranspose>⁠ form.

• template-extend (⁠template = <vrt>⁠ or template = TRUE) - for mdim-recognised sources (NetCDF/Zarr/HDF5) that may carry several arrays (OISST: anom/err/ice/sst). A small "first-n" mdim VRT is taken as the recipe, and the full file set is fleshed out by re-parameterising only the concat dimension and the per-array ⁠<Source>⁠ lists. Everything else (coordinate arrays, slab geometry, scale/offset/unit/attributes) is carried verbatim from the template.

The template is the GDAL-native analogue of the "recipe-not-payload" principle: GDAL's own ⁠gdal mdim mosaic⁠ writes the recipe for the first n files; we only template the axis that varies with the file set.

Value

Path to the VRT (invisibly) or the XML string if output is NULL.

Examples

#vrt <- vrtstack(<files>, concat = "(\d{8})", parse_format = "%Y%m%d",
#origin = "1978-01-01", unit = "days",   # must match the template
#template = TRUE, concat_dim = "time")

---