Open modes

Classic rasters (bands as a dimension or as variables) and multidimensional arrays.

gdalxarray opens a source in one of two GDAL data models, chosen by multidim=:

	multidim=False (classic)	multidim=True (default)
GDAL API	gdal.Open, raster bands	gdal.OpenEx(..., OF_MULTIDIM_RASTER), groups and arrays
Dimensions	(band, y, x)	whatever the source names: time, depth, lat, lon, …
Georeferencing	geotransform + CRS, as a RasterIndex and CRSIndex	coordinate arrays from the source
Natural for	GeoTIFF/COG, JPEG2000, classic VRT, warped VRT	netCDF, HDF5, Zarr, GRIB, multidim VRT, Icechunk

Classic, bands as a dimension

The default for multidim=False: one band_data variable with dimensions (band, y, x), the layout rioxarray uses. The example is a cloud-free Sentinel-2 summer scene over Hobart, from kunanyi / Mount Wellington to the Derwent estuary: four bands (blue, green, red, near infrared) on the scene’s own UTM zone 55S grid, averaged to 20 m.

Data: Copernicus Sentinel-2 Level-2A surface reflectance; licence: Copernicus Sentinel data terms (free and open); attribution: contains modified Copernicus Sentinel data; day 2024-12-04. Built by scripts/fixtures/make_sentinel2.py; details in docs/fixtures/sentinel2/provenance.json.

src = data("sentinel2/" + data_json("sentinel2/provenance.json")["file"])
ds = xr.open_dataset(src, engine="gdalxarray", multidim=False)
ds

<xarray.Dataset> Size: 15MB
Dimensions:    (band: 4, y: 625, x: 750)
Coordinates:
  * band       (band) int64 32B 1 2 3 4
  * y          (y) float64 5kB 5.255e+06 5.255e+06 ... 5.243e+06 5.243e+06
  * x          (x) float64 6kB 5.15e+05 5.15e+05 5.15e+05 ... 5.3e+05 5.3e+05
  * crs        int64 8B 0
Data variables:
    band_data  (band, y, x) float64 15MB ...
Indexes:
  ┌ x        RasterIndex (crs=EPSG:32755)
  └ y
    crs      CRSIndex (crs=EPSG:32755)

xarray.Dataset

• Dimensions:
  • band: 4
  • y: 625
  • x: 750
• Coordinates: (4)
  • band
    (band)
    int64
    1 2 3 4

    array([1, 2, 3, 4])

  • y
    (y)
    float64
    5.255e+06 5.255e+06 ... 5.243e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    [625 values with dtype=float64]

  • x
    (x)
    float64
    5.15e+05 5.15e+05 ... 5.3e+05

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    [750 values with dtype=float64]

  • crs
    ()
    int64
    0

    array(0)

• Data variables: (1)
  • band_data
    (band, y, x)
    float64
    ...

    descriptions :

    ['blue', 'green', 'red', 'nir']

    [1875000 values with dtype=float64]

• Indexes: (2)
  • x
    y
    RasterIndex (crs=EPSG:32755)

    RasterIndex(crs=EPSG:32755)
        AxisAffineTransformIndex(AxisAffineTransform(a=20, b=0, c=5.15e+05, d=0, e=-20, f=5.255e+06, axis=X, dim='x'))
        AxisAffineTransformIndex(AxisAffineTransform(a=20, b=0, c=5.15e+05, d=0, e=-20, f=5.255e+06, axis=Y, dim='y'))

  • crs
    CRSIndex (crs=EPSG:32755)

    CRSIndex
    <Projected CRS: EPSG:32755>
    Name: WGS 84 / UTM zone 55S
    Axis Info [cartesian]:
    - E[east]: Easting (metre)
    - N[north]: Northing (metre)
    Area of Use:
    - undefined
    Coordinate Operation:
    - name: UTM zone 55S
    - method: Transverse Mercator
    Datum: World Geodetic System 1984
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

The cell coordinates come from the geotransform, and the CRS rides on the crs coordinate’s index:

ds.xindexes["crs"]

CRSIndex
<Projected CRS: EPSG:32755>
Name: WGS 84 / UTM zone 55S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: UTM zone 55S
- method: Transverse Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The bands are stored as integers with a scale and offset (Sentinel-2 L2A digital numbers); with the default mask_and_scale=True they arrive as surface reflectance, with nodata masked:

ds["band_data"].encoding.get("scale_factor"), ds["band_data"].encoding.get("add_offset"), \
    float(ds["band_data"].mean())

(0.0001, -0.1, 0.08045086261333333)

Ordinary xarray operations apply. A true-colour composite is a selection and reordering along band:

rgb = ds["band_data"].sel(band=[3, 2, 1])               # red, green, blue
rgb = (rgb / rgb.quantile(0.98)).clip(0, 1).transpose("y", "x", "band")
fig, ax = plt.subplots(figsize=(8, 6.7))
rgb.plot.imshow(ax=ax, rgb="band")
ax.set_aspect("equal")
ax.set_title("true colour, " + data_json("sentinel2/provenance.json")["day"])
ax.set_xlabel("easting (m)"); ax.set_ylabel("northing (m)")
plt.show()

