From cd967ab31ebc15887cc3b37cd248f50defb16d71 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Mon, 15 Jun 2026 15:26:39 +0000
Subject: [PATCH 01/11] phase 1 of issue #29

---
 src/zagg/config.py   | 134 +++++++++++++++++++---
 tests/test_config.py | 256 +++++++++++++++++++++++++++++++++++--------
 2 files changed, 331 insertions(+), 59 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index 9acd02a5..e4f56a1a 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -147,16 +147,12 @@ def validate_config(config: PipelineConfig) -> None:
         if grid["type"] == "rectilinear":
             for field in ("crs", "resolution", "bounds"):
                 if field not in grid:
-                    raise ValueError(
-                        f"output.grid.{field} is required for rectilinear grid"
-                    )
+                    raise ValueError(f"output.grid.{field} is required for rectilinear grid")
             if len(grid["bounds"]) != 4:
                 raise ValueError("output.grid.bounds must be [xmin, ymin, xmax, ymax]")
         layout = grid.get("layout")
         if layout is not None and layout not in ("dense", "fullsphere"):
-            raise ValueError(
-                f"output.grid.layout must be 'dense' or 'fullsphere' (got {layout!r})"
-            )
+            raise ValueError(f"output.grid.layout must be 'dense' or 'fullsphere' (got {layout!r})")
 
     # Validate bounds structure (optional)
     if config.bounds is not None:
@@ -184,16 +180,12 @@ def validate_config(config: PipelineConfig) -> None:
 
         # Must have one (count via function:len is allowed)
         if not has_func and not has_expr:
-            raise ValueError(
-                f"Variable '{name}': must specify 'function' or 'expression'"
-            )
+            raise ValueError(f"Variable '{name}': must specify 'function' or 'expression'")
 
         # Validate source references
         source = meta.get("source")
         if source is not None and source not in ds_vars:
-            raise ValueError(
-                f"Variable '{name}': source '{source}' not in data_source.variables"
-            )
+            raise ValueError(f"Variable '{name}': source '{source}' not in data_source.variables")
 
         # Validate function resolves
         if has_func:
@@ -217,6 +209,85 @@ def validate_config(config: PipelineConfig) -> None:
         if has_expr:
             _validate_expression_columns(name, meta["expression"], ds_vars)
 
+        # Validate the output-kind declaration (kind + trailing_shape + dtype)
+        _validate_output_kind(name, meta)
+
+
+# Recognized per-field output kinds. ``ragged`` (CSR) is Tier 2 and not yet
+# accepted; see issue #29.
+OUTPUT_KINDS = ("scalar", "vector")
+
+
+def _validate_output_kind(name: str, meta: dict) -> None:
+    """Validate a variable's non-scalar output declaration.
+
+    A field may declare ``kind`` (``scalar`` default, or ``vector``),
+    ``trailing_shape`` (required for ``vector``), and a ``fill``/``pad_weight``
+    sentinel. ``scalar`` fields need none of these and stay the default path.
+
+    Parameters
+    ----------
+    name : str
+        Variable name (for error messages).
+    meta : dict
+        The variable's aggregation metadata.
+
+    Raises
+    ------
+    ValueError
+        On any invalid output-kind declaration.
+    """
+    kind = meta.get("kind", "scalar")
+    if kind not in OUTPUT_KINDS:
+        raise ValueError(
+            f"Variable '{name}': output kind '{kind}' is not supported "
+            f"(allowed: {OUTPUT_KINDS}; 'ragged' is planned but not yet implemented)"
+        )
+
+    # dtype, when declared, must name a real numpy dtype (applies to all kinds).
+    if "dtype" in meta:
+        try:
+            np.dtype(meta["dtype"])
+        except TypeError as e:
+            raise ValueError(
+                f"Variable '{name}': dtype {meta['dtype']!r} is not a valid numpy dtype ({e})"
+            ) from e
+
+    has_trailing = "trailing_shape" in meta
+
+    if kind == "scalar":
+        if has_trailing:
+            raise ValueError(
+                f"Variable '{name}': 'trailing_shape' is only valid for kind 'vector', not 'scalar'"
+            )
+        return
+
+    # kind == "vector": trailing_shape is required and must be positive ints.
+    if not has_trailing:
+        raise ValueError(f"Variable '{name}': kind 'vector' requires 'trailing_shape'")
+    _validate_trailing_shape(name, meta["trailing_shape"])
+
+
+def _validate_trailing_shape(name: str, trailing_shape) -> None:
+    """Check a vector field's trailing_shape is a tuple of positive ints."""
+    if isinstance(trailing_shape, int):
+        dims: tuple = (trailing_shape,)
+    elif isinstance(trailing_shape, (list, tuple)):
+        dims = tuple(trailing_shape)
+    else:
+        raise ValueError(
+            f"Variable '{name}': 'trailing_shape' must be an int or a "
+            f"sequence of ints (got {trailing_shape!r})"
+        )
+    if not dims:
+        raise ValueError(f"Variable '{name}': 'trailing_shape' must have at least one dimension")
+    for dim in dims:
+        if not isinstance(dim, int) or isinstance(dim, bool) or dim < 1:
+            raise ValueError(
+                f"Variable '{name}': 'trailing_shape' entries must be positive "
+                f"integers (got {dim!r})"
+            )
+
 
 def _is_numeric(s: str) -> bool:
     """Check if a string is a numeric literal."""
@@ -304,11 +375,48 @@ def get_agg_fields(config: PipelineConfig) -> dict:
     Returns
     -------
     dict
-        ``{name: {function/expression, source, params, dtype, fill_value, ...}}``
+        ``{name: {function/expression, source, params, dtype, fill_value, ...}}``.
+        A field may also declare a non-scalar output (issue #29) via ``kind``
+        (``scalar`` default, or ``vector``) and ``trailing_shape``; use
+        :func:`get_output_signature` to read the normalized declaration.
     """
     return dict(config.aggregation.get("variables", {}))
 
 
+def get_output_signature(meta: dict) -> dict:
+    """Return the normalized non-scalar output signature for one agg field.
+
+    This is the single read point for a field's Option B declaration (issue
+    #29): its output ``kind``, the per-cell ``trailing_shape``, and ``dtype``.
+    Later phases (statistic eval, the per-shard container, and the grid
+    ``signature()``) consume this rather than re-parsing the raw metadata.
+
+    Parameters
+    ----------
+    meta : dict
+        A single variable's aggregation metadata (a value of
+        :func:`get_agg_fields`).
+
+    Returns
+    -------
+    dict
+        ``{"kind": str, "trailing_shape": tuple[int, ...], "dtype": str}``.
+        ``trailing_shape`` is ``()`` for scalar fields. ``dtype`` is the
+        declared dtype string, or ``None`` if unset.
+    """
+    kind = meta.get("kind", "scalar")
+    if kind == "vector":
+        ts = meta["trailing_shape"]
+        trailing_shape = (ts,) if isinstance(ts, int) else tuple(ts)
+    else:
+        trailing_shape = ()
+    return {
+        "kind": kind,
+        "trailing_shape": trailing_shape,
+        "dtype": meta.get("dtype"),
+    }
+
+
 def get_coords(config: PipelineConfig) -> list[str]:
     """Return coordinate column names from the aggregation config.
 
diff --git a/tests/test_config.py b/tests/test_config.py
index bea099bb..51bf948b 100644
--- a/tests/test_config.py
+++ b/tests/test_config.py
@@ -14,6 +14,7 @@
     get_child_order,
     get_coords,
     get_data_vars,
+    get_output_signature,
     get_store_path,
     load_config,
     load_config_from_dict,
@@ -26,12 +27,14 @@
 # Fixtures
 # ---------------------------------------------------------------------------
 
+
 @pytest.fixture
 def atl06_yaml(tmp_path):
     """Path to the built-in atl06.yaml (copied to tmp for load_config tests)."""
     from importlib import resources
 
     import zagg.configs
+
     ref = resources.files(zagg.configs).joinpath("atl06.yaml")
     text = ref.read_text(encoding="utf-8")
     p = tmp_path / "atl06.yaml"
@@ -46,16 +49,19 @@ def atl06_config():
 
 @pytest.fixture
 def synthetic_df():
-    return pd.DataFrame({
-        "h_li": np.array([120.5, 118.3, 122.1, 119.7, 121.0], dtype=np.float32),
-        "s_li": np.array([0.05, 0.10, 0.03, 0.08, 0.06], dtype=np.float32),
-    })
+    return pd.DataFrame(
+        {
+            "h_li": np.array([120.5, 118.3, 122.1, 119.7, 121.0], dtype=np.float32),
+            "s_li": np.array([0.05, 0.10, 0.03, 0.08, 0.06], dtype=np.float32),
+        }
+    )
 
 
 # ---------------------------------------------------------------------------
 # Loading
 # ---------------------------------------------------------------------------
 
+
 class TestLoading:
     def test_load_yaml(self, atl06_yaml):
         cfg = load_config(atl06_yaml)
@@ -77,6 +83,7 @@ def test_all_sections_present(self, atl06_config):
 # Default config
 # ---------------------------------------------------------------------------
 
+
 class TestDefaultConfig:
     def test_default_atl06(self):
         cfg = default_config("atl06")
@@ -92,6 +99,7 @@ def test_nonexistent_raises(self):
 # Function resolution
 # ---------------------------------------------------------------------------
 
+
 class TestResolveFunction:
     def test_min(self):
         assert resolve_function("min") is np.min
@@ -117,6 +125,7 @@ def test_nonexistent_raises(self):
 # Expression evaluation
 # ---------------------------------------------------------------------------
 
+
 class TestEvaluateExpression:
     def test_simple_expression(self):
         cols = {"x": np.array([1.0, 2.0, 3.0])}
@@ -149,6 +158,7 @@ def test_undefined_column(self):
 # Validation
 # ---------------------------------------------------------------------------
 
+
 class TestValidation:
     def test_atl06_validates(self, atl06_config):
         # Should not raise
@@ -157,9 +167,11 @@ def test_atl06_validates(self, atl06_config):
     def test_missing_source(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {
-                "bad": {"function": "min", "source": "nonexistent", "dtype": "float32"},
-            }},
+            aggregation={
+                "variables": {
+                    "bad": {"function": "min", "source": "nonexistent", "dtype": "float32"},
+                }
+            },
             output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
         )
         with pytest.raises(ValueError, match="source.*nonexistent"):
@@ -168,14 +180,16 @@ def test_missing_source(self):
     def test_missing_weights_column(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {
-                "bad": {
-                    "function": "average",
-                    "source": "h_li",
-                    "params": {"weights": "missing_col"},
-                    "dtype": "float32",
-                },
-            }},
+            aggregation={
+                "variables": {
+                    "bad": {
+                        "function": "average",
+                        "source": "h_li",
+                        "params": {"weights": "missing_col"},
+                        "dtype": "float32",
+                    },
+                }
+            },
             output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
         )
         with pytest.raises(ValueError, match="missing_col"):
@@ -184,12 +198,14 @@ def test_missing_weights_column(self):
     def test_expression_unknown_column(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {
-                "bad": {
-                    "expression": "unknown_col + 1",
-                    "dtype": "float32",
-                },
-            }},
+            aggregation={
+                "variables": {
+                    "bad": {
+                        "expression": "unknown_col + 1",
+                        "dtype": "float32",
+                    },
+                }
+            },
             output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
         )
         with pytest.raises(ValueError, match="unknown_col"):
@@ -198,14 +214,16 @@ def test_expression_unknown_column(self):
     def test_function_and_expression_mutual_exclusion(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {
-                "bad": {
-                    "function": "min",
-                    "expression": "np.min(h_li)",
-                    "source": "h_li",
-                    "dtype": "float32",
-                },
-            }},
+            aggregation={
+                "variables": {
+                    "bad": {
+                        "function": "min",
+                        "expression": "np.min(h_li)",
+                        "source": "h_li",
+                        "dtype": "float32",
+                    },
+                }
+            },
             output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
         )
         with pytest.raises(ValueError, match="mutually exclusive"):
@@ -214,16 +232,20 @@ def test_function_and_expression_mutual_exclusion(self):
     def test_neither_function_nor_expression(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {
-                "bad": {"source": "h_li", "dtype": "float32"},
-            }},
+            aggregation={
+                "variables": {
+                    "bad": {"source": "h_li", "dtype": "float32"},
+                }
+            },
             output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
         )
         with pytest.raises(ValueError, match="must specify"):
             validate_config(cfg)
 
     def test_missing_required_section(self):
