Foreign GeoZarr stores — reading what other tools write¶

Zarr stores written by rioxarray.to_zarr, odc-geo, or GDAL's Zarr driver follow the same CF / GeoZarr convention pyramids does — but they're allowed to vary in two annoying ways:

1. The primary data array is named after the variable, e.g. "elevation" or "sst", not "data".
2. The CRS variable is the only place the geo info lives — there's no GeoTransform attr on the data array, just 1-D x / y pixel-centre coordinates plus a spatial_ref grid-mapping variable.

Dataset.from_zarr handles both. It auto-detects the primary array via the CF grid_mapping link (with a data_name= override for ambiguous stores) and derives the transform from the x/y coords when no explicit GeoTransform is present.

This notebook builds two foreign-style stores by hand (no pyramids writer involved) so you can see exactly which on-disk fields pyramids needs to read them back.

Setup¶

In [1]:

import tempfile
from pathlib import Path

import numpy as np
import zarr
from pyproj import CRS

from pyramids.dataset import Dataset
from pyramids.plot import ColorBar

workspace = Path(tempfile.mkdtemp(prefix="pyramids-zarr-foreign-"))
print("workspace:", workspace)

ROWS, COLS = 16, 24
CELL = 30.0  # 30 m pixels
EPSG = 32633  # UTM zone 33N
ULX, ULY = 500000.0, 4_500_000.0  # arbitrary anchor in projected coords

elevation = np.arange(ROWS * COLS, dtype=np.float32).reshape(1, ROWS, COLS) + 100.0
# When the source has a single band, Dataset.read_array() returns a 2-D array,
# so we compare against the squeezed view below.
elevation_2d = elevation[0]
print("on-disk shape (band, y, x):", elevation.shape, "dtype:", elevation.dtype)

import tempfile from pathlib import Path import numpy as np import zarr from pyproj import CRS from pyramids.dataset import Dataset from pyramids.plot import ColorBar workspace = Path(tempfile.mkdtemp(prefix="pyramids-zarr-foreign-")) print("workspace:", workspace) ROWS, COLS = 16, 24 CELL = 30.0 # 30 m pixels EPSG = 32633 # UTM zone 33N ULX, ULY = 500000.0, 4_500_000.0 # arbitrary anchor in projected coords elevation = np.arange(ROWS * COLS, dtype=np.float32).reshape(1, ROWS, COLS) + 100.0 # When the source has a single band, Dataset.read_array() returns a 2-D array, # so we compare against the squeezed view below. elevation_2d = elevation[0] print("on-disk shape (band, y, x):", elevation.shape, "dtype:", elevation.dtype)

workspace: /tmp/pyramids-zarr-foreign-zzvxz66d
on-disk shape (band, y, x): (1, 16, 24) dtype: float32

1. Build a foreign GeoZarr store by hand¶

Lay out the array as something like odc-geo or rioxarray would produce it:

• a named data array elevation with _ARRAY_DIMENSIONS and grid_mapping="spatial_ref"
• a scalar spatial_ref variable carrying the WKT + GeoTransform
• 1-D x / y pixel-centre coordinates

No pyramids_zarr_version attribute, no data array.

In [2]:

store_a = workspace / "named-array.zarr"

root = zarr.open_group(store_a, mode="w")

# x and y at pixel centres
x_coords = ULX + (np.arange(COLS) + 0.5) * CELL
y_coords = ULY - (np.arange(ROWS) + 0.5) * CELL
root.create_array("x", shape=x_coords.shape, dtype=x_coords.dtype).attrs.update(
    {"_ARRAY_DIMENSIONS": ["x"]}
)
root.create_array("y", shape=y_coords.shape, dtype=y_coords.dtype).attrs.update(
    {"_ARRAY_DIMENSIONS": ["y"]}
)
root["x"][:] = x_coords
root["y"][:] = y_coords

