#!/usr/bin/env python3 """Measure Xarray fill masking on newly created local synthetic Zarr stores.""" import argparse import hashlib import importlib.metadata import json import math from pathlib import Path import platform import tempfile import warnings import numpy as np import xarray as xr import zarr from zarr.codecs import BytesCodec, Crc32cCodec from xarray.backends.zarr import ZarrStore VALUES = list(range(10, 18)) FILL = 0 def require(condition, message): if not condition: raise RuntimeError(message) def scalar(value): value = float(value) return value if math.isfinite(value) else None def inventory(path): return {p.relative_to(path).as_posix(): hashlib.sha256(p.read_bytes()).hexdigest() for p in sorted(path.rglob("*")) if p.is_file()} def store_bytes(path): return {p.relative_to(path).as_posix(): p.read_bytes() for p in path.rglob("*") if p.is_file()} def create_store(path, version, cf_fill=None, state="intact", values=VALUES): require(version in (2, 3), "Unsupported format") require(state in ("intact", "deleted", "sparse"), "Unsupported state") require(len(values) == 8, "The fixture requires eight values") group = zarr.open_group(store=path, mode="w", zarr_format=version) attrs = {} if cf_fill is None else {"_FillValue": cf_fill} kwargs = {"compressors": None} if version == 2: attrs["_ARRAY_DIMENSIONS"] = ["sample"] else: kwargs = {"dimension_names": ("sample",), "serializer": BytesCodec(endian="little"), "compressors": [Crc32cCodec()]} array = group.create_array( "signal", shape=(8,), chunks=(4,), dtype="int32", fill_value=FILL, attributes=attrs, config={"write_empty_chunks": False}, **kwargs ) if state == "sparse": array[:4] = values[:4] else: array[:] = values before = inventory(path) key = "signal/1" if version == 2 else "signal/c/1" if state == "deleted": require((path / key).is_file(), "Expected initialized chunk is absent before deletion") (path / key).unlink() return sorted(before.keys() - inventory(path).keys()) def read_summary(path, policy=None, mask_and_scale=True, route="public"): # Disable Dask and consolidation so this study isolates fill decoding. with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always") if route == "public": reader = xr.open_zarr( str(path), chunks=None, consolidated=False, use_zarr_fill_value_as_mask=policy, mask_and_scale=mask_and_scale ) else: require(route == "backend_control", "Unknown reader route") store = ZarrStore.open_group( str(path), mode="r", consolidated=False, use_zarr_fill_value_as_mask=policy ) reader = xr.open_dataset( store, engine="store", chunks=None, mask_and_scale=mask_and_scale ) with reader as dataset: data = dataset["signal"].load() result = { "values": [scalar(v) for v in data.values], "dtype": str(data.dtype), "size": data.size, "nonmissing_count": int(data.count().item()), "default_sum": scalar(data.sum().item()), "default_mean": scalar(data.mean().item()), "mean_skipna_false": scalar(data.mean(skipna=False).item()), "null_indices": np.flatnonzero(data.isnull().values).tolist(), } result["warnings"] = [ {"category": type(w.message).__name__, "message": str(w.message)} for w in caught ] return result def expected_values(version, cf_fill, policy, state, mask_and_scale=True, values=VALUES): raw = list(values) if state == "intact" else list(values[:4]) + [FILL] * 4 effective_policy = version == 2 if policy is None else policy sentinel = FILL if effective_policy else cf_fill return [None if mask_and_scale and sentinel is not None and v == sentinel else v for v in raw] def check_summary(result, expected): require(result["values"] == expected, f"Decoded values differ: {result}") valid = [v for v in expected if v is not None] require(result["size"] == len(expected), "Array extent changed") require(result["nonmissing_count"] == len(valid), "Nonmissing count differs") require(result["default_sum"] == sum(valid), "Default sum differs") require(result["default_mean"] == sum(valid) / len(valid), "Default mean differs") strict = None if len(valid) < len(expected) else sum(valid) / len(valid) require(result["mean_skipna_false"] == strict, "Strict mean differs") require(result["null_indices"] == [i for i, v in enumerate(expected) if v is None], "Missing-value positions differ") require(result["warnings"] == [], "Unexpected warning") def one_repetition(root, repetition): rows = [] def record(case, path, version, cf_fill, policy, state, mask_and_scale=True, values=VALUES, missing_keys=None, route="public"): observed = read_summary(path, policy, mask_and_scale, route) documented = expected_values(version, cf_fill, policy, state, mask_and_scale, values) # Both pinned public entrypoints reset the requested option to None. forwarded = None if route == "public" else policy check_summary(observed, expected_values( version, cf_fill, forwarded, state, mask_and_scale, values )) rows.append({ "repetition": repetition, "case": case, "zarr_format": version, "reader_route": route, "documented_expected_values": documented, "matches_documented_policy": observed["values"] == documented, "cf_fill_attribute": cf_fill, "storage_fill_value": FILL, "use_zarr_fill_value_as_mask": policy, "mask_and_scale": mask_and_scale, "state": state, "source_values": list(values), "missing_from_pre_mutation_inventory": missing_keys or [], **observed, }) for version in (2, 3): for cf_fill in (None, 0, -999): paths = {} missing = {} for state in ("intact", "deleted", "sparse"): path = root / f"v{version}-cf{cf_fill}-{state}" paths[state] = path missing[state] = create_store(path, version, cf_fill, state) for policy in (None, True, False): record("fill_policy_matrix", path, version, cf_fill, policy, state, missing_keys=missing[state]) require(store_bytes(paths["deleted"]) == store_bytes(paths["sparse"]), "Sparse and deleted stores must be byte-identical") record("decoding_disabled", paths["deleted"], version, cf_fill, None, "deleted", mask_and_scale=False, missing_keys=missing["deleted"]) for policy in (None, True, False): record("backend_option_control", paths["deleted"], version, cf_fill, policy, "deleted", missing_keys=missing["deleted"], route="backend_control") path = root / f"v{version}-written-zero" values = list(range(8)) missing = create_store(path, version, values=values) for policy in (None, True, False): record("written_zero", path, version, None, policy, "intact", values=values, missing_keys=missing) require(len(rows) == 84, "Unexpected case count") require(sum(not row["matches_documented_policy"] for row in rows) == 10, "Unexpected documented-policy mismatch count") return rows def run(repetitions=3): require(type(repetitions) is int and repetitions > 0, "Repetitions must be a positive integer") require(xr.__version__ in ("2025.10.1", "2026.9.0"), "Review policy forwarding before adding another Xarray version") rows = [] reference = None with tempfile.TemporaryDirectory(prefix="hesela-zarr-xarray-") as directory: for repetition in range(1, repetitions + 1): root = Path(directory) / str(repetition) root.mkdir() current = one_repetition(root, repetition) comparison = [{k: v for k, v in row.items() if k != "repetition"} for row in current] if reference is None: reference = comparison require(comparison == reference, "Deterministic repetitions disagree") rows.extend(current) return { "schema_version": "1.0.0", "evidence_type": "controlled software observations on synthetic inputs", "study": "Xarray Zarr fill masking and aggregate denominators", "environment": { "python": platform.python_version(), "system": platform.system(), "machine": platform.machine(), "packages": {p: importlib.metadata.version(p) for p in ("xarray", "zarr", "numpy", "pandas", "numcodecs", "crc32c")}, }, "case_count": 84, "repetitions": repetitions, "reader_options": {"chunks": None, "consolidated": False}, "observations": rows, "limitations": [ "Two pinned Xarray releases over the same Zarr-Python backend; not independent Zarr implementations.", "The direct backend route is a diagnostic control, not a supported-API migration recommendation.", "Local unsharded synthetic int32 arrays only; no incidence, cloud, power-loss or Dask claim.", "Repeated deterministic runs are a reproducibility check, not independent statistical samples.", "Storage fill is zero; the prior Zarr-only experiment used -999.", "No scaling factors, coordinate indexes, floating-point stored types or nullable fill are tested.", "The expected pre-mutation object inventory is trusted in these controls.", ], } if __name__ == "__main__": parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--repetitions", type=int, default=3) args = parser.parse_args() result = run(args.repetitions) args.output.write_text(json.dumps(result, indent=2, sort_keys=True, allow_nan=False) + "\n") print(f"PASS: {result['case_count']} cases x {result['repetitions']} repetitions")