Raster sources: GeoTIFF and COG

Classic raster files (GeoTIFF, Cloud Optimized GeoTIFF) have bands, not named arrays. xvrt stacks them with <SourceBand> and a <SourceTranspose> that maps the destination’s (time, y, x) onto the band’s (y, x), with -1 marking the new length-1 axis.

import numpy as np
import xarray as xr
from osgeo import gdal
from xvrt import write_mdim_vrt
from _fixtures import geotiffs, WKT_4326

gdal.UseExceptions()
tifs = geotiffs(n=3)

# Describe the cube: here directly from what we know about the files.
raw = np.stack([gdal.Open(str(p)).ReadAsArray() for p in tifs])
da = xr.DataArray(raw, dims=("time", "y", "x"), name="analysed_sst")
da.encoding.update({"dtype": np.dtype("int16"), "_FillValue": -32768,
                    "scale_factor": 0.001, "add_offset": 25.0})
ds = da.to_dataset().assign_coords(
    time=("time", np.array(["2026-01-01", "2026-01-02", "2026-01-03"], "datetime64[ns]")),
    y=-42.5 - np.arange(3.0), x=145.5 + np.arange(4.0))

out = write_mdim_vrt(ds, tifs, tifs[0].parent / "cogs.vrt",
                     composition="stack", concat_dim="time", crs=WKT_4326)
print("\n".join(out.read_text().splitlines()[-14:]))

        <SourceTranspose>-1,0,1</SourceTranspose>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="1,0,0" />
      </Source>
      <Source>
        <SourceFilename relativeToVRT="1">20260103_analysed_sst.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SourceTranspose>-1,0,1</SourceTranspose>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="2,0,0" />
      </Source>
    </Array>
  </Group>
</VRTDataset>

Paths ending .tif or .tiff default to band 1; dict-form sources take "band": N. Band sources only make sense in a stack (a 2D band has no dimension to concatenate along).

Reading only what is needed

Every source carries <SourceSlab offset count step>, stating its extent in the destination’s terms. Without it, GDAL opens every source, in order, to learn its shape before any read: for thousands of remote COGs that is a HEAD request and a header read per file, on every access. With it, a read of one day or one pixel series opens only the files it touches.

Classic VRT or multidimensional VRT?

For maps of one date, overviews matter, and a classic raster VRT (or the COG itself) uses them. A multidimensional VRT does not use overviews; its strength is the time (or depth) axis: pixel and window series across many files, with the metadata of every array carried through. The two are complementary.