# grid-mapping variable: scalar, WKT + GeoTransform on attrs
wkt = CRS.from_epsg(EPSG).to_wkt()
gt = [ULX, CELL, 0.0, ULY, 0.0, -CELL]
root.create_array("spatial_ref", shape=(), dtype="int32")
root["spatial_ref"].attrs.update(
    {
        "crs_wkt": wkt,
        "GeoTransform": " ".join(str(v) for v in gt),
    }
)

# data array — named after the variable, not "data"
root.create_array("elevation", shape=elevation.shape, dtype=elevation.dtype)
root["elevation"][:] = elevation
root["elevation"].attrs.update(
    {
        "_ARRAY_DIMENSIONS": ["band", "y", "x"],
        "grid_mapping": "spatial_ref",
    }
)

print("arrays:", sorted(root.array_keys()))

store_a = workspace / "named-array.zarr" root = zarr.open_group(store_a, mode="w") # x and y at pixel centres x_coords = ULX + (np.arange(COLS) + 0.5) * CELL y_coords = ULY - (np.arange(ROWS) + 0.5) * CELL root.create_array("x", shape=x_coords.shape, dtype=x_coords.dtype).attrs.update( {"_ARRAY_DIMENSIONS": ["x"]} ) root.create_array("y", shape=y_coords.shape, dtype=y_coords.dtype).attrs.update( {"_ARRAY_DIMENSIONS": ["y"]} ) root["x"][:] = x_coords root["y"][:] = y_coords # grid-mapping variable: scalar, WKT + GeoTransform on attrs wkt = CRS.from_epsg(EPSG).to_wkt() gt = [ULX, CELL, 0.0, ULY, 0.0, -CELL] root.create_array("spatial_ref", shape=(), dtype="int32") root["spatial_ref"].attrs.update( { "crs_wkt": wkt, "GeoTransform": " ".join(str(v) for v in gt), } ) # data array — named after the variable, not "data" root.create_array("elevation", shape=elevation.shape, dtype=elevation.dtype) root["elevation"][:] = elevation root["elevation"].attrs.update( { "_ARRAY_DIMENSIONS": ["band", "y", "x"], "grid_mapping": "spatial_ref", } ) print("arrays:", sorted(root.array_keys()))

arrays: ['elevation', 'spatial_ref', 'x', 'y']

2. Read it back via Dataset.from_zarr¶

Two ways to point pyramids at the right array:

• Auto-detect: pyramids picks the array whose grid_mapping attr names the CRS variable. Works whenever exactly one data array sets grid_mapping.
• Explicit: pass data_name="elevation". Works regardless of attrs.

In [3]:

auto = Dataset.from_zarr(store_a)
explicit = Dataset.from_zarr(store_a, data_name="elevation")

for label, d in (("auto-detected", auto), ("data_name=elevation", explicit)):
    print(f"{label:<22} shape={d.shape}  epsg={d.epsg}  cell_size={d.cell_size}")
    np.testing.assert_array_equal(d.read_array(), elevation_2d)
print("both reads match the original array")

auto = Dataset.from_zarr(store_a) explicit = Dataset.from_zarr(store_a, data_name="elevation") for label, d in (("auto-detected", auto), ("data_name=elevation", explicit)): print(f"{label:<22} shape={d.shape} epsg={d.epsg} cell_size={d.cell_size}") np.testing.assert_array_equal(d.read_array(), elevation_2d) print("both reads match the original array")

auto-detected          shape=(1, 16, 24)  epsg=32633  cell_size=30.0
data_name=elevation    shape=(1, 16, 24)  epsg=32633  cell_size=30.0
both reads match the original array

In [4]:

# --- Visualise the foreign-store elevation ---
# We hand-built the store the way odc-geo / rioxarray would, with the data array
# named `elevation` (not `data`). Rendering `auto` (the auto-detected read) via
# `Dataset.plot()` confirms two things at once: (1) the auto-detect picked the right
# array, and (2) row/col orientation came back the right way up — top-left of the
# image is the high-y / low-x corner that the y-coord (descending) defined.
auto.plot(
    band=0,
    figsize=(7, 4),
    title=f"elevation  (EPSG:{auto.epsg}, {auto.cell_size} m pixels)",
    colorbar=ColorBar(label="elevation (m)"),
    cmap="terrain",
)