Sentinel-2 true colour (red, green, blue reflectance, stretched) over Hobart, kunanyi / Mount Wellington and the Derwent estuary.

Classic, bands as variables

band_as_dim=False makes each band its own variable, named from the band description (or band_N when there is none). This suits rasters whose bands are different quantities, and makes band arithmetic read like the quantities it combines:

dsv = xr.open_dataset(src, engine="gdalxarray", multidim=False, band_as_dim=False)
list(dsv.data_vars), dsv["nir"].dims

(['blue', 'green', 'red', 'nir'], ('y', 'x'))

ndvi = (dsv["nir"] - dsv["red"]) / (dsv["nir"] + dsv["red"])
fig, ax = plt.subplots(figsize=(8, 6.7))
ndvi.plot(ax=ax, cmap="RdYlGn", vmin=-0.5, vmax=0.9, cbar_kwargs={"label": "NDVI"})
ax.set_aspect("equal"); ax.set_title("NDVI")
ax.set_xlabel("easting (m)"); ax.set_ylabel("northing (m)")
plt.show()

Normalised difference vegetation index, (nir - red) / (nir + red), computed from the named band variables.

Multidimensional

For sources with their own dimensions. The fixture below is a ROMS-style ocean model history file written as Zarr v3: a vertical s_rho axis, two index dimensions, 2D longitude and latitude arrays, and the scalar parameters of the vertical transform.

roms = xr.open_dataset(fixture("index/roms-fan.zarr"), engine="gdalxarray")
roms

<xarray.Dataset> Size: 80kB
Dimensions:      (eta_rho: 46, xi_rho: 24, s_rho: 8)
Coordinates:
  * eta_rho      (eta_rho) int64 368B 0 1 2 3 4 5 6 7 ... 39 40 41 42 43 44 45
  * xi_rho       (xi_rho) int64 192B 0 1 2 3 4 5 6 7 ... 16 17 18 19 20 21 22 23
    lon_rho      (eta_rho, xi_rho) float64 9kB ...
    lat_rho      (eta_rho, xi_rho) float64 9kB ...
  * s_rho        (s_rho) float64 64B -0.9375 -0.8125 -0.6875 ... -0.1875 -0.0625
Data variables:
    Vstretching  float64 8B ...
    temp         (s_rho, eta_rho, xi_rho) float32 35kB ...
    h            (eta_rho, xi_rho) float64 9kB ...
    Vtransform   float64 8B ...
    hc           float64 8B ...
    mask_rho     (eta_rho, xi_rho) float64 9kB ...
    Cs_r         (s_rho) float64 64B ...
    theta_b      float64 8B ...
    theta_s      float64 8B ...
    zeta         (eta_rho, xi_rho) float64 9kB ...
Attributes:
    Conventions:  CF-1.8
    title:        wxylib fixtures I1 roms-fan
    type:         ROMS/TOMS history file (synthetic)

xarray.Dataset

• Dimensions:
  • eta_rho: 46
  • xi_rho: 24
  • s_rho: 8
• Coordinates: (5)
  • eta_rho
    (eta_rho)
    int64
    0 1 2 3 4 5 6 ... 40 41 42 43 44 45

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39, 40, 41, 42, 43, 44, 45])

  • xi_rho
    (xi_rho)
    int64
    0 1 2 3 4 5 6 ... 18 19 20 21 22 23

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23])

  • lon_rho
    (eta_rho, xi_rho)
    float64
    ...

    long_name :

    longitude of RHO-points

    standard_name :

    longitude

    units :

    degree_east

    [1104 values with dtype=float64]

  • lat_rho
    (eta_rho, xi_rho)
    float64
    ...

    long_name :

    latitude of RHO-points

    standard_name :

    latitude

    units :

    degree_north

    [1104 values with dtype=float64]

  • s_rho
    (s_rho)
    float64
    -0.9375 -0.8125 ... -0.1875 -0.0625

    formula_terms :

    s: s_rho C: Cs_r eta: zeta depth: h depth_c: hc

    long_name :

    S-coordinate at RHO-points

    positive :

    up

    standard_name :

    ocean_s_coordinate_g2

    valid_max :

    0.0

    valid_min :

    -1.0

    array([-0.9375, -0.8125, -0.6875, -0.5625, -0.4375, -0.3125, -0.1875, -0.0625])