-        cfg = PipelineConfig(data_source={}, aggregation={"variables": {}}, output={"grid": {"type": "x"}})
+        cfg = PipelineConfig(
+            data_source={}, aggregation={"variables": {}}, output={"grid": {"type": "x"}}
+        )
         with pytest.raises(ValueError, match="Missing required section"):
             validate_config(cfg)
 
@@ -232,6 +254,7 @@ def test_missing_required_section(self):
 # Helper accessors
 # ---------------------------------------------------------------------------
 
+
 class TestAccessors:
     def test_get_agg_fields(self, atl06_config):
         fields = get_agg_fields(atl06_config)
@@ -254,6 +277,7 @@ def test_get_data_vars(self, atl06_config):
 # Equivalence with calculate_cell_statistics
 # ---------------------------------------------------------------------------
 
+
 def _dispatch_config_stat(name, meta, df):
     """Compute a single statistic using config metadata, mirroring calculate_cell_statistics."""
     if "function" in meta:
@@ -273,8 +297,12 @@ def _dispatch_config_stat(name, meta, df):
             if isinstance(v, str) and v in df.columns:
                 resolved[k] = df[v].values
             elif isinstance(v, str) and any(c in v for c in df.columns):
-                ns = {"__builtins__": {}, "np": np, "numpy": np,
-                      **{c: df[c].values for c in df.columns}}
+                ns = {
+                    "__builtins__": {},
+                    "np": np,
+                    "numpy": np,
+                    **{c: df[c].values for c in df.columns},
+                }
                 resolved[k] = eval(v, ns)  # noqa: S307
             else:
                 resolved[k] = v
@@ -318,6 +346,7 @@ def test_config_matches_calculate_cell_statistics(self, atl06_config, synthetic_
 # Roundtrip: YAML -> PipelineConfig -> dict -> PipelineConfig
 # ---------------------------------------------------------------------------
 
+
 class TestRoundtrip:
     def test_dict_roundtrip(self, atl06_config):
         d = asdict(atl06_config)
@@ -328,10 +357,16 @@ def test_dict_roundtrip(self, atl06_config):
 
     def test_catalog_and_bounds_roundtrip(self):
         d = {
-            "data_source": {"variables": {"h_li": "/path"}, "reader": "h5coro",
-                            "groups": ["gt1l"], "coordinates": {"latitude": "/lat", "longitude": "/lon"}},
-            "aggregation": {"variables": {"count": {"function": "len", "source": "h_li", "dtype": "int32"}},
-                            "coordinates": {"cell_ids": {"dtype": "uint64"}}},
+            "data_source": {
+                "variables": {"h_li": "/path"},
+                "reader": "h5coro",
+                "groups": ["gt1l"],
+                "coordinates": {"latitude": "/lat", "longitude": "/lon"},
+            },
+            "aggregation": {
+                "variables": {"count": {"function": "len", "source": "h_li", "dtype": "int32"}},
+                "coordinates": {"cell_ids": {"dtype": "uint64"}},
+            },
             "output": {"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
             "catalog": "my_catalog.json",
             "bounds": {"temporal": {"start_date": "2024-01-01", "end_date": "2024-06-01"}},
@@ -345,6 +380,7 @@ def test_catalog_and_bounds_roundtrip(self):
 # Output config helpers
 # ---------------------------------------------------------------------------
 
+
 class TestOutputHelpers:
     def test_get_child_order(self, atl06_config):
         assert get_child_order(atl06_config) == 12
@@ -355,14 +391,24 @@ def test_get_child_order_missing(self):
             get_child_order(cfg)
 
     def test_get_store_path(self):
-        cfg = PipelineConfig(output={"store": "./test.zarr", "grid": {"type": "healpix", "parent_order": 6, "child_order": 12}})
+        cfg = PipelineConfig(
+            output={
+                "store": "./test.zarr",
+                "grid": {"type": "healpix", "parent_order": 6, "child_order": 12},
+            }
+        )
         assert get_store_path(cfg) == "./test.zarr"
 
     def test_get_store_path_none(self, atl06_config):
         assert get_store_path(atl06_config) is None
 
     def test_get_store_path_s3(self):
-        cfg = PipelineConfig(output={"store": "s3://bucket/prefix.zarr", "grid": {"type": "healpix", "parent_order": 6, "child_order": 12}})
+        cfg = PipelineConfig(
+            output={
+                "store": "s3://bucket/prefix.zarr",
+                "grid": {"type": "healpix", "parent_order": 6, "child_order": 12},
+            }
+        )
         assert get_store_path(cfg) == "s3://bucket/prefix.zarr"
 
 
@@ -370,11 +416,14 @@ def test_get_store_path_s3(self):
 # Output grid validation
 # ---------------------------------------------------------------------------
 
+
 class TestOutputGridValidation:
     def test_grid_must_be_dict(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}},
+            aggregation={
+                "variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}
+            },
             output={"grid": "healpix"},
         )
         with pytest.raises(ValueError, match="output.grid must be a mapping"):
@@ -383,7 +432,9 @@ def test_grid_must_be_dict(self):
     def test_grid_missing_type(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}},
+            aggregation={
+                "variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}
+            },
             output={"grid": {"child_order": 12}},
         )
         with pytest.raises(ValueError, match="output.grid.type"):
@@ -392,7 +443,9 @@ def test_grid_missing_type(self):
     def test_healpix_missing_child_order(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}},
+            aggregation={
+                "variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}
+            },
             output={"grid": {"type": "healpix"}},
         )
         with pytest.raises(ValueError, match="child_order"):