# --- Visualise the foreign-store elevation --- # We hand-built the store the way odc-geo / rioxarray would, with the data array # named `elevation` (not `data`). Rendering `auto` (the auto-detected read) via # `Dataset.plot()` confirms two things at once: (1) the auto-detect picked the right # array, and (2) row/col orientation came back the right way up — top-left of the # image is the high-y / low-x corner that the y-coord (descending) defined. auto.plot( band=0, figsize=(7, 4), title=f"elevation (EPSG:{auto.epsg}, {auto.cell_size} m pixels)", colorbar=ColorBar(label="elevation (m)"), cmap="terrain", )

Out[4]:

<cleopatra.glyphs.gridded.array_glyph.ArrayGlyph at 0x7f0f141f3cb0>

3. Derive the transform from x / y when GeoTransform is absent¶

Plenty of foreign stores don't bother emitting an explicit GeoTransform. As long as the x and y pixel-centre coordinates are there, pyramids reconstructs the transform from them.

In [5]:

store_b = workspace / "no-geotransform.zarr"

root = zarr.open_group(store_b, mode="w")

root.create_array("x", shape=x_coords.shape, dtype=x_coords.dtype).attrs.update(
    {"_ARRAY_DIMENSIONS": ["x"]}
)
root.create_array("y", shape=y_coords.shape, dtype=y_coords.dtype).attrs.update(
    {"_ARRAY_DIMENSIONS": ["y"]}
)
root["x"][:] = x_coords
root["y"][:] = y_coords

# CRS only — NO GeoTransform attr
root.create_array("spatial_ref", shape=(), dtype="int32")
root["spatial_ref"].attrs.update({"crs_wkt": wkt})

root.create_array("elevation", shape=elevation.shape, dtype=elevation.dtype)
root["elevation"][:] = elevation
root["elevation"].attrs.update(
    {
        "_ARRAY_DIMENSIONS": ["band", "y", "x"],
        "grid_mapping": "spatial_ref",
    }
)

rt = Dataset.from_zarr(store_b)
print("epsg      :", rt.epsg)
print("cell_size :", rt.cell_size, " (recovered from x/y spacing)")
print("top-left  :", rt.top_left_corner)
np.testing.assert_array_equal(rt.read_array(), elevation_2d)
print("transform derivation OK")

store_b = workspace / "no-geotransform.zarr" root = zarr.open_group(store_b, mode="w") root.create_array("x", shape=x_coords.shape, dtype=x_coords.dtype).attrs.update( {"_ARRAY_DIMENSIONS": ["x"]} ) root.create_array("y", shape=y_coords.shape, dtype=y_coords.dtype).attrs.update( {"_ARRAY_DIMENSIONS": ["y"]} ) root["x"][:] = x_coords root["y"][:] = y_coords # CRS only — NO GeoTransform attr root.create_array("spatial_ref", shape=(), dtype="int32") root["spatial_ref"].attrs.update({"crs_wkt": wkt}) root.create_array("elevation", shape=elevation.shape, dtype=elevation.dtype) root["elevation"][:] = elevation root["elevation"].attrs.update( { "_ARRAY_DIMENSIONS": ["band", "y", "x"], "grid_mapping": "spatial_ref", } ) rt = Dataset.from_zarr(store_b) print("epsg :", rt.epsg) print("cell_size :", rt.cell_size, " (recovered from x/y spacing)") print("top-left :", rt.top_left_corner) np.testing.assert_array_equal(rt.read_array(), elevation_2d) print("transform derivation OK")

epsg      : 32633
cell_size : 30.0  (recovered from x/y spacing)
top-left  : (500000.0, 4500000.0)
transform derivation OK

Where to next¶

• Single-raster zarr basics: zarr-basics.ipynb.
• Cube workflow with append / region writes: zarr-cube.ipynb.
• Multiscale pyramids for fast previews: zarr-pyramid-preview.ipynb.
• Full API + on-disk layout: Zarr reference.