• Data variables: (10)
  • Vstretching
    ()
    float64
    ...

    long_name :

    vertical terrain-following stretching function

    [1 values with dtype=float64]

  • temp
    (s_rho, eta_rho, xi_rho)
    float32
    ...

    long_name :

    potential temperature

    units :

    Celsius

    [8832 values with dtype=float32]

  • h
    (eta_rho, xi_rho)
    float64
    ...

    long_name :

    bathymetry at RHO-points

    units :

    meter

    [1104 values with dtype=float64]

  • Vtransform
    ()
    float64
    ...

    long_name :

    vertical terrain-following transformation equation

    [1 values with dtype=float64]

  • hc
    ()
    float64
    ...

    long_name :

    S-coordinate parameter, critical depth

    units :

    meter

    [1 values with dtype=float64]

  • mask_rho
    (eta_rho, xi_rho)
    float64
    ...

    flag_meanings :

    land water

    flag_values :

    (0.0, 1.0)

    long_name :

    mask on RHO-points

    [1104 values with dtype=float64]

  • Cs_r
    (s_rho)
    float64
    ...

    long_name :

    S-coordinate stretching curves at RHO-points

    valid_max :

    0.0

    valid_min :

    -1.0

    [8 values with dtype=float64]

  • theta_b
    ()
    float64
    ...

    long_name :

    S-coordinate bottom control parameter

    [1 values with dtype=float64]

  • theta_s
    ()
    float64
    ...

    long_name :

    S-coordinate surface control parameter

    [1 values with dtype=float64]

  • zeta
    (eta_rho, xi_rho)
    float64
    ...

    long_name :

    free-surface

    units :

    meter

    [1104 values with dtype=float64]

• Attributes: (3)

  Conventions :

  CF-1.8

  title :

  wxylib fixtures I1 roms-fan

  type :

  ROMS/TOMS history file (synthetic)

Dimensions, coordinates and attributes come from the store, so selection reads like the model’s own description:

fig, axs = plt.subplots(1, 2, figsize=(7, 4.5), sharey=True, layout="constrained")
for ax, k, name in [(axs[0], -1, "surface"), (axs[1], 0, "bottom")]:
    roms["temp"].isel(s_rho=k).plot(ax=ax, cmap="RdYlBu_r", vmin=float(roms.temp.min()),
                                     vmax=float(roms.temp.max()), add_colorbar=(k == 0))
    ax.set_title(f"{name} (s_rho = {float(roms.s_rho[k]):.2f})")
plt.show()

Temperature on the surface level and the deepest level of the fan grid, by index.

Placing these cells on the earth uses the 2D lon_rho and lat_rho arrays; see Curvilinear grids.

Checked against GDAL

The three Zarr fixtures in wxylib/fixtures hold the same pixels as three of the COGs, and their manifest lists windows whose values GDAL read with osgeo.gdal. gdalxarray must give the same numbers through xarray indexing, in both modes:

man = fixture_json("zarr/manifest.json")
rows = []
for case in man["cases"]:
    path = fixture("zarr/" + case["store"])
    md = xr.open_dataset(path, engine="gdalxarray", mask_and_scale=False)[case["array"]]
    cl = xr.open_dataset(f'ZARR:"{path}":/{case["array"]}', engine="gdalxarray",
                         multidim=False, mask_and_scale=False)["band_data"]
    for w in case["windows"]:
        win = dict(y=slice(w["row"], w["row"] + w["height"]),
                   x=slice(w["col"], w["col"] + w["width"]))
        expect = fixture_array("zarr/" + w["file"], case["dtype"])
        a = md.isel(**win).values
        b = cl.isel(**win).values
        ok_md = np.array_equal(a.ravel(), expect)
        ok_cl = np.array_equal(b.ravel(), expect)
        rows.append((case["store"], w["file"].rsplit("/", 1)[-1], a.shape, ok_md, ok_cl))
        assert ok_md and ok_cl, (case["store"], w["file"])

import pandas as pd
pd.DataFrame(rows, columns=["store", "oracle window", "shape", "multidim", "classic"])

	store	oracle window	shape	multidim	classic
0	world-int16-spatial.zarr	world-int16-deflate-w0.bin	(40, 50)	True	True
1	world-int16-spatial.zarr	world-int16-deflate-w1.bin	(77, 130)	True	True
2	world-int16-spatial.zarr	world-int16-deflate-w2.bin	(45, 70)	True	True
3	polar-f32-cf.zarr	polar-f32-zstd-bigtiff-w0.bin	(40, 50)	True	True
4	polar-f32-cf.zarr	polar-f32-zstd-bigtiff-w1.bin	(77, 130)	True	True
5	polar-f32-cf.zarr	polar-f32-zstd-bigtiff-w2.bin	(45, 70)	True	True
6	utm-u16-3band-node.zarr	utm-u16-lzw-pred2-3band-w0.bin	(3, 40, 50)	True	True
7	utm-u16-3band-node.zarr	utm-u16-lzw-pred2-3band-w1.bin	(3, 77, 130)	True	True
8	utm-u16-3band-node.zarr	utm-u16-lzw-pred2-3band-w2.bin	(3, 45, 70)	True	True

Rendered with GDAL 3.13.3, xarray 2026.4.1.dev32+g71cfb13ac, numpy 2.5.4, gdalxarray 0.5.1, rasterix 0.2.4.