@@ -401,7 +454,9 @@ def test_healpix_missing_child_order(self):
     def test_healpix_missing_parent_order(self):
         cfg = PipelineConfig(
             data_source={"variables": {"h_li": "/path"}},
-            aggregation={"variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}},
+            aggregation={
+                "variables": {"c": {"function": "len", "source": "h_li", "dtype": "int32"}}
+            },
             output={"grid": {"type": "healpix", "child_order": 12}},
         )
         with pytest.raises(ValueError, match="parent_order"):
@@ -412,6 +467,7 @@ def test_healpix_missing_parent_order(self):
 # Bounds validation
 # ---------------------------------------------------------------------------
 
+
 class TestBoundsValidation:
     def test_valid_bounds(self, atl06_config):
         atl06_config.bounds = {
@@ -429,7 +485,10 @@ def test_spatial_only(self, atl06_config):
         validate_config(atl06_config)
 
     def test_unknown_bounds_key(self, atl06_config):
-        atl06_config.bounds = {"temporal": {"start_date": "2024-01-01", "end_date": "2024-06-01"}, "foo": "bar"}
+        atl06_config.bounds = {
+            "temporal": {"start_date": "2024-01-01", "end_date": "2024-06-01"},
+            "foo": "bar",
+        }
         with pytest.raises(ValueError, match="Unknown bounds keys"):
             validate_config(atl06_config)
 
@@ -441,3 +500,108 @@ def test_temporal_missing_dates(self, atl06_config):
     def test_none_bounds_ok(self, atl06_config):
         atl06_config.bounds = None
         validate_config(atl06_config)
+
+
+# ---------------------------------------------------------------------------
+# Non-scalar output kind declaration (issue #29, phase 1)
+# ---------------------------------------------------------------------------
+
+
+def _vector_config(var_meta: dict) -> PipelineConfig:
+    """Build a minimal config whose single agg variable carries ``var_meta``."""
+    cfg = PipelineConfig(
+        data_source={
+            "reader": "h5coro",
+            "groups": ["gt1l"],
+            "coordinates": {"latitude": "/lat", "longitude": "/lon"},
+            "variables": {"h_li": "/path"},
+        },
+        aggregation={"variables": {"hist": {"source": "h_li", **var_meta}}},
+        output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}},
+    )
+    return cfg
+
+
+class TestOutputKind:
+    def test_scalar_default_backward_compatible(self, atl06_config):
+        """Existing scalar configs declare no kind and still validate."""
+        validate_config(atl06_config)
+        for meta in get_agg_fields(atl06_config).values():
+            assert "kind" not in meta
+
+    def test_vector_int_trailing_shape(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": 64})
+        validate_config(cfg)
+
+    def test_vector_list_trailing_shape(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": [16, 2]})
+        validate_config(cfg)
+
+    def test_unknown_kind_rejected(self):
+        cfg = _vector_config({"function": "min", "kind": "matrix"})
+        with pytest.raises(ValueError, match="output kind 'matrix' is not supported"):
+            validate_config(cfg)
+
+    def test_ragged_not_yet_supported(self):
+        cfg = _vector_config({"function": "min", "kind": "ragged"})
+        with pytest.raises(ValueError, match="ragged.*not yet implemented"):
+            validate_config(cfg)
+
+    def test_vector_requires_trailing_shape(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector"})
+        with pytest.raises(ValueError, match="requires 'trailing_shape'"):
+            validate_config(cfg)
+
+    def test_scalar_rejects_trailing_shape(self):
+        cfg = _vector_config({"function": "min", "trailing_shape": 8})
+        with pytest.raises(ValueError, match="only valid for kind 'vector'"):
+            validate_config(cfg)
+
+    def test_trailing_shape_must_be_positive(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": 0})
+        with pytest.raises(ValueError, match="positive"):
+            validate_config(cfg)
+
+    def test_trailing_shape_rejects_bool(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": True})
+        with pytest.raises(ValueError, match="positive"):
+            validate_config(cfg)
+
+    def test_trailing_shape_rejects_empty(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": []})
+        with pytest.raises(ValueError, match="at least one dimension"):
+            validate_config(cfg)
+
+    def test_trailing_shape_bad_type(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": "64"})
+        with pytest.raises(ValueError, match="int or a sequence of ints"):
+            validate_config(cfg)
+
+    def test_invalid_dtype_rejected(self):
+        cfg = _vector_config({"function": "min", "dtype": "not_a_dtype"})
+        with pytest.raises(ValueError, match="not a valid"):
+            validate_config(cfg)
+
+    def test_valid_dtypes_accepted(self):
+        for dt in ("float32", "int32", "uint64", "float64"):
+            cfg = _vector_config({"function": "min", "dtype": dt})
+            validate_config(cfg)
+
+
+class TestGetOutputSignature:
+    def test_scalar_signature(self):
+        sig = get_output_signature({"function": "min", "dtype": "float32"})
+        assert sig == {"kind": "scalar", "trailing_shape": (), "dtype": "float32"}
+
+    def test_scalar_default_dtype_none(self):
+        sig = get_output_signature({"function": "min"})
+        assert sig == {"kind": "scalar", "trailing_shape": (), "dtype": None}
+
+    def test_vector_int_signature(self):
+        sig = get_output_signature({"kind": "vector", "trailing_shape": 64, "dtype": "float32"})
+        assert sig == {"kind": "vector", "trailing_shape": (64,), "dtype": "float32"}
+
+    def test_vector_list_signature(self):
+        sig = get_output_signature({"kind": "vector", "trailing_shape": [16, 2]})
+        assert sig["kind"] == "vector"
+        assert sig["trailing_shape"] == (16, 2)

From 77f4838b4584368c76ea1be9908316c4194c2395 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Mon, 15 Jun 2026 15:30:42 +0000
Subject: [PATCH 02/11] address self-review on issue #29 phase 1

---
 src/zagg/config.py   | 12 +++++++-----
 tests/test_config.py | 19 +++++++++++++++++++
 2 files changed, 26 insertions(+), 5 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index e4f56a1a..0dbfeac7 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -221,9 +221,10 @@ def validate_config(config: PipelineConfig) -> None:
 def _validate_output_kind(name: str, meta: dict) -> None:
     """Validate a variable's non-scalar output declaration.
 
-    A field may declare ``kind`` (``scalar`` default, or ``vector``),
-    ``trailing_shape`` (required for ``vector``), and a ``fill``/``pad_weight``
-    sentinel. ``scalar`` fields need none of these and stay the default path.
+    A field may declare ``kind`` (``scalar`` default, or ``vector``) and
+    ``trailing_shape`` (required for ``vector``). ``scalar`` fields need
+    neither and stay the default path. (A per-cell padding sentinel for
+    vectors is a later phase; see issue #29.)
 
     Parameters
     ----------
@@ -239,16 +240,17 @@ def _validate_output_kind(name: str, meta: dict) -> None:
     """
     kind = meta.get("kind", "scalar")
     if kind not in OUTPUT_KINDS:
+        allowed = ", ".join(OUTPUT_KINDS)
         raise ValueError(
             f"Variable '{name}': output kind '{kind}' is not supported "
-            f"(allowed: {OUTPUT_KINDS}; 'ragged' is planned but not yet implemented)"
+            f"(allowed: {allowed}; 'ragged' is planned but not yet implemented)"
         )
 
     # dtype, when declared, must name a real numpy dtype (applies to all kinds).
     if "dtype" in meta:
         try:
             np.dtype(meta["dtype"])
-        except TypeError as e:
+        except (TypeError, ValueError) as e:
             raise ValueError(
                 f"Variable '{name}': dtype {meta['dtype']!r} is not a valid numpy dtype ({e})"
             ) from e
diff --git a/tests/test_config.py b/tests/test_config.py
index 51bf948b..504a722f 100644
--- a/tests/test_config.py
+++ b/tests/test_config.py
@@ -577,6 +577,25 @@ def test_trailing_shape_bad_type(self):
         with pytest.raises(ValueError, match="int or a sequence of ints"):
             validate_config(cfg)
 
+    def test_trailing_shape_list_bad_element_type(self):
+        cfg = _vector_config(
+            {"function": "histogram", "kind": "vector", "trailing_shape": [16, "x"]}
+        )
+        with pytest.raises(ValueError, match="positive"):
+            validate_config(cfg)
+
+    def test_trailing_shape_list_zero_element(self):
+        cfg = _vector_config({"function": "histogram", "kind": "vector", "trailing_shape": [16, 0]})
+        with pytest.raises(ValueError, match="positive"):
+            validate_config(cfg)
+
+    def test_trailing_shape_list_bool_element(self):
+        cfg = _vector_config(
+            {"function": "histogram", "kind": "vector", "trailing_shape": [16, True]}
+        )
+        with pytest.raises(ValueError, match="positive"):
+            validate_config(cfg)
+
     def test_invalid_dtype_rejected(self):
         cfg = _vector_config({"function": "min", "dtype": "not_a_dtype"})
         with pytest.raises(ValueError, match="not a valid"):

From 6b0d35f090354e5a2e5fb9919f309348c814613d Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 00:40:48 +0000
Subject: [PATCH 03/11] phase 2 of issue #29

---
 src/zagg/processing.py   | 74 ++++++++++++++++++++++++++++++++++----
 tests/test_processing.py | 77 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 144 insertions(+), 7 deletions(-)

diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index 7b30f610..f90d96bb 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -16,7 +16,13 @@
 from zarr import config, open_array
 from zarr.abc.store import Store
 
-from zagg.config import PipelineConfig, default_config, get_agg_fields, get_data_vars
+from zagg.config import (
+    PipelineConfig,
+    default_config,
+    get_agg_fields,
+    get_data_vars,
+    get_output_signature,
+)
 from zagg.schema import ProcessingMetadata
 
 logger = logging.getLogger(__name__)
@@ -34,6 +40,18 @@ def _make_url_rewriter(driver: str | None):
     return lambda url: url.replace("s3://", "", 1)
 
 
+def _field_sentinel(meta: dict) -> float:
+    """Per-cell fill value for an agg field's empty/unused slots.
+
+    Mirrors how ``process_shard`` / :func:`_kernel_aggregate` seed their output
+    arrays: the schema-declared ``fill_value`` (default ``"NaN"`` -> ``np.nan``,
+    else the literal numeric fill). Used both for scalar empty cells and for the
+    padding of ``vector`` fields (issue #29 Option B).
+    """
+    fill_value = meta.get("fill_value", "NaN")
+    return np.nan if fill_value == "NaN" else fill_value
+
+
 def _group_columns(
     col_dict: dict[str, np.ndarray],
     cell_col: np.ndarray,
@@ -244,10 +262,7 @@ def calculate_cell_statistics(
 
     n_obs = len(next(iter(cell_data.values()))) if cell_data else 0
     if n_obs == 0:
-        return {
-            name: (0 if meta.get("function") in ("len", "count") else np.nan)
-            for name, meta in agg_fields.items()
-        }
+        return {name: _empty_cell_value(meta) for name, meta in agg_fields.items()}
 
     result = {}
     for name, meta in agg_fields.items():
@@ -255,9 +270,18 @@ def calculate_cell_statistics(
         expression = meta.get("expression")
         source = meta.get("source") or value_col
         params = dict(meta.get("params", {}))
+        sig = get_output_signature(meta)
 
-        # Expression-based aggregation (e.g. h_sigma)
+        # Expression-based aggregation (e.g. h_sigma). Expressions are scalar-only
+        # in Tier 1: ``evaluate_expression`` returns a Python float. A vector
+        # expression field would need a non-scalar eval contract; that's deferred
+        # (issue #29) — declaring ``kind: vector`` with ``expression`` is rejected.
         if expression:
+            if sig["kind"] == "vector":
+                raise ValueError(
+                    f"Variable '{name}': vector output from an expression is not yet "
+                    f"supported; use 'function' (issue #29 Tier 1)"
+                )
             result[name] = evaluate_expression(expression, cell_data)
             continue
 
@@ -285,11 +309,47 @@ def calculate_cell_statistics(
                 resolved_params[pkey] = pval
 
         func = resolve_function(func_name)
-        result[name] = float(func(values, **resolved_params))
+        out = func(values, **resolved_params)
+        # Scalar fields stay byte-for-byte identical to the pre-#29 path; only a
+        # declared ``vector`` field is allowed to return an ndarray (issue #29).
+        result[name] = _coerce_field_value(out, sig) if sig["kind"] == "vector" else float(out)
 
     return result
 
 
+def _empty_cell_value(meta: dict):
+    """Value emitted for a single agg field when its cell has no observations.
+
+    Scalar fields keep the pre-#29 contract: ``0`` for ``len``/``count``,
+    ``np.nan`` otherwise. A ``vector`` field (issue #29) instead gets a full
+    ``trailing_shape`` array filled with its schema-declared sentinel
+    (:func:`_field_sentinel`), so empty and populated cells emit the same shape.
+    """
+    sig = get_output_signature(meta)
+    if sig["kind"] == "vector":
+        dtype = np.dtype(sig["dtype"]) if sig["dtype"] is not None else np.dtype("float32")
+        return np.full(sig["trailing_shape"], _field_sentinel(meta), dtype=dtype)
+    return 0 if meta.get("function") in ("len", "count") else np.nan
+
+
+def _coerce_field_value(value, sig: dict) -> np.ndarray:
+    """Coerce a ``vector`` field's aggregation output to its declared signature.
+
+    The field's ``function`` must yield exactly ``trailing_shape`` values (issue
+    #29 Tier-1 fixed-width vectors; ragged/CSR is Tier 2; vector ``expression``
+    fields are deferred and rejected upstream). Returns
+    a contiguous array of the declared dtype (default ``float32``), so every cell
+    emits an identically-shaped slab the dense writer (phase 5) can stack.
+    """
+    dtype = np.dtype(sig["dtype"]) if sig["dtype"] is not None else np.dtype("float32")
+    arr = np.asarray(value, dtype=dtype)
+    if arr.shape != sig["trailing_shape"]:
+        raise ValueError(
+            f"vector field produced shape {arr.shape}, expected {sig['trailing_shape']}"
+        )
+    return arr
+
+
 # EXPERIMENTAL (phase 2b of #30) -----------------------------------------------
 # Dual aggregation contract
 # -------------------------
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 88fe881d..46bec01f 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -145,6 +145,83 @@ def test_numpy_nan_aware_functions(self):
         assert result["h_nanmin"] == 1.0
 
 
+class TestVectorOutputs:
+    """Issue #29 phase 2: a ``kind: vector`` field yields a per-cell ndarray of
+    its declared ``trailing_shape``/``dtype``; scalar fields are unchanged."""
+
+    @staticmethod
+    def _hist_config(bins=4, dtype="int64", fill_value=None):
+        from zagg.config import PipelineConfig
+
+        meta = {
+            "function": "np.bincount",
+            "source": "b",
+            "kind": "vector",
+            "trailing_shape": bins,
+            "dtype": dtype,
+            "params": {"minlength": bins},
+        }
+        if fill_value is not None:
+            meta["fill_value"] = fill_value
+        return PipelineConfig(aggregation={"variables": {"hist": meta}})
+
+    def test_vector_field_returns_declared_shape(self):
+        cfg = self._hist_config(bins=4)
+        data = {"b": np.array([0, 1, 1, 3])}
+        result = calculate_cell_statistics(data, config=cfg)
+        hist = result["hist"]
+        assert isinstance(hist, np.ndarray)
+        assert hist.shape == (4,)
+        assert hist.dtype == np.dtype("int64")
+        np.testing.assert_array_equal(hist, [1, 2, 0, 1])
+
+    def test_vector_empty_cell_gets_sentinel(self):
+        cfg = self._hist_config(bins=4, dtype="float32")
+        result = calculate_cell_statistics({"b": np.array([])}, config=cfg)
+        hist = result["hist"]
+        assert hist.shape == (4,)
+        assert np.all(np.isnan(hist))  # default fill_value "NaN"
+
+    def test_vector_empty_cell_numeric_sentinel(self):
+        cfg = self._hist_config(bins=3, dtype="int64", fill_value=0)
+        result = calculate_cell_statistics({"b": np.array([])}, config=cfg)
+        np.testing.assert_array_equal(result["hist"], [0, 0, 0])
+
+    def test_vector_wrong_width_raises(self):
+        cfg = self._hist_config(bins=2)  # but bincount yields width 4 below
+        with pytest.raises(ValueError, match="expected"):
+            calculate_cell_statistics({"b": np.array([0, 3])}, config=cfg)
+
+    def test_vector_expression_rejected(self):
+        """Tier 1 vectors come from ``function``; a vector ``expression`` field is
+        deferred (issue #29) and must fail loudly rather than silently scalarize."""
+        from zagg.config import PipelineConfig
+
+        cfg = PipelineConfig(
+            aggregation={
+                "variables": {
+                    "edges": {
+                        "expression": "np.array([np.min(h), np.max(h)])",
+                        "kind": "vector",
+                        "trailing_shape": 2,
+                        "dtype": "float32",
+                    }
+                }
+            }
+        )
+        with pytest.raises(ValueError, match="vector output from an expression"):
+            calculate_cell_statistics({"h": np.array([1.0, 5.0, 3.0])}, config=cfg)
+
+    def test_scalar_fields_unchanged_alongside_vector(self):
+        """Adding a vector field must not perturb scalar outputs in the same dict."""
+        scalar = calculate_cell_statistics(
+            {"h_li": np.array([1.0, 2.0, 3.0]), "s_li": np.array([0.1, 0.1, 0.1])}
+        )
+        assert isinstance(scalar["h_min"], float)
+        assert scalar["h_min"] == 1.0
+        assert scalar["count"] == 3
+
+
 class TestBuildGroups:
     def test_slice_counts_match_per_cell_mask(self):
         """_build_groups produces identical cell populations as the old boolean-mask loop."""

From 3fda3c4afe8c60221ef91f9c44f20b993cf77a3d Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 10:15:06 +0000
Subject: [PATCH 04/11] vector expression fields for issue #29

---
 src/zagg/config.py       | 44 ++++++++++++++++++++++++++++++-----
 src/zagg/processing.py   | 29 +++++++++++------------
 tests/test_config.py     | 22 ++++++++++++++++++
 tests/test_processing.py | 50 ++++++++++++++++++++++++++++++++++++----
 4 files changed, 119 insertions(+), 26 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index 0dbfeac7..4b8205d3 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -223,8 +223,9 @@ def _validate_output_kind(name: str, meta: dict) -> None:
 
     A field may declare ``kind`` (``scalar`` default, or ``vector``) and
     ``trailing_shape`` (required for ``vector``). ``scalar`` fields need
-    neither and stay the default path. (A per-cell padding sentinel for
-    vectors is a later phase; see issue #29.)
+    neither and stay the default path. A ``vector`` field may be driven by
+    either ``function`` or ``expression``; ``len``/``count`` are rejected for
+    ``vector`` (they short-circuit to a scalar count). See issue #29.
 
     Parameters
     ----------
@@ -269,6 +270,15 @@ def _validate_output_kind(name: str, meta: dict) -> None:
         raise ValueError(f"Variable '{name}': kind 'vector' requires 'trailing_shape'")
     _validate_trailing_shape(name, meta["trailing_shape"])
 
+    # ``len``/``count`` short-circuit to a scalar obs count in
+    # ``calculate_cell_statistics``; pairing them with kind 'vector' would
+    # silently emit a scalar, so reject the nonsensical combination.
+    if meta.get("function") in ("len", "count"):
+        raise ValueError(
+            f"Variable '{name}': function {meta['function']!r} produces a scalar "
+            f"count and cannot be combined with kind 'vector'"
+        )
+
 
 def _validate_trailing_shape(name: str, trailing_shape) -> None:
     """Check a vector field's trailing_shape is a tuple of positive ints."""
@@ -568,8 +578,12 @@ def get_output_region(config: PipelineConfig) -> str | None:
     return config.output.get("region")
 
 
-def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> float:
-    """Evaluate an expression string in a restricted namespace.
+def eval_expression_raw(expression: str, columns: dict[str, np.ndarray]):
+    """Evaluate an expression string in a restricted namespace, uncoerced.
+
+    Returns the expression's native value (a scalar, an ndarray, ...). Used by
+    vector ``expression`` fields (issue #29), which coerce the result through
+    ``_coerce_field_value`` rather than casting to ``float``.
 
     Parameters
     ----------
@@ -580,7 +594,8 @@ def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> floa
 
     Returns
     -------
-    float
+    Any
+        Whatever the expression evaluates to.
     """
     ns = {
         "__builtins__": {},
@@ -593,4 +608,21 @@ def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> floa
         "sum": sum,
         **columns,
     }
-    return float(eval(expression, ns))  # noqa: S307
+    return eval(expression, ns)  # noqa: S307
+
+
+def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> float:
+    """Evaluate an expression string in a restricted namespace.
+
+    Parameters
+    ----------
+    expression : str
+        Python expression using numpy and column variables.
+    columns : dict[str, np.ndarray]
+        Mapping of column names to arrays.
+
+    Returns
+    -------
+    float
+    """
+    return float(eval_expression_raw(expression, columns))
diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index f90d96bb..12af9919 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -254,7 +254,7 @@ def calculate_cell_statistics(
     dict
         Dictionary of statistics keyed by aggregation variable name.
     """
-    from zagg.config import evaluate_expression, resolve_function
+    from zagg.config import eval_expression_raw, evaluate_expression, resolve_function
 
     if config is None:
         config = default_config()
@@ -272,17 +272,16 @@ def calculate_cell_statistics(
         params = dict(meta.get("params", {}))
         sig = get_output_signature(meta)
 
-        # Expression-based aggregation (e.g. h_sigma). Expressions are scalar-only
-        # in Tier 1: ``evaluate_expression`` returns a Python float. A vector
-        # expression field would need a non-scalar eval contract; that's deferred
-        # (issue #29) — declaring ``kind: vector`` with ``expression`` is rejected.
+        # Expression-based aggregation (e.g. h_sigma). A scalar expression casts
+        # to a Python float; a ``kind: vector`` expression is coerced through the
+        # same ``_coerce_field_value``/``trailing_shape``/dtype path as a vector
+        # ``function`` field (issue #29).
         if expression:
             if sig["kind"] == "vector":
-                raise ValueError(
-                    f"Variable '{name}': vector output from an expression is not yet "
-                    f"supported; use 'function' (issue #29 Tier 1)"
-                )
-            result[name] = evaluate_expression(expression, cell_data)
+                out = eval_expression_raw(expression, cell_data)
+                result[name] = _coerce_field_value(out, sig)
+            else:
+                result[name] = evaluate_expression(expression, cell_data)
             continue
 
         values = cell_data[source]
@@ -335,11 +334,11 @@ def _empty_cell_value(meta: dict):
 def _coerce_field_value(value, sig: dict) -> np.ndarray:
     """Coerce a ``vector`` field's aggregation output to its declared signature.
 
-    The field's ``function`` must yield exactly ``trailing_shape`` values (issue
-    #29 Tier-1 fixed-width vectors; ragged/CSR is Tier 2; vector ``expression``
-    fields are deferred and rejected upstream). Returns
-    a contiguous array of the declared dtype (default ``float32``), so every cell
-    emits an identically-shaped slab the dense writer (phase 5) can stack.
+    The field's ``function`` or ``expression`` must yield exactly
+    ``trailing_shape`` values (issue #29 Tier-1 fixed-width vectors; ragged/CSR
+    is Tier 2). Returns a contiguous array of the declared dtype (default
+    ``float32``), so every cell emits an identically-shaped slab the dense
+    writer (phase 5) can stack.
     """
     dtype = np.dtype(sig["dtype"]) if sig["dtype"] is not None else np.dtype("float32")
     arr = np.asarray(value, dtype=dtype)
diff --git a/tests/test_config.py b/tests/test_config.py
index 48cfa7ec..da1fdc9d 100644
--- a/tests/test_config.py
+++ b/tests/test_config.py
@@ -597,6 +597,28 @@ def test_trailing_shape_list_bool_element(self):
         with pytest.raises(ValueError, match="positive"):
             validate_config(cfg)
 
+    def test_vector_expression_allowed(self):
+        """A vector field may be driven by an expression (issue #29)."""
+        cfg = _vector_config(
+            {
+                "expression": "np.array([np.min(h_li), np.max(h_li)])",
+                "kind": "vector",
+                "trailing_shape": 2,
+            }
+        )
+        validate_config(cfg)
+
+    def test_vector_len_rejected(self):
+        """``len`` short-circuits to a scalar count; kind 'vector' is nonsensical."""
+        cfg = _vector_config({"function": "len", "kind": "vector", "trailing_shape": 4})
+        with pytest.raises(ValueError, match="cannot be combined with kind 'vector'"):
+            validate_config(cfg)
+
+    def test_vector_count_rejected(self):
+        cfg = _vector_config({"function": "count", "kind": "vector", "trailing_shape": 4})
+        with pytest.raises(ValueError, match="cannot be combined with kind 'vector'"):
+            validate_config(cfg)
+
     def test_invalid_dtype_rejected(self):
         cfg = _vector_config({"function": "min", "dtype": "not_a_dtype"})
         with pytest.raises(ValueError, match="not a valid"):
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 46bec01f..8492d778 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -192,16 +192,56 @@ def test_vector_wrong_width_raises(self):
         with pytest.raises(ValueError, match="expected"):
             calculate_cell_statistics({"b": np.array([0, 3])}, config=cfg)
 
-    def test_vector_expression_rejected(self):
-        """Tier 1 vectors come from ``function``; a vector ``expression`` field is
-        deferred (issue #29) and must fail loudly rather than silently scalarize."""
+    @staticmethod
+    def _edges_config(fill_value=None):
+        """A ``kind: vector`` field driven by an ``expression`` (issue #29)."""
+        from zagg.config import PipelineConfig
+
+        meta = {
+            "expression": "np.array([np.min(h), np.max(h)])",
+            "source": "h",
+            "kind": "vector",
+            "trailing_shape": 2,
+            "dtype": "float32",
+        }
+        if fill_value is not None:
+            meta["fill_value"] = fill_value
+        return PipelineConfig(aggregation={"variables": {"edges": meta}})
+
+    def test_vector_expression_returns_declared_shape(self):
+        """A vector ``expression`` field coerces to its declared shape/dtype, the
+        same path a vector ``function`` field uses (issue #29)."""
+        cfg = self._edges_config()
+        result = calculate_cell_statistics({"h": np.array([1.0, 5.0, 3.0])}, config=cfg)
+        edges = result["edges"]
+        assert isinstance(edges, np.ndarray)
+        assert edges.shape == (2,)
+        assert edges.dtype == np.dtype("float32")
+        np.testing.assert_array_equal(edges, [1.0, 5.0])
+
+    def test_vector_expression_empty_cell_gets_sentinel(self):
+        """An empty cell emits the same shape filled with the fill_value sentinel."""
+        cfg = self._edges_config()
+        result = calculate_cell_statistics({"h": np.array([])}, config=cfg)
+        edges = result["edges"]
+        assert edges.shape == (2,)
+        assert np.all(np.isnan(edges))  # default fill_value "NaN"
+
+    def test_vector_expression_empty_cell_numeric_sentinel(self):
+        cfg = self._edges_config(fill_value=0)
+        result = calculate_cell_statistics({"h": np.array([])}, config=cfg)
+        np.testing.assert_array_equal(result["edges"], [0, 0])
+
+    def test_vector_expression_wrong_width_raises(self):
+        """An expression yielding the wrong width fails loudly, like the function case."""
         from zagg.config import PipelineConfig
 
         cfg = PipelineConfig(
             aggregation={
                 "variables": {
                     "edges": {
-                        "expression": "np.array([np.min(h), np.max(h)])",
+                        "expression": "np.array([np.min(h), np.max(h), np.mean(h)])",
+                        "source": "h",
                         "kind": "vector",
                         "trailing_shape": 2,
                         "dtype": "float32",
@@ -209,7 +249,7 @@ def test_vector_expression_rejected(self):
                 }
             }
         )
-        with pytest.raises(ValueError, match="vector output from an expression"):
+        with pytest.raises(ValueError, match="expected"):
             calculate_cell_statistics({"h": np.array([1.0, 5.0, 3.0])}, config=cfg)
 
     def test_scalar_fields_unchanged_alongside_vector(self):

From 79b89605a89bca00f1318a745575eaa0fc256703 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 10:20:59 +0000
Subject: [PATCH 05/11] fold review: privatize _eval_expression_raw + clarify
 empty-cell tests (#29)

---
 src/zagg/config.py       | 6 +++---
 src/zagg/processing.py   | 4 ++--
 tests/test_processing.py | 6 +++++-
 3 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index 4b8205d3..6685e67e 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -5,7 +5,7 @@
 from collections.abc import Callable
 from dataclasses import dataclass, field
 from importlib import resources
-from typing import NotRequired, TypedDict
+from typing import Any, NotRequired, TypedDict
 
 import numpy as np
 import yaml
@@ -578,7 +578,7 @@ def get_output_region(config: PipelineConfig) -> str | None:
     return config.output.get("region")
 
 
-def eval_expression_raw(expression: str, columns: dict[str, np.ndarray]):
+def _eval_expression_raw(expression: str, columns: dict[str, np.ndarray]) -> Any:
     """Evaluate an expression string in a restricted namespace, uncoerced.
 
     Returns the expression's native value (a scalar, an ndarray, ...). Used by
@@ -625,4 +625,4 @@ def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> floa
     -------
     float
     """
-    return float(eval_expression_raw(expression, columns))
+    return float(_eval_expression_raw(expression, columns))
diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index 12af9919..c937004f 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -254,7 +254,7 @@ def calculate_cell_statistics(
     dict
         Dictionary of statistics keyed by aggregation variable name.
     """
-    from zagg.config import eval_expression_raw, evaluate_expression, resolve_function
+    from zagg.config import _eval_expression_raw, evaluate_expression, resolve_function
 
     if config is None:
         config = default_config()
@@ -278,7 +278,7 @@ def calculate_cell_statistics(
         # ``function`` field (issue #29).
         if expression:
             if sig["kind"] == "vector":
-                out = eval_expression_raw(expression, cell_data)
+                out = _eval_expression_raw(expression, cell_data)
                 result[name] = _coerce_field_value(out, sig)
             else:
                 result[name] = evaluate_expression(expression, cell_data)
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 8492d778..abc10920 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -220,7 +220,10 @@ def test_vector_expression_returns_declared_shape(self):
         np.testing.assert_array_equal(edges, [1.0, 5.0])
 
     def test_vector_expression_empty_cell_gets_sentinel(self):
-        """An empty cell emits the same shape filled with the fill_value sentinel."""
+        """An empty cell short-circuits to the fill_value sentinel WITHOUT
+        evaluating the expression. ``_edges_config``'s expression is
+        ``np.array([np.min(h), np.max(h)])`` and ``np.min([])`` raises, so this
+        passing proves the empty-cell path never reaches the eval (issue #29)."""
         cfg = self._edges_config()
         result = calculate_cell_statistics({"h": np.array([])}, config=cfg)
         edges = result["edges"]
@@ -228,6 +231,7 @@ def test_vector_expression_empty_cell_gets_sentinel(self):
         assert np.all(np.isnan(edges))  # default fill_value "NaN"
 
     def test_vector_expression_empty_cell_numeric_sentinel(self):
+        """Empty cell short-circuits to a numeric sentinel (no expression eval)."""
         cfg = self._edges_config(fill_value=0)
         result = calculate_cell_statistics({"h": np.array([])}, config=cfg)
         np.testing.assert_array_equal(result["edges"], [0, 0])

From 4d5a8443438096a3c9dbe4a7676c0e762d1b0a7a Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 17:22:58 +0000
Subject: [PATCH 06/11] phase 3 of issue #29

---
 src/zagg/processing.py   | 147 +++++++++++++++++++++++++++++++++------
 src/zagg/runner.py       |  13 ++--
 tests/test_processing.py | 128 ++++++++++++++++++++++++++++++++++
 3 files changed, 262 insertions(+), 26 deletions(-)

diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index c937004f..fc3ed9d7 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -158,21 +158,88 @@ def _concat_and_group(all_reads, grid, handoff: str):
     return col_arrays, cell_to_slice, n_obs_total
 
 
+def _has_vector_fields(config: PipelineConfig) -> bool:
+    """Whether any aggregation field declares a non-scalar (``vector``) output.
+
+    A pure-scalar config keeps the unchanged pandas carrier; any ``vector`` field
+    (issue #29) routes the whole cell->table handoff through Arrow (see
+    :func:`_arrow_column`).
+    """
+    return any(
+        get_output_signature(meta)["kind"] == "vector"
+        for meta in get_agg_fields(config).values()
+    )
+
+
+def _arrow_column(block: np.ndarray, sig: dict):
+    """Build the Arrow column for one agg field from its per-cell stats block.
+
+    A scalar field's ``(n_cells,)`` block becomes a plain Arrow array (values
+    byte-for-byte identical to the pandas carrier). A ``vector`` field's
+    ``(n_cells, *trailing_shape)`` block becomes a ``FixedSizeList<C>`` column
+    (``C = prod(trailing_shape)``), so every cell carries an identically-sized
+    list. Keeping the vector path a list-carrier (rather than a bespoke 2-D
+    column) is what lets the future ragged t-digest slot in as a variable-length
+    ``List<FixedSizeList<2>>`` through the same seam (issue #29 Tier 2).
+    """
+    import pyarrow as pa
+
+    if sig["kind"] != "vector":
+        return pa.array(block)
+    width = int(np.prod(sig["trailing_shape"]))
+    flat = np.ascontiguousarray(block).reshape(-1)
+    return pa.FixedSizeListArray.from_arrays(pa.array(flat), width)
+
+
+def _build_output(stats_arrays, data_vars, agg_fields, grid, shard_key, use_arrow: bool):
+    """Assemble the per-shard output carrier from the per-cell stats blocks.
+
+    Returns a ``pandas.DataFrame`` for a pure-scalar config (unchanged) or a
+    ``pyarrow.Table`` when any ``vector`` field is present, in both cases with the
+    data-variable columns followed by the grid's per-cell coord columns.
+    """
+    if not use_arrow:
+        df_out = pd.DataFrame({var: stats_arrays[var] for var in data_vars})
+        for col_name, vals in grid.chunk_coords(shard_key).items():
+            df_out[col_name] = vals
+        return df_out
+
+    import pyarrow as pa
+
+    columns = {
+        var: _arrow_column(stats_arrays[var], get_output_signature(agg_fields[var]))
+        for var in data_vars
+    }
+    for col_name, vals in grid.chunk_coords(shard_key).items():
+        columns[col_name] = pa.array(np.asarray(vals))
+    return pa.table(columns)
+
+
+def _carrier_empty(carrier) -> bool:
+    """Whether a process_shard output carrier (DataFrame or Arrow table) is empty."""
+    if isinstance(carrier, pd.DataFrame):
+        return carrier.empty
+    return carrier.num_rows == 0
+
+
 def write_dataframe_to_zarr(
-    df_out: pd.DataFrame,
+    df_out,
     store: Store,
     *,
     grid,
     chunk_idx: tuple,
 ) -> Store:
-    """Write a per-shard DataFrame to an existing Zarr template.
+    """Write a per-shard output carrier to an existing Zarr template.
 
     Parameters
     ----------
-    df_out : pd.DataFrame
-        Coordinate + data-variable columns. Row count must equal
-        ``prod(grid.chunk_shape)``; rows are in the grid's canonical
-        chunk order (``grid.children(shard_key)``).
+    df_out : pandas.DataFrame or pyarrow.Table
+        Coordinate + data-variable columns. A ``pyarrow.Table`` is used when the
+        config declares any ``vector`` field (issue #29): its ``FixedSizeList``
+        columns carry the per-cell ``trailing_shape`` payload, written to a
+        Zarr array with a trailing dimension. Cell count must equal
+        ``prod(grid.chunk_shape)``; cells are in the grid's canonical chunk order
+        (``grid.children(shard_key)``).
     store : Store
         Zarr-compatible store with the template already written.
     grid : OutputGrid
@@ -187,20 +254,24 @@ def write_dataframe_to_zarr(
     Store
         The same store, with data written.
     """
-    if df_out.empty:
+    if _carrier_empty(df_out):
         return store
 
     expected_count = int(np.prod(grid.chunk_shape))
-    if len(df_out) != expected_count:
+    n_cells = len(df_out) if isinstance(df_out, pd.DataFrame) else df_out.num_rows
+    if n_cells != expected_count:
         raise ValueError(
-            f"Expected {expected_count} rows for chunk_shape={grid.chunk_shape}, got {len(df_out)}"
+            f"Expected {expected_count} rows for chunk_shape={grid.chunk_shape}, got {n_cells}"
         )
 
     chunk_idx = tuple(int(i) for i in chunk_idx)
-    for name, series in df_out.items():
-        values = series.values
-        if values.shape != grid.chunk_shape:
-            values = values.reshape(grid.chunk_shape)
+    for name, values in _iter_carrier_columns(df_out):
+        # Scalar columns reshape to the grid's chunk_shape; a vector column keeps
+        # its trailing payload dim(s), so the block (and target array) is
+        # (*chunk_shape, *trailing_shape). The cell count invariant is unchanged.
+        trailing = values.shape[1:]
+        values = values.reshape((*grid.chunk_shape, *trailing))
+        block_idx = chunk_idx + (0,) * len(trailing)
         with config.set({"async.concurrency": 128}):
             array = open_array(
                 store,
@@ -208,11 +279,36 @@ def write_dataframe_to_zarr(
                 zarr_format=3,
                 consolidated=False,
             )
-            array.set_block_selection(chunk_idx, values)
+            array.set_block_selection(block_idx, values)
 
     return store
 
 
+def _iter_carrier_columns(carrier):
+    """Yield ``(name, ndarray)`` for each column of a DataFrame or Arrow table.
+
+    Scalar columns yield a 1-D array; a ``FixedSizeList<C>`` Arrow column yields a
+    2-D ``(n_cells, C)`` array (the per-cell vector block), so the writer can map
+    it onto the Zarr trailing payload dimension (issue #29).
+    """
+    if isinstance(carrier, pd.DataFrame):
+        for name, series in carrier.items():
+            yield name, series.values
+        return
+
+    import pyarrow as pa
+
+    n_rows = carrier.num_rows
+    for name in carrier.column_names:
+        col = carrier.column(name).combine_chunks()
+        if pa.types.is_fixed_size_list(col.type):
+            width = col.type.list_size
+            flat = col.values.to_numpy(zero_copy_only=False)
+            yield name, flat.reshape(n_rows, width)
+        else:
+            yield name, col.to_numpy(zero_copy_only=False)
+
+
 def calculate_cell_statistics(
     cell_data: dict[str, np.ndarray],
     value_col: str = "h_li",
@@ -815,12 +911,16 @@ def process_shard(
         stats_arrays = {}
         for name in data_vars:
             meta = agg_fields[name]
+            # Vector fields (issue #29) get a per-cell (n_cells, *trailing_shape)
+            # block; scalars keep the 1-D (n_cells,) layout, unchanged. Either way
+            # ``stats_arrays[name][i] = value`` assigns the cell's result row.
+            shape = (n_cells, *get_output_signature(meta)["trailing_shape"])
             zarr_dtype = np.dtype(meta.get("dtype", "float32"))
             fill_value = meta.get("fill_value", "NaN")
             if fill_value == "NaN":
-                stats_arrays[name] = np.full(n_cells, np.nan, dtype=zarr_dtype)
+                stats_arrays[name] = np.full(shape, np.nan, dtype=zarr_dtype)
             else:
-                stats_arrays[name] = np.zeros(n_cells, dtype=zarr_dtype)
+                stats_arrays[name] = np.zeros(shape, dtype=zarr_dtype)
 
         # Per-cell observation slices (grouped above, carrier-agnostic).
         _empty: dict[str, np.ndarray] = {col: arr[:0] for col, arr in col_arrays.items()}
@@ -843,10 +943,17 @@ def process_shard(
 
     logger.info(f"  Statistics: {cells_with_data}/{n_cells} cells with data")
 
-    # Create output DataFrame: data_vars + grid-specific per-cell coord columns
-    df_out = pd.DataFrame({var: stats_arrays[var] for var in data_vars})
-    for col_name, vals in grid.chunk_coords(shard_key).items():
-        df_out[col_name] = vals
+    # Assemble the output carrier: a plain DataFrame for a pure-scalar config
+    # (unchanged), or a pyarrow.Table with FixedSizeList vector columns when any
+    # field declares a non-scalar output (issue #29). Scalars stay byte-identical.
+    df_out = _build_output(
+        stats_arrays,
+        data_vars,
+        get_agg_fields(config),
+        grid,
+        shard_key,
+        use_arrow=_has_vector_fields(config),
+    )
 
     duration = (datetime.now() - start_time).total_seconds()
     logger.info(f"Completed shard {parent_morton} in {duration:.1f}s")
diff --git a/src/zagg/runner.py b/src/zagg/runner.py
index 95fe86dc..6d7efc80 100644
--- a/src/zagg/runner.py
+++ b/src/zagg/runner.py
@@ -284,12 +284,13 @@ def _process_and_write(shard_key, chunk_idx, records, grid,
         config=config,
         driver=driver,
     )
-    if not df_out.empty:
-        write_dataframe_to_zarr(
-            df_out, zarr_store,
-            grid=grid,
-            chunk_idx=chunk_idx,
-        )
+    # write_dataframe_to_zarr no-ops on an empty carrier (DataFrame or Arrow
+    # table), so no carrier-specific emptiness check is needed here.
+    write_dataframe_to_zarr(
+        df_out, zarr_store,
+        grid=grid,
+        chunk_idx=chunk_idx,
+    )
     return metadata
 
 
diff --git a/tests/test_processing.py b/tests/test_processing.py
index abc10920..792d29fe 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -8,9 +8,13 @@
 from zagg.grids import HealpixGrid
 from zagg.processing import (
     KERNEL_RTOL,
+    _arrow_column,
     _build_groups,
+    _build_output,
     _concat_and_group,
     _group_columns,
+    _has_vector_fields,
+    _iter_carrier_columns,
     _kernel_able,
     _kernel_aggregate,
     calculate_cell_statistics,
@@ -737,6 +741,130 @@ def test_invalid_handoff_rejected(self):
             process_shard(grid, 0, ["s3://x"], s3_credentials={}, handoff="bogus")
 
 
+class TestVectorCarrier:
+    """Issue #29 phase 3: a config with any ``vector`` field routes the
+    cell->table handoff through Arrow (FixedSizeList vector columns), while a
+    pure-scalar config keeps the unchanged pandas carrier with byte-identical
+    scalar outputs."""
+
+    @staticmethod
+    def _scalar_cfg():
+        from zagg.config import PipelineConfig
+
+        return PipelineConfig(
+            data_source={"groups": ["g"]},
+            aggregation={
+                "variables": {
+                    "count": {"function": "len"},
+                    "h_min": {"function": "min", "source": "h_li", "dtype": "float32"},
+                }
+            },
+        )
+
+    @staticmethod
+    def _vector_cfg():
+        """``_scalar_cfg`` plus a vector ``hist`` field (FixedSizeList<3>)."""
+        from zagg.config import PipelineConfig
+
+        return PipelineConfig(
+            data_source={"groups": ["g"]},
+            aggregation={
+                "variables": {
+                    "count": {"function": "len"},
+                    "h_min": {"function": "min", "source": "h_li", "dtype": "float32"},
+                    "hist": {
+                        "function": "np.bincount",
+                        "source": "b",
+                        "kind": "vector",
+                        "trailing_shape": 3,
+                        "dtype": "int64",
+                        "fill_value": 0,
+                        "params": {"minlength": 3},
+                    },
+                }
+            }
+        )
+
+    def test_has_vector_fields(self):
+        assert not _has_vector_fields(self._scalar_cfg())
+        assert _has_vector_fields(self._vector_cfg())
+
+    def _run(self, monkeypatch, cfg):
+        """Drive process_shard on a canned read via the default (pandas) handoff;
+        the output carrier (pandas vs Arrow) is chosen by the config's field kinds,
+        independent of the input handoff."""
+        pytest.importorskip("pyarrow")
+        leaf_to_cell = {1: 10, 2: 10, 3: 20}
+        children = [10, 20]
+        grid = _KernelShardGrid(children, leaf_to_cell)
+        df = pd.DataFrame(
+            {
+                "h_li": np.array([1.0, 2.0, 5.0], dtype=np.float32),
+                "b": np.array([0, 2, 1], dtype=np.int64),
+                "leaf_id": np.array([1, 1, 3], dtype=np.int64),
+            }
+        )
+        TestProcessShardKernelBranch()._patch_reads(monkeypatch, [df])
+        return process_shard(grid, 0, ["s3://x"], s3_credentials={}, config=cfg), children
+
+    def test_scalar_config_returns_dataframe(self, monkeypatch):
+        (df, _meta), _children = self._run(monkeypatch, self._scalar_cfg())
+        assert isinstance(df, pd.DataFrame)
+
+    def test_vector_config_returns_arrow_table(self, monkeypatch):
+        pa = pytest.importorskip("pyarrow")
+        (tbl, _meta), _children = self._run(monkeypatch, self._vector_cfg())
+        assert isinstance(tbl, pa.Table)
+        assert pa.types.is_fixed_size_list(tbl.column("hist").type)
+        assert tbl.column("hist").type.list_size == 3
+
+    def test_scalar_columns_byte_identical_with_and_without_vector(self, monkeypatch):
+        """The hard #29/#30 criterion: adding a vector field must not perturb the
+        scalar columns. Run the same canned input through both configs and assert
+        the shared scalar columns match exactly."""
+        (df, _m1), _c = self._run(monkeypatch, self._scalar_cfg())
+        (tbl, _m2), _c = self._run(monkeypatch, self._vector_cfg())
+        for name in ("count", "h_min"):
+            np.testing.assert_array_equal(
+                df[name].to_numpy(),
+                tbl.column(name).to_numpy(zero_copy_only=False),
+                err_msg=name,
+            )
+
+    def test_vector_column_values(self, monkeypatch):
+        """The FixedSizeList payload holds each cell's per-cell vector. cell 10 has
+        b=[0,2] -> bincount(minlength=3)=[1,0,1]; cell 20 has b=[1] -> [0,1,0]."""
+        (tbl, _meta), children = self._run(monkeypatch, self._vector_cfg())
+        hist = tbl.column("hist").combine_chunks()
+        block = hist.values.to_numpy(zero_copy_only=False).reshape(len(children), 3)
+        idx = {c: i for i, c in enumerate(children)}
+        np.testing.assert_array_equal(block[idx[10]], [1, 0, 1])
+        np.testing.assert_array_equal(block[idx[20]], [0, 1, 0])
+
+    def test_arrow_column_roundtrips_through_iter(self):
+        """_arrow_column -> _iter_carrier_columns recovers the (n_cells, C) block,
+        the seam the dense vector writer consumes (phase 5)."""
+        pa = pytest.importorskip("pyarrow")
+        sig = {"kind": "vector", "trailing_shape": (3,), "dtype": "int64"}
+        block = np.array([[1, 0, 1], [0, 1, 0]], dtype=np.int64)
+        col = _arrow_column(block, sig)
+        assert pa.types.is_fixed_size_list(col.type)
+        tbl = pa.table({"hist": col})
+        recovered = dict(_iter_carrier_columns(tbl))["hist"]
+        np.testing.assert_array_equal(recovered, block)
+
+    def test_build_output_scalar_is_plain_dataframe(self):
+        """_build_output(use_arrow=False) is the unchanged pandas assembly."""
+        grid = _KernelShardGrid([10, 20], {1: 10})
+        stats = {"count": np.array([2, 1]), "h_min": np.array([1.0, 5.0], dtype=np.float32)}
+        cfg = self._scalar_cfg()
+        out = _build_output(
+            stats, ["count", "h_min"], get_agg_fields(cfg), grid, 0, use_arrow=False
+        )
+        assert isinstance(out, pd.DataFrame)
+        np.testing.assert_array_equal(out["count"].to_numpy(), [2, 1])
+
+
 class TestDataSource:
     """Test data_source section of default config (replaces old DataSourceConfig tests)."""
 

From 93ee562b3222928a45e68b6625efccab72e91571 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 19:16:13 +0000
Subject: [PATCH 07/11] phase 4 of issue #29

---
 src/zagg/config.py            | 34 +++++++++++++
 src/zagg/grids/base.py        |  6 ++-
 src/zagg/grids/healpix.py     | 18 +++++--
 src/zagg/grids/rectilinear.py | 12 ++++-
 tests/test_grids.py           | 90 +++++++++++++++++++++++++++++++++++
 5 files changed, 154 insertions(+), 6 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index 6685e67e..2703cec7 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -429,6 +429,40 @@ def get_output_signature(meta: dict) -> dict:
     }
 
 
+def output_field_signature(config: PipelineConfig) -> list[dict]:
+    """Return the Option-B output-field signature for a config (issue #29).
+
+    A canonical, JSON-serializable list of ``{"name", "kind", "trailing_shape",
+    "dtype"}`` for every aggregation variable, sorted by ``name``. Recorded in a
+    grid's :meth:`signature` so a shard map can never be silently paired with a
+    grid whose output schema (scalar vs vector, trailing shape, dtype) differs,
+    and compared in ``nests_with`` so co-aggregated grids must share a field set.
+
+    ``trailing_shape`` is rendered as a ``list`` (``()``  for scalar fields) so
+    the structure round-trips through JSON unchanged.
+
+    Parameters
+    ----------
+    config : PipelineConfig
+
+    Returns
+    -------
+    list of dict
+    """
+    fields = []
+    for name, meta in get_agg_fields(config).items():
+        sig = get_output_signature(meta)
+        fields.append(
+            {
+                "name": name,
+                "kind": sig["kind"],
+                "trailing_shape": list(sig["trailing_shape"]),
+                "dtype": sig["dtype"],
+            }
+        )
+    return sorted(fields, key=lambda f: f["name"])
+
+
 def get_coords(config: PipelineConfig) -> list[str]:
     """Return coordinate column names from the aggregation config.
 
diff --git a/src/zagg/grids/base.py b/src/zagg/grids/base.py
index d3c40d63..b1ed7bef 100644
--- a/src/zagg/grids/base.py
+++ b/src/zagg/grids/base.py
@@ -124,8 +124,10 @@ def nests_with(self, other: "OutputGrid") -> bool:
         """Whether this grid and ``other`` tile compatibly (align + nest).
 
         The primitive for cross-aggregator compatibility validation: same
-        family, aligned, whole-number resolution ratios. Cross-family
-        (HEALPix vs rectilinear) always returns False.
+        family, aligned, whole-number resolution ratios, and the same Option-B
+        output-field set (issue #29 — same scalar/vector kinds, trailing
+        shapes, dtypes). Cross-family (HEALPix vs rectilinear) always returns
+        False.
         """
         ...
 
diff --git a/src/zagg/grids/healpix.py b/src/zagg/grids/healpix.py
index 92d0f45c..2138f2ba 100644
--- a/src/zagg/grids/healpix.py
+++ b/src/zagg/grids/healpix.py
@@ -8,7 +8,12 @@
 from zarr import config as zarr_config
 from zarr.abc.store import Store
 
-from zagg.config import PipelineConfig, default_config, get_agg_fields
+from zagg.config import (
+    PipelineConfig,
+    default_config,
+    get_agg_fields,
+    output_field_signature,
+)
 from zagg.grids.base import InconsistentShardError
 
 HEALPIX_BASE_CELLS: int = 12
@@ -204,15 +209,22 @@ def signature(self) -> dict:
             "parent_order": self.parent_order,
             "child_order": self.child_order,
             "layout": self.layout,
+            "output_fields": output_field_signature(self.config),
         }
 
     def nests_with(self, other) -> bool:
         """Whether ``self`` and ``other`` tile compatibly.
 
         Any two HEALPix grids nest (the nested hierarchy subdivides 4-for-1 at
-        every order). Cross-family (e.g. rectilinear) never nests.
+        every order), provided they declare the same Option-B output-field set
+        (issue #29) — co-aggregated grids must produce the same scalar/vector
+        schema. Cross-family (e.g. rectilinear) never nests.
         """
-        return isinstance(other, HealpixGrid)
+        if not isinstance(other, HealpixGrid):
+            return False
+        return output_field_signature(self.config) == output_field_signature(
+            other.config
+        )
 
     def emit_template(self, store: Store, *, overwrite: bool = False) -> Store:
         """Write the Zarr template (group + arrays) to ``store``."""
diff --git a/src/zagg/grids/rectilinear.py b/src/zagg/grids/rectilinear.py
index 9dee5520..5a3a9e2d 100644
--- a/src/zagg/grids/rectilinear.py
+++ b/src/zagg/grids/rectilinear.py
@@ -41,7 +41,12 @@
 from zarr import config as zarr_config
 from zarr.abc.store import Store
 
-from zagg.config import PipelineConfig, default_config, get_agg_fields
+from zagg.config import (
+    PipelineConfig,
+    default_config,
+    get_agg_fields,
+    output_field_signature,
+)
 
 OOB_SENTINEL: int = -1
 
@@ -152,6 +157,7 @@ def signature(self) -> dict:
             "affine": [a.a, a.b, a.c, a.d, a.e, a.f],
             "shape": [self.height, self.width],
             "chunk_shape": [self.chunk_h, self.chunk_w],
+            "output_fields": output_field_signature(self.config),
         }
 
     def nests_with(self, other) -> bool:
@@ -165,6 +171,10 @@ def nests_with(self, other) -> bool:
             return False
         if self._geobox.crs != other._geobox.crs:
             return False
+        if output_field_signature(self.config) != output_field_signature(other.config):
+            # Co-aggregated grids must declare the same Option-B output-field
+            # set (issue #29): same scalar/vector kinds, trailing shapes, dtypes.
+            return False
         if not (_whole_ratio(self.res_x, other.res_x)
                 and _whole_ratio(self.res_y, other.res_y)):
             return False
diff --git a/tests/test_grids.py b/tests/test_grids.py
index 871b6931..448e4eea 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -245,3 +245,93 @@ def test_with_n_parent_cells_means_dense(self, cfg):
         )
         group = open_group(store, path="8", mode="r")
         assert group["count"].shape == (4 ** (8 - 6) * 3,)
+
+
+def _vector_config(bins=4, dtype="int64"):
+    """A config with one scalar (``count``) and one ``kind: vector`` field
+    (issue #29), reusing the atl06 coordinates so grids build normally."""
+    from zagg.config import PipelineConfig
+
+    base = default_config("atl06")
+    agg = {
+        "coordinates": base.aggregation.get("coordinates", {}),
+        "variables": {
+            "count": {"function": "len", "source": "h_li"},
+            "hist": {
+                "function": "np.bincount",
+                "source": "b",
+                "kind": "vector",
+                "trailing_shape": bins,
+                "dtype": dtype,
+            },
+        },
+    }
+    return PipelineConfig(
+        data_source=base.data_source, aggregation=agg, output=base.output
+    )
+
+
+class TestOutputFieldSignature:
+    """Issue #29 phase 4: ``signature()`` carries the Option-B output-field set
+    and ``nests_with()`` requires a matching set."""
+
+    def test_signature_includes_output_fields(self, cfg):
+        g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
+        sig = g.signature()
+        assert "output_fields" in sig
+        names = {f["name"] for f in sig["output_fields"]}
+        assert "count" in names
+        # Each entry carries the Option-B keys.
+        for f in sig["output_fields"]:
+            assert set(f) == {"name", "kind", "trailing_shape", "dtype"}
+
+    def test_signature_marks_vector_field(self):
+        g = HealpixGrid(
+            parent_order=6, child_order=8, layout="fullsphere", config=_vector_config()
+        )
+        by_name = {f["name"]: f for f in g.signature()["output_fields"]}
+        assert by_name["hist"]["kind"] == "vector"
+        assert by_name["hist"]["trailing_shape"] == [4]
+        assert by_name["count"]["kind"] == "scalar"
+        assert by_name["count"]["trailing_shape"] == []
+
+    def test_signature_is_json_serializable(self):
+        import json
+
+        g = HealpixGrid(
+            parent_order=6, child_order=8, layout="fullsphere", config=_vector_config()
+        )
+        # Round-trips unchanged (recorded in a ShardMap as JSON).
+        assert json.loads(json.dumps(g.signature())) == json.loads(
+            json.dumps(g.signature())
+        )
+
+    def test_nests_with_same_field_set(self, cfg):
+        a = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
+        b = HealpixGrid(parent_order=4, child_order=8, layout="fullsphere", config=cfg)
+        assert a.nests_with(b) and b.nests_with(a)
+
+    def test_nests_with_differing_field_kind_rejected(self, cfg):
+        scalar = HealpixGrid(
+            parent_order=6, child_order=8, layout="fullsphere", config=cfg
+        )
+        vector = HealpixGrid(
+            parent_order=6, child_order=8, layout="fullsphere", config=_vector_config()
+        )
+        assert not scalar.nests_with(vector)
+        assert not vector.nests_with(scalar)
+
+    def test_nests_with_differing_trailing_shape_rejected(self):
+        a = HealpixGrid(
+            parent_order=6,
+            child_order=8,
+            layout="fullsphere",
+            config=_vector_config(bins=4),
+        )
+        b = HealpixGrid(
+            parent_order=6,
+            child_order=8,
+            layout="fullsphere",
+            config=_vector_config(bins=8),
+        )
+        assert not a.nests_with(b)

From 9c736f9b8b84c88fbe0f0a05dc55cb4ee085dd25 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 19:20:57 +0000
Subject: [PATCH 08/11] phase 5 of issue #29

---
 src/zagg/config.py            |  2 +-
 src/zagg/grids/base.py        | 62 ++++++++++++++++++++++++++++++++++-
 src/zagg/grids/healpix.py     | 13 ++++++--
 src/zagg/grids/rectilinear.py | 12 ++++++-
 src/zagg/processing.py        |  5 +++
 tests/test_grids.py           | 55 ++++++++++++++++++++++++++++---
 tests/test_rectilinear.py     | 27 +++++++++++++++
 7 files changed, 167 insertions(+), 9 deletions(-)

diff --git a/src/zagg/config.py b/src/zagg/config.py
index 2703cec7..3d860e07 100644
--- a/src/zagg/config.py
+++ b/src/zagg/config.py
@@ -438,7 +438,7 @@ def output_field_signature(config: PipelineConfig) -> list[dict]:
     grid whose output schema (scalar vs vector, trailing shape, dtype) differs,
     and compared in ``nests_with`` so co-aggregated grids must share a field set.
 
-    ``trailing_shape`` is rendered as a ``list`` (``()``  for scalar fields) so
+    ``trailing_shape`` is rendered as a ``list`` (``()`` for scalar fields) so
     the structure round-trips through JSON unchanged.
 
     Parameters
diff --git a/src/zagg/grids/base.py b/src/zagg/grids/base.py
index b1ed7bef..11724fed 100644
--- a/src/zagg/grids/base.py
+++ b/src/zagg/grids/base.py
@@ -30,6 +30,61 @@
 ShardKey = Any  # int for HEALPix, tuple[int,int] for rectilinear, etc.
 
 
+def vector_array_spec(base, sig, *, base_dims, base_chunk_shape):
+    """Extend a scalar data-var ``ArraySpec`` with a trailing payload dim.
+
+    Issue #29 phase 5: a ``kind: vector`` field is stored as a dense array of
+    shape ``(*spatial_shape, *trailing_shape)`` — the per-cell vector rides on
+    one or more trailing dimensions appended to the grid's spatial axes.
+
+    **Single-trailing-chunk invariant.** The trailing payload dimension(s) are
+    chunked *whole* (one chunk spans the full ``trailing_shape``). This is what
+    lets :func:`zagg.processing.write_dataframe_to_zarr` address a shard's
+    payload with ``block_idx = chunk_idx + (0,) * len(trailing_shape)`` — the
+    trailing block index is always ``0``. Do not chunk the trailing dim; the
+    writer assumes block 0.
+
+    Parameters
+    ----------
+    base : ArraySpec
+        The scalar data-var spec for this grid (already carries the field's
+        ``dtype``/``fill_value``); its ``shape`` is the spatial shape and its
+        chunk grid the spatial chunk shape.
+    sig : dict
+        Output signature from :func:`zagg.config.get_output_signature`
+        (``kind``/``trailing_shape``/``dtype``). For a scalar field (empty
+        ``trailing_shape``) ``base`` is returned unchanged.
+    base_dims : tuple of str
+        Spatial dimension names (e.g. ``("cells",)`` or ``("y", "x")``).
+    base_chunk_shape : tuple of int
+        Spatial chunk shape (the grid's ``chunk_shape``).
+
+    Returns
+    -------
+    ArraySpec
+    """
+    from pydantic_zarr.experimental.v3 import NamedConfig
+
+    trailing = tuple(int(t) for t in sig["trailing_shape"])
+    if not trailing:
+        return base
+    shape = (*base.shape, *trailing)
+    dim_names = (*base_dims, *(f"{base_dims[-1]}_v{i}" for i in range(len(trailing))))
+    chunk_shape = (*base_chunk_shape, *trailing)
+    chunk_grid = NamedConfig(name="regular", configuration={"chunk_shape": list(chunk_shape)})
+    # Set shape + dimension_names together: ArraySpec validates their ranks
+    # match on construction, so a chained ``with_shape`` then
+    # ``with_dimension_names`` would transiently mismatch and raise.
+    return type(base)(
+        **{
+            **base.model_dump(),
+            "shape": shape,
+            "dimension_names": dim_names,
+            "chunk_grid": chunk_grid,
+        }
+    )
+
+
 class InconsistentShardError(ValueError):
     """Raised by shard_of when input cells don't all share a shard."""
 
@@ -132,4 +187,9 @@ def nests_with(self, other: "OutputGrid") -> bool:
         ...
 
 
-__all__ = ["OutputGrid", "ShardKey", "InconsistentShardError"]
+__all__ = [
+    "OutputGrid",
+    "ShardKey",
+    "InconsistentShardError",
+    "vector_array_spec",
+]
diff --git a/src/zagg/grids/healpix.py b/src/zagg/grids/healpix.py
index 2138f2ba..7e89c0cd 100644
--- a/src/zagg/grids/healpix.py
+++ b/src/zagg/grids/healpix.py
@@ -12,9 +12,10 @@
     PipelineConfig,
     default_config,
     get_agg_fields,
+    get_output_signature,
     output_field_signature,
 )
-from zagg.grids.base import InconsistentShardError
+from zagg.grids.base import InconsistentShardError, vector_array_spec
 
 HEALPIX_BASE_CELLS: int = 12
 HEALPIX_REF_ORDER: int = 18  # mortie's clip2order reference order; do not change
@@ -269,7 +270,15 @@ def _spec(self) -> GroupSpec:
         for name, meta in get_agg_fields(self.config).items():
             dtype = meta.get("dtype", "float32")
             fill = meta.get("fill_value", "NaN")
-            members[name] = base.with_data_type(dtype).with_fill_value(fill)
+            spec = base.with_data_type(dtype).with_fill_value(fill)
+            # A vector field (issue #29) gets a trailing payload dim chunked
+            # whole; scalars are returned unchanged.
+            members[name] = vector_array_spec(
+                spec,
+                get_output_signature(meta),
+                base_dims=("cells",),
+                base_chunk_shape=(self.n_children,),
+            )
 
         return GroupSpec(members=members, attributes=self._dggs_attrs())
 
diff --git a/src/zagg/grids/rectilinear.py b/src/zagg/grids/rectilinear.py
index 5a3a9e2d..86698d59 100644
--- a/src/zagg/grids/rectilinear.py
+++ b/src/zagg/grids/rectilinear.py
@@ -45,8 +45,10 @@
     PipelineConfig,
     default_config,
     get_agg_fields,
+    get_output_signature,
     output_field_signature,
 )
+from zagg.grids.base import vector_array_spec
 
 OOB_SENTINEL: int = -1
 
@@ -368,7 +370,15 @@ def _spec(self) -> GroupSpec:
         for name, meta in get_agg_fields(self.config).items():
             dtype = meta.get("dtype", "float32")
             fill = meta.get("fill_value", "NaN")
-            members[name] = base.with_data_type(dtype).with_fill_value(fill)
+            spec = base.with_data_type(dtype).with_fill_value(fill)
+            # A vector field (issue #29) gets a trailing payload dim chunked
+            # whole; scalars are returned unchanged.
+            members[name] = vector_array_spec(
+                spec,
+                get_output_signature(meta),
+                base_dims=("y", "x"),
+                base_chunk_shape=self.chunk_shape,
+            )
 
         return GroupSpec(members=members, attributes=self._geozarr_attrs())
 
diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index fc3ed9d7..7c34121e 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -269,6 +269,11 @@ def write_dataframe_to_zarr(
         # Scalar columns reshape to the grid's chunk_shape; a vector column keeps
         # its trailing payload dim(s), so the block (and target array) is
         # (*chunk_shape, *trailing_shape). The cell count invariant is unchanged.
+        #
+        # Single-trailing-chunk invariant (issue #29): the template
+        # (``grids.base.vector_array_spec``) chunks the trailing payload dim
+        # *whole*, so the trailing block index is always 0 and a shard's payload
+        # lands in one Zarr block via ``chunk_idx + (0,) * len(trailing)``.
         trailing = values.shape[1:]
         values = values.reshape((*grid.chunk_shape, *trailing))
         block_idx = chunk_idx + (0,) * len(trailing)
diff --git a/tests/test_grids.py b/tests/test_grids.py
index 448e4eea..d4a0d9ce 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -301,10 +301,9 @@ def test_signature_is_json_serializable(self):
         g = HealpixGrid(
             parent_order=6, child_order=8, layout="fullsphere", config=_vector_config()
         )
-        # Round-trips unchanged (recorded in a ShardMap as JSON).
-        assert json.loads(json.dumps(g.signature())) == json.loads(
-            json.dumps(g.signature())
-        )
+        fields = g.signature()["output_fields"]
+        # Round-trips through JSON unchanged (recorded in a ShardMap as JSON).
+        assert json.loads(json.dumps(fields)) == fields
 
     def test_nests_with_same_field_set(self, cfg):
         a = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
@@ -335,3 +334,51 @@ def test_nests_with_differing_trailing_shape_rejected(self):
             config=_vector_config(bins=8),
         )
         assert not a.nests_with(b)
+
+
+class TestVectorTemplate:
+    """Issue #29 phase 5: a ``kind: vector`` field's template array gets a
+    trailing payload dim chunked whole (single-trailing-chunk invariant)."""
+
+    def test_healpix_vector_array_has_trailing_dim(self):
+        cfg = _vector_config(bins=4, dtype="int64")
+        # int vector field needs an int fill_value (NaN is invalid on int dtype).
+        cfg.aggregation["variables"]["hist"]["fill_value"] = 0
+        g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
+        store = MemoryStore()
+        g.emit_template(store)
+        grp = open_group(store, path="8", mode="r")
+        n_pix = HEALPIX_BASE_CELLS * 4**8
+        assert grp["count"].shape == (n_pix,)  # scalar unchanged
+        assert grp["hist"].shape == (n_pix, 4)  # trailing payload dim
+        # Trailing dim is ONE chunk (block_idx invariant): chunk == full width.
+        assert grp["hist"].chunks == (4 ** (8 - 6), 4)
+        assert grp["hist"].dtype == np.dtype("int64")
+
+    def test_healpix_dimension_names_extend(self):
+        cfg = _vector_config(bins=3, dtype="int64")
+        cfg.aggregation["variables"]["hist"]["fill_value"] = 0
+        g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
+        spec = g.spec()
+        names = spec.members["hist"].dimension_names
+        assert names[0] == "cells"
+        assert len(names) == 2  # spatial + one trailing dim
+
+    def test_rectilinear_vector_array_has_trailing_dim(self):
+        from zagg.grids import RectilinearGrid
+
+        cfg = _vector_config(bins=4, dtype="int64")
+        cfg.aggregation["variables"]["hist"]["fill_value"] = 0
+        g = RectilinearGrid(
+            "EPSG:3031",
+            1000.0,
+            (-1e6, -1e6, 1e6, 1e6),
+            chunk_shape=(64, 64),
+            config=cfg,
+        )
+        store = MemoryStore()
+        g.emit_template(store)
+        grp = open_group(store, path="rectilinear", mode="r")
+        assert grp["count"].shape == (g.height, g.width)
+        assert grp["hist"].shape == (g.height, g.width, 4)
+        assert grp["hist"].chunks == (64, 64, 4)  # trailing dim whole
diff --git a/tests/test_rectilinear.py b/tests/test_rectilinear.py
index 8d70c862..48d522d8 100644
--- a/tests/test_rectilinear.py
+++ b/tests/test_rectilinear.py
@@ -295,6 +295,33 @@ def test_cross_family(self, grid):
 
         assert not grid.nests_with(HealpixGrid(6, 12, layout="fullsphere"))
 
+    def test_differing_output_field_set_rejected(self, cfg):
+        """Otherwise-nesting grids with a different Option-B output-field set
+        (issue #29 phase 4) must not nest, symmetrically."""
+        from zagg.config import PipelineConfig
+
+        agg = {
+            "coordinates": cfg.aggregation.get("coordinates", {}),
+            "variables": {
+                "count": {"function": "len", "source": "h_li"},
+                "hist": {
+                    "function": "np.bincount",
+                    "source": "b",
+                    "kind": "vector",
+                    "trailing_shape": 4,
+                    "dtype": "int64",
+                },
+            },
+        }
+        vcfg = PipelineConfig(
+            data_source=cfg.data_source, aggregation=agg, output=cfg.output
+        )
+        scalar = _grid(cfg, 5000, (256, 256))
+        vector = _grid(vcfg, 5000, (256, 256))
+        # Same CRS / whole-ratio / aligned — only the field set differs.
+        assert not scalar.nests_with(vector)
+        assert not vector.nests_with(scalar)
+
 
 class TestValidateCompatible:
     def test_compatible_pass(self, cfg):

From 1511853564cf9bd45c6c6dc518d51ce524118a10 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 19:23:38 +0000
Subject: [PATCH 09/11] phase 6 of issue #29

---
 tests/test_processing.py | 81 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 80 insertions(+), 1 deletion(-)

diff --git a/tests/test_processing.py b/tests/test_processing.py
index 792d29fe..0308857f 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -5,7 +5,7 @@
 from zarr.storage import MemoryStore
 
 from zagg.config import default_config, get_agg_fields, get_coords, get_data_vars
-from zagg.grids import HealpixGrid
+from zagg.grids import HEALPIX_BASE_CELLS, HealpixGrid
 from zagg.processing import (
     KERNEL_RTOL,
     _arrow_column,
@@ -894,3 +894,82 @@ def test_group_template_substitution(self):
         ds = default_config().data_source
         path = ds["coordinates"]["latitude"].format(group="gt2r")
         assert path == "/gt2r/land_ice_segments/latitude"
+
+
+class TestVectorRoundTrip:
+    """Issue #29 phase 6: a vector field written to a real Zarr template reads
+    back through the trailing-dim block, and NaN-padded empty cells are skipped
+    by a NaN-aware reducer."""
+
+    @staticmethod
+    def _vector_cfg():
+        cfg = default_config("atl06")
+        agg = {
+            "coordinates": cfg.aggregation.get("coordinates", {}),
+            "variables": {
+                "count": {"function": "len", "source": "h_li"},
+                "edges": {
+                    "expression": "np.array([np.min(h), np.max(h)])",
+                    "source": "h",
+                    "kind": "vector",
+                    "trailing_shape": 2,
+                    "dtype": "float32",
+                },
+            },
+        }
+        from zagg.config import PipelineConfig
+
+        return PipelineConfig(
+            data_source=cfg.data_source, aggregation=agg, output=cfg.output
+        )
+
+    def test_vector_leaf_to_zarr_to_read(self):
+        pytest.importorskip("pyarrow")
+        from mortie import geo2mort
+
+        cfg = self._vector_cfg()
+        parent_order, child_order = 2, 4
+        grid = HealpixGrid(parent_order, child_order, layout="fullsphere", config=cfg)
+        store = MemoryStore()
+        grid.emit_template(store)
+
+        parent = int(geo2mort(-78.5, -132.0, order=parent_order)[0])
+        children = grid.children(parent)
+        n = len(children)  # 4 ** (child_order - parent_order)
+        assert n == 4 ** (child_order - parent_order)
+
+        # Two populated cells; the rest stay NaN-padded (the empty-cell sentinel).
+        stats = {
+            "count": np.zeros(n, dtype="float32"),
+            "edges": np.full((n, 2), np.nan, dtype="float32"),
+        }
+        stats["count"][0] = 5
+        stats["edges"][0] = [1.0, 9.0]
+        stats["count"][3] = 2
+        stats["edges"][3] = [-2.0, 4.0]
+
+        carrier = _build_output(
+            stats, get_data_vars(cfg), get_agg_fields(cfg), grid, parent, use_arrow=True
+        )
+        # The vector column is carried as a FixedSizeList (issue #29 B').
+        assert carrier.column_names[:2] == ["count", "edges"]
+
+        chunk_idx = grid.block_index(parent)
+        write_dataframe_to_zarr(carrier, store, grid=grid, chunk_idx=chunk_idx)
+
+        group = open_group(store=store, mode="r", path=str(child_order))
+        assert group["edges"].shape == (HEALPIX_BASE_CELLS * 4**child_order, 2)
+        block_start = chunk_idx[0] * n
+        got = group["edges"][block_start : block_start + n]
+
+        # Populated cells round-trip exactly through the trailing-dim selection.
+        np.testing.assert_array_equal(got[0], [1.0, 9.0])
+        np.testing.assert_array_equal(got[3], [-2.0, 4.0])
+        # Empty cells carry the NaN padding sentinel.
+        assert np.all(np.isnan(got[1]))
+        assert np.all(np.isnan(got[2]))
+
+        # A NaN-aware reducer skips the padding: the per-edge mean over cells is
+        # taken only over the two populated rows.
+        reduced = np.nanmean(got, axis=0)
+        np.testing.assert_allclose(reduced, [(1.0 - 2.0) / 2, (9.0 + 4.0) / 2])

From 91cce63c1bb123d1cc243d98573fb0e16c1417fc Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 19:25:43 +0000
Subject: [PATCH 10/11] fold review: enforce single-trailing-chunk invariant +
 name payload axis (#29)

---
 src/zagg/grids/base.py   |  8 +++++++-
 src/zagg/processing.py   | 11 +++++++++++
 tests/test_grids.py      |  3 +--
 tests/test_processing.py | 30 ++++++++++++++++++++++++++++++
 4 files changed, 49 insertions(+), 3 deletions(-)

diff --git a/src/zagg/grids/base.py b/src/zagg/grids/base.py
index 11724fed..579127f0 100644
--- a/src/zagg/grids/base.py
+++ b/src/zagg/grids/base.py
@@ -69,7 +69,13 @@ def vector_array_spec(base, sig, *, base_dims, base_chunk_shape):
     if not trailing:
         return base
     shape = (*base.shape, *trailing)
-    dim_names = (*base_dims, *(f"{base_dims[-1]}_v{i}" for i in range(len(trailing))))
+    # Name trailing payload axes ``vector`` (one dim) or ``vector_0``/``vector_1``
+    # (multi-dim, e.g. a t-digest ``(k, 2)``) — distinct from the spatial axes.
+    if len(trailing) == 1:
+        trailing_names: tuple[str, ...] = ("vector",)
+    else:
+        trailing_names = tuple(f"vector_{i}" for i in range(len(trailing)))
+    dim_names = (*base_dims, *trailing_names)
     chunk_shape = (*base_chunk_shape, *trailing)
     chunk_grid = NamedConfig(name="regular", configuration={"chunk_shape": list(chunk_shape)})
     # Set shape + dimension_names together: ArraySpec validates their ranks
diff --git a/src/zagg/processing.py b/src/zagg/processing.py
index 7c34121e..87634f10 100644
--- a/src/zagg/processing.py
+++ b/src/zagg/processing.py
@@ -284,6 +284,17 @@ def write_dataframe_to_zarr(
                 zarr_format=3,
                 consolidated=False,
             )
+            if trailing:
+                # Enforce the single-trailing-chunk invariant: the target array's
+                # trailing chunk must span the whole payload, or set_block_selection
+                # at block 0 would silently write only part of it (issue #29).
+                target_trailing_chunks = array.chunks[len(grid.chunk_shape) :]
+                if target_trailing_chunks != trailing:
+                    raise ValueError(
+                        f"vector field {name!r}: trailing chunk "
+                        f"{target_trailing_chunks} must equal trailing shape "
+                        f"{trailing} (the payload dim must be one whole chunk)"
+                    )
             array.set_block_selection(block_idx, values)
 
     return store
diff --git a/tests/test_grids.py b/tests/test_grids.py
index d4a0d9ce..b562c439 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -361,8 +361,7 @@ def test_healpix_dimension_names_extend(self):
         g = HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg)
         spec = g.spec()
         names = spec.members["hist"].dimension_names
-        assert names[0] == "cells"
-        assert len(names) == 2  # spatial + one trailing dim
+        assert names == ("cells", "vector")  # spatial + the trailing payload axis
 
     def test_rectilinear_vector_array_has_trailing_dim(self):
         from zagg.grids import RectilinearGrid
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 0308857f..36313cf2 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -973,3 +973,33 @@ def test_vector_leaf_to_zarr_to_read(self):
         # taken only over the two populated rows.
         reduced = np.nanmean(got, axis=0)
         np.testing.assert_allclose(reduced, [(1.0 - 2.0) / 2, (9.0 + 4.0) / 2])
+
+    def test_split_trailing_chunk_rejected(self):
+        """The writer enforces the single-trailing-chunk invariant: if the target
+        array chunks the trailing payload dim, ``set_block_selection`` at block 0
+        would drop the rest, so the write must raise instead (issue #29)."""
+        pa = pytest.importorskip("pyarrow")
+        from zarr import create_array
+
+        class _OneChunkGrid:
+            group_path = "g"
+            chunk_shape = (2,)
+
+        store = MemoryStore()
+        # Trailing dim of width 4 deliberately split into two chunks of 2.
+        create_array(
+            store,
+            name="g/edges",
+            shape=(2, 4),
+            chunks=(2, 2),
+            dtype="float32",
+            fill_value=np.float32("nan"),
+        )
+        edges = pa.FixedSizeListArray.from_arrays(
+            pa.array(np.arange(8.0, dtype="float32")), 4
+        )
+        table = pa.table({"edges": edges})
+        with pytest.raises(ValueError, match="one whole chunk"):
+            write_dataframe_to_zarr(
+                table, store, grid=_OneChunkGrid(), chunk_idx=(0,)
+            )

From ea17405fe3d53c139281b38593c33f8b46d0242e Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Tue, 16 Jun 2026 19:26:47 +0000
Subject: [PATCH 11/11] fold review: clarify round-trip test isolates the
 write/read path (#29)

---
 tests/test_processing.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/test_processing.py b/tests/test_processing.py
index 36313cf2..d4af2392 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -938,6 +938,9 @@ def test_vector_leaf_to_zarr_to_read(self):
         n = len(children)  # 4 ** (child_order - parent_order)
         assert n == 4 ** (child_order - parent_order)
 
+        # This test isolates the carrier->writer->Zarr->reader half of #29, so it
+        # fabricates the per-cell stats blocks directly rather than running the
+        # ``edges`` expression (the stat-eval path is covered by TestVectorOutputs).
         # Two populated cells; the rest stay NaN-padded (the empty-cell sentinel).
         stats = {
             "count": np.zeros(n, dtype="float32"),