Get started

This page takes a set of daily NetCDF files, opens them with xarray, writes a VRT describing them as one cube, and reads it back through GDAL.

Some files

Four daily files of sea surface temperature, each (time, lat, lon) with one time step, stored as packed int16 (the usual layout for gridded daily products). The helper writes them to a temporary folder.

import xarray as xr
from _fixtures import daily_netcdf

paths = daily_netcdf(n=4)
[p.name for p in paths]

['sst_20260101.nc', 'sst_20260102.nc', 'sst_20260103.nc', 'sst_20260104.nc']

Open them as usual

ds = xr.open_mfdataset(paths, combine="nested", concat_dim="time")
ds

<xarray.Dataset> Size: 472B
Dimensions:  (time: 4, lat: 3, lon: 4)
Coordinates:
  * time     (time) datetime64[ns] 32B 2026-01-01T12:00:00 ... 2026-01-04T12:...
  * lat      (lat) float64 24B -44.5 -43.5 -42.5
  * lon      (lon) float64 32B 145.5 146.5 147.5 148.5
Data variables:
    sst      (time, lat, lon) float64 384B dask.array<chunksize=(1, 3, 4), meta=np.ndarray>
Attributes:
    title:    documentation fixture

xarray has decoded the packing (the values are floats) and the times (datetime64), and remembers how they were stored in .encoding.

Write the VRT

write_mdim_vrt needs the Dataset, the source files in order, an output path, and how the files combine: here each file contributes its slice of an existing time dimension, which is composition="concat".

from xvrt import write_mdim_vrt

out = write_mdim_vrt(
    ds, paths, paths[0].parent / "sst.vrt",
    composition="concat", concat_dim="time",
)
print(out.read_text())

<?xml version='1.0' encoding='utf-8'?>
<VRTDataset>
  <Group name="/">
    <Dimension name="time" size="4" indexingVariable="time" type="TEMPORAL" />
    <Dimension name="lat" size="3" indexingVariable="lat" type="HORIZONTAL_Y" direction="NORTH" />
    <Dimension name="lon" size="4" indexingVariable="lon" type="HORIZONTAL_X" direction="EAST" />
    <Attribute name="title">
      <DataType>String</DataType>
      <Value>documentation fixture</Value>
    </Attribute>
    <Array name="time">
      <DataType>Float64</DataType>
      <DimensionRef ref="time" />
      <Unit>days since 2026-01-01 12:00:00</Unit>
      <RegularlySpacedValues start="0.0" increment="1.0" />
    </Array>
    <Array name="lat">
      <DataType>Float64</DataType>
      <DimensionRef ref="lat" />
      <Unit>degrees_north</Unit>
      <RegularlySpacedValues start="-44.5" increment="1.0" />
    </Array>
    <Array name="lon">
      <DataType>Float64</DataType>
      <DimensionRef ref="lon" />
      <Unit>degrees_east</Unit>
      <RegularlySpacedValues start="145.5" increment="1.0" />
    </Array>
    <Array name="sst">
      <DataType>Int16</DataType>
      <DimensionRef ref="time" />
      <DimensionRef ref="lat" />
      <DimensionRef ref="lon" />
      <BlockSize>1,3,4</BlockSize>
      <Unit>degC</Unit>
      <NoDataValue>-32768.0</NoDataValue>
      <Offset>0.0</Offset>
      <Scale>0.01</Scale>
      <Source>
        <SourceFilename relativeToVRT="1">sst_20260101.nc</SourceFilename>
        <SourceArray>/sst</SourceArray>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="0,0,0" />
      </Source>
      <Source>
        <SourceFilename relativeToVRT="1">sst_20260102.nc</SourceFilename>
        <SourceArray>/sst</SourceArray>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="1,0,0" />
      </Source>
      <Source>
        <SourceFilename relativeToVRT="1">sst_20260103.nc</SourceFilename>
        <SourceArray>/sst</SourceArray>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="2,0,0" />
      </Source>
      <Source>
        <SourceFilename relativeToVRT="1">sst_20260104.nc</SourceFilename>
        <SourceArray>/sst</SourceArray>
        <SourceSlab offset="0,0,0" count="1,3,4" step="1,1,1" />
        <DestSlab offset="3,0,0" />
      </Source>
      <Attribute name="long_name">
        <DataType>String</DataType>
        <Value>daily sst</Value>
      </Attribute>
    </Array>
  </Group>
</VRTDataset>

A few things to notice:

The data array is declared Int16 with <Scale>, <Offset> and <NoDataValue>: the VRT reads the stored bytes, and GDAL applies the decoding, so nothing is lost or double-applied. See Packed data.
Regularly spaced coordinates are written as <RegularlySpacedValues> (start and step), not as lists.
Time is written back in its original CF units, origin time of day included.
Each <Source> states its extent (<SourceSlab>) and where it lands (<DestSlab>), so GDAL opens only the files a read needs.
No CRS was found on this Dataset, so there is no <SRS> (xvrt warns). Pass crs= or see CRS.

Read it back

from osgeo import gdal

gdal.UseExceptions()
vrt = gdal.OpenEx(str(out), gdal.OF_MULTIDIM_RASTER)
sst = vrt.GetRootGroup().OpenMDArray("sst")
print([d.GetName() for d in sst.GetDimensions()], [d.GetSize() for d in sst.GetDimensions()])
print("scale", sst.GetScale(), "nodata", sst.GetNoDataValueAsDouble(), "unit", sst.GetUnit())

['time', 'lat', 'lon'] [4, 3, 4]
scale 0.01 nodata -32768.0 unit degC

The unscaled values match what xarray shows:

import numpy as np

from_vrt = sst.GetUnscaled().ReadAsArray()
np.allclose(from_vrt, ds.sst.values)

True

The same file opens in any tool built on GDAL’s multidimensional API, for example gdal mdim info sst.vrt on the command line, or an xarray backend for GDAL:

xr.open_dataset(out, engine="gdalxarray")