Remove ndpyramid spec & add ability to pass None to sharding encoding

README.md

@@ -1,5 +1,7 @@
 # topozarr - lightweight multiscale zarr pyramids
 
+Python library to create multiscale zarr pyramids for usage with [zarr-layer](https://zarr-layer.demo.carbonplan.org/).
+Attempts to follow the WIP [zarr-multiscales spec](https://github.com/zarr-conventions/multiscales).
 
 **Warning: experimental**
 
@@ -31,25 +33,38 @@ pyramid = create_pyramid(
     ds, 
     levels=2, 
     x_dim="lon", 
-    y_dim="lat", 
-    spec="ndpyramid" # or "zarr-multiscales"
-)
+    y_dim="lat")
 print(pyramid.encoding)
 print(pyramid.dt)
 ```
 
 ```python
 # Optional: Write to Zarr
-!pip install obstore zarr
+# !pip install obstore zarr
 from obstore.store import from_url
 from zarr.storage import ObjectStore
 
 
 store = from_url(url = "<add_your_bucket_url>", region="<add_your_region>")
 zstore = ObjectStore(store) 
-result_GEOG.dt.to_zarr(zstore, mode="w", encoding = pyramid.encoding, zarr_format=3)
+pyramid.dt.to_zarr(zstore, mode="w", encoding = pyramid.encoding, zarr_format=3)
 ```
 
+```python
+# Optional: Write to Icechunk
+# !pip install icechunk 
+import icechunk
+
+storage = icechunk.s3_storage(bucket="<add_your_bucket_name>", prefix="<add_your_prefix>", from_env=True)
+repo = icechunk.Repository.create(storage)
+session = repo.writable_session("main")
+
+store = from_url(url = "<add_your_bucket_url>", region="<add_your_region>")
+zstore = ObjectStore(store) 
+pyramid.dt.to_zarr(session.store, mode="w", encoding = pyramid.encoding, consolidated=False)
+```
+
+
 
 ## Development
 

pyproject.toml

@@ -5,7 +5,7 @@ description = "Lightweight multiscale zarr"
 readme = "README.md"
 requires-python = ">=3.11"
 dependencies = [
-    "dask>=2025.12.0",
+    "dask>=2025.9.0",
     "xarray>=2025.9.0",
     "xproj>=0.2.1",
 ]
@@ -25,3 +25,4 @@ test = [
     "pytest-xdist>=3.8.0",
     "ruff>=0.14.11",
 ]
+tests = []

src/topozarr/coarsen.py

@@ -2,7 +2,7 @@
 import xarray as xr
 from xarray import DataTree
 import xproj  # noqa ignore
-from .metadata import create_level_encoding, create_multiscale_metadata, SpecType
+from .metadata import create_level_encoding, create_multiscale_metadata
 from .pyramid import Pyramid
 from .chunking import DEFAULT_CHUNK_BYTES, DEFAULT_SHARD_BYTES
 
@@ -24,7 +24,6 @@ def build_coarsened_levels(
     x_dim: str,
     y_dim: str,
     method: CoarseningMethod,
-    spec: SpecType = "ndpyramid",
 ) -> dict[int, xr.Dataset]:
     levels = [ds]
     for lvl in range(num_levels - 1):
@@ -34,10 +33,8 @@ def build_coarsened_levels(
         coarsened = curr.coarsen({x_dim: 2, y_dim: 2}, boundary="trim")
         levels.insert(0, getattr(coarsened, method)())
 
-    if spec == "ndpyramid":
-        return dict(enumerate(levels))
-    else:
-        return dict(enumerate(reversed(levels)))
+    # zarr-multiscales: lowest levels = highest resolution (level 0 = highest res)
+    return dict(enumerate(reversed(levels)))
 
 
 def create_pyramid(
@@ -46,12 +43,11 @@ def create_pyramid(
     x_dim: str = "x",
     y_dim: str = "y",
     method: CoarseningMethod = "mean",
-    spec: SpecType = "ndpyramid",
     target_chunk_bytes: int = DEFAULT_CHUNK_BYTES,
-    target_shard_bytes: int = DEFAULT_SHARD_BYTES,
+    target_shard_bytes: int | None = DEFAULT_SHARD_BYTES,
 ) -> Pyramid:
     crs_str = get_crs(ds)
-    level_datasets = build_coarsened_levels(ds, levels, x_dim, y_dim, method, spec=spec)
+    level_datasets = build_coarsened_levels(ds, levels, x_dim, y_dim, method)
 
     dt = DataTree(name="root")
     full_encoding = {}
@@ -71,20 +67,20 @@ def create_pyramid(
         dim_chunks = {}
         for var_name, var_enc in level_encoding.items():
             if var_name in ds_level.data_vars and "chunks" in var_enc:
-                target_shards = var_enc["shards"]
+                dask_chunks = var_enc.get("shards", var_enc["chunks"])
                 da = ds_level[var_name]
 
-                for dim, shard_size in zip(da.dims, target_shards):
+                for dim, chunk_size in zip(da.dims, dask_chunks):
                     if dim not in dim_chunks:
-                        dim_chunks[dim] = shard_size
+                        dim_chunks[dim] = chunk_size
                     else:
-                        dim_chunks[dim] = min(dim_chunks[dim], shard_size)
+                        dim_chunks[dim] = min(dim_chunks[dim], chunk_size)
 
         if dim_chunks:
             ds_level = ds_level.chunk(dim_chunks)
 
         dt[path] = DataTree(ds_level, name=name)
         full_encoding[path] = level_encoding
 
-    dt.attrs = create_multiscale_metadata(levels, crs_str, method, spec=spec)
+    dt.attrs = create_multiscale_metadata(levels, crs_str, method)
     return Pyramid(datatree=dt, encoding=full_encoding)

src/topozarr/metadata.py

@@ -1,4 +1,4 @@
-from typing import Any, Literal
+from typing import Any
 import xarray as xr
 from .chunking import (
     calculate_chunk_size,
@@ -9,76 +9,103 @@
 )
 
 
-SpecType = Literal["ndpyramid", "zarr-multiscales"]
-
-
 def create_level_encoding(
     ds: xr.Dataset,
     x_dim: str,
     y_dim: str,
     target_chunk_bytes: int = DEFAULT_CHUNK_BYTES,
-    target_shard_bytes: int = DEFAULT_SHARD_BYTES,
+    target_shard_bytes: int | None = DEFAULT_SHARD_BYTES,
 ) -> dict[str, Any]:
-    encoding = {}
-    for var_name, da in ds.data_vars.items():
-        if x_dim not in da.dims or y_dim not in da.dims:
-            continue
+    spatial_vars = {
+        var_name: da
+        for var_name, da in ds.data_vars.items()
+        if x_dim in da.dims and y_dim in da.dims
+    }
 
-        itemsize = da.dtype.itemsize
-        ideal_chunk = get_ideal_dim(itemsize, target_chunk_bytes)
-        ideal_shard = get_ideal_dim(itemsize, target_shard_bytes)
+    return {
+        var_name: _create_var_encoding(
+            da, x_dim, y_dim, target_chunk_bytes, target_shard_bytes
+        )
+        for var_name, da in spatial_vars.items()
+    }
 
-        y_idx, x_idx = da.get_axis_num(y_dim), da.get_axis_num(x_dim)
 
-        chunks = list(da.shape)
-        shards = list(da.shape)
+def _create_var_encoding(
+    da: xr.DataArray,
+    x_dim: str,
+    y_dim: str,
+    target_chunk_bytes: int,
+    target_shard_bytes: int | None,
+) -> dict[str, Any]:
+    itemsize = da.dtype.itemsize
+    ideal_chunk = get_ideal_dim(itemsize, target_chunk_bytes)
+
+    y_idx, x_idx = da.get_axis_num(y_dim), da.get_axis_num(x_dim)
 
-        for idx, dim_name in [(y_idx, y_dim), (x_idx, x_dim)]:
-            c = calculate_chunk_size(da.shape[idx], ideal_chunk)
-            chunks[idx] = c
+    chunks = list(da.shape)
+    shards = list(da.shape) if target_shard_bytes is not None else None
+
+    for idx, dim_name in [(y_idx, y_dim), (x_idx, x_dim)]:
+        c = calculate_chunk_size(da.shape[idx], ideal_chunk)
+        chunks[idx] = c
+
+        if shards is not None:
+            ideal_shard = get_ideal_dim(itemsize, target_shard_bytes)
             shards[idx] = calculate_shard_size(da.shape[idx], c, ideal_shard)
 
-        for i, dim in enumerate(da.dims):
-            if dim not in [x_dim, y_dim]:
-                chunks[i] = 1
+    for i, dim in enumerate(da.dims):
+        if dim not in [x_dim, y_dim]:
+            chunks[i] = 1
+            if shards is not None:
                 shards[i] = 1
 
-        encoding[var_name] = {"chunks": tuple(chunks), "shards": tuple(shards)}
-    return encoding
+    var_encoding = {"chunks": tuple(chunks)}
+    if shards is not None:
+        var_encoding["shards"] = tuple(shards)
+
+    return var_encoding
 
 
 def create_multiscale_metadata(
-    levels: int, crs: str, method: str, spec: SpecType = "ndpyramid"
+    levels: int,
+    crs: str,
+    method: str,
 ) -> dict[str, Any]:
-    indices = list(range(levels))
-
-    if spec == "ndpyramid":
-        # ndpyramid-ish (highest levels = highest resolution)
-        datasets = [{"path": str(i), "level": i, "crs": crs} for i in reversed(indices)]
-        return {
-            "multiscales": [
-                {
-                    "datasets": datasets,
-                    "type": "reduce",
-                    "metadata": {"method": "coarsen", "coarsening_method": method},
-                }
-            ]
-        }
+    layout = []
 
-    # zarr-multiscales (as of Jan 8th 2026) (lowest levels = highest resolution)
-    layout = [
-        {
+    for i in range(levels):
+        entry = {
             "asset": str(i),
             "transform": {
                 "scale": [float(2**i), float(2**i)],
-                "translation": [0.0, 0.0],
+                "translation": [0.5, 0.5] if i > 0 else [0.0, 0.0],
             },
-            **({"derived_from": str(i - 1)} if i > 0 else {}),
         }
-        for i in indices
-    ]
 
+        if i > 0:
+            entry["derived_from"] = str(i - 1)
+            entry["resampling_method"] = method
+
+        layout.append(entry)
+
+    # attempting to match this example: https://github.com/zarr-conventions/multiscales/blob/main/examples/array-based-pyramid.json
     return {
+        "zarr_conventions": [
+            {
+                "schema_url": "https://raw.githubusercontent.com/zarr-conventions/multiscales/refs/tags/v1/schema.json",
+                "spec_url": "https://github.com/zarr-conventions/multiscales/blob/v1/README.md",
+                "uuid": "d35379db-88df-4056-af3a-620245f8e347",
+                "name": "multiscales",
+                "description": "Multiscale layout of zarr datasets",
+            },
+            {
+                "schema_url": "https://raw.githubusercontent.com/zarr-experimental/geo-proj/refs/tags/v1/schema.json",
+                "spec_url": "https://github.com/zarr-experimental/geo-proj/blob/v1/README.md",
+                "uuid": "f17cb550-5864-4468-aeb7-f3180cfb622f",
+                "name": "proj:",
+                "description": "Coordinate reference system information for geospatial data",
+            },
+        ],
         "multiscales": {"layout": layout, "resampling_method": method},
         "proj:code": crs,
     }

tests/test_chunking.py

@@ -77,6 +77,17 @@ def test_dask_chunks_match_shard_encoding(create_dataset):
                 )
 
 
+def test_disable_sharding(create_dataset):
+    from topozarr.coarsen import create_pyramid
+
+    ds = create_dataset(nx=1000, ny=1000)
+    pyramid = create_pyramid(ds, levels=1, target_shard_bytes=None)
+    enc = pyramid.encoding["/0"]["elevation"]
+
+    assert "chunks" in enc
+    assert "shards" not in enc
+
+
 def test_calculate_shard_size():
     """calculate_shard_size must return values divisible by chunk_size and less than or equal to the dim_size."""
     test_cases = [

tests/test_hypothesis.py

@@ -1,21 +1,20 @@
 import numpy as np
 import pytest
 import xarray as xr
-from hypothesis import given, strategies as st
+from hypothesis import given, settings, strategies as st
 from topozarr.coarsen import create_pyramid
 
-spatial_names = st.sampled_from(["x", "y", "lon", "lat"])
-extra_names = st.sampled_from(["time", "band"])
+spatial_names = st.sampled_from(["x", "y", "lon", "lat", "X", "Y"])
+extra_names = st.sampled_from(["time", "band", "z"])
 
 
 @st.composite
 def heterogeneous_datasets(draw):
-    """Generates datasets with spatial and extra dims ( 1 to 10)."""
     x_n = draw(spatial_names)
     y_n = draw(spatial_names.filter(lambda x: x != x_n))
     nx, ny = draw(st.integers(1, 10)), draw(st.integers(1, 10))
 
-    extras = draw(st.dictionaries(extra_names, st.integers(1, 3), max_size=1))
+    extras = draw(st.dictionaries(extra_names, st.integers(1, 3), max_size=2))
 
     all_dims = list(extras.keys()) + [y_n, x_n]
     shape = tuple(list(extras.values()) + [ny, nx])
@@ -31,6 +30,25 @@ def heterogeneous_datasets(draw):
     return ds.proj.assign_crs(spatial_ref="EPSG:4326"), x_n, y_n
 
 
+@st.composite
+def multi_variable_datasets(draw):
+    x_n = draw(spatial_names)
+    y_n = draw(spatial_names.filter(lambda x: x != x_n))
+    nx, ny = draw(st.integers(4, 16)), draw(st.integers(4, 16))
+
+    num_vars = draw(st.integers(2, 4))
+    data_vars = {}
+    for i in range(num_vars):
+        data_vars[f"var_{i}"] = ([y_n, x_n], np.zeros((ny, nx), dtype="f4"))
+
+    ds = xr.Dataset(
+        data_vars,
+        coords={x_n: np.arange(nx), y_n: np.arange(ny)},
+    )
+    return ds.proj.assign_crs(spatial_ref="EPSG:4326"), x_n, y_n
+
+
+@settings(deadline=1000)
 @given(ds_info=heterogeneous_datasets(), levels=st.integers(1, 5))
 def test_pyramid_integration_robustness(ds_info, levels):
     ds, x_dim, y_dim = ds_info
@@ -46,7 +64,19 @@ def test_pyramid_integration_robustness(ds_info, levels):
         assert len(pyramid.dt.children) == levels
 
         for path in pyramid.encoding:
-            enc = pyramid.encoding[path]["elevation"]
-            for c, s in zip(enc["chunks"], enc["shards"]):
-                assert c >= 1
-                assert s % c == 0
+            for var_name, enc in pyramid.encoding[path].items():
+                for c, s in zip(enc["chunks"], enc["shards"]):
+                    assert c >= 1
+                    assert s % c == 0
+
+
+@settings(deadline=1000)
+@given(ds_info=multi_variable_datasets(), levels=st.integers(1, 3))
+def test_multi_variable_encoding(ds_info, levels):
+    ds, x_dim, y_dim = ds_info
+
+    pyramid = create_pyramid(ds, levels=levels, x_dim=x_dim, y_dim=y_dim)
+
+    for level_idx in range(levels):
+        level_encoding = pyramid.encoding[f"/{level_idx}"]
+        assert len(level_encoding) == len(ds.data_vars)