From 9cc3a2ce27bce9bdf3eaab43284e4da95d856dd0 Mon Sep 17 00:00:00 2001 From: james Date: Wed, 7 Oct 2026 19:49:37 +0100 Subject: [PATCH 01/12] Add a streaming CoverageJSON encoder over the polytope-mars block stream CovjsonStreamEncoder writes the collection opening in begin(), MultiPoint coverages as coordinate and value blocks arrive, and buffers the small PointSeries, VerticalProfile and Trajectory collections until end(). The bytes equal json.dumps of the legacy encoders' output, with NaN written as null. Blocks are read by attribute only: no polytope-feature, polytope-mars or covjson-pydantic imports. --- covjsonkit/stream.py | 453 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 453 insertions(+) create mode 100644 covjsonkit/stream.py diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py new file mode 100644 index 0000000..a80081a --- /dev/null +++ b/covjsonkit/stream.py @@ -0,0 +1,453 @@ +"""Streaming CoverageJSON encoder over a block stream. + +``CovjsonStreamEncoder`` consumes the block IR produced by polytope-mars (``polytope_mars.blocks``) and +returns the CoverageJSON document in pieces: ``begin(header)`` -> collection opening, ``encode(block)`` +per block, ``end()`` -> collection closing. Concatenating the pieces gives exactly the bytes the legacy +path produced with ``json.dumps(encoder.from_polytope*(tree)).encode()``, except that missing values +(NaN) are written as ``null``. + +Blocks are read structurally (attributes only), so this module imports neither polytope-feature nor +polytope-mars: + +* header: ``feature_type``, ``domain_type``, ``time_axis``, ``parameters`` (each with ``id``, + ``shortname``, ``name``, ``unit``, ``description``), ``extra``; +* group (``block.group``): ``index``, ``t``, ``params``, ``levels``, ``n_points``, ``n_bands``, + ``mars_metadata``; +* coordinates block: ``lat``, ``lon`` (float arrays), ``band``; values block: ``param``, ``level``, + ``band``, ``values`` (float array, NaN = missing); group end: only ``group``. + +MultiPoint coverages are written as the blocks arrive (memory independent of the coverage size, apart +from the coordinates of the current coverage when it has several levels). PointSeries, +VerticalProfile and Trajectory collections are small and their legacy layout is point-major across field +groups, so their groups are buffered and written by ``end()``. +""" + +from __future__ import annotations + +import json +import math + +import numpy as np +import orjson + +__all__ = ["CovjsonStreamEncoder", "legacy_referencing", "pointseries_coverages"] + +_CRS = {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"} +_LATLON = ["latitude", "longitude", "levelist"] +_XYZ = ["x", "y", "z"] + + +def _dumps(obj) -> bytes: + """``json.dumps`` bytes (default separators, ASCII escapes) for small structural fragments.""" + return json.dumps(obj).encode("utf-8") + + +def format_floats(values) -> list: + """Each value as ``json.dumps`` writes a Python float, with NaN/inf as ``null``. + + orjson and CPython both print the shortest round-trip representation; they differ only in the + exponent form of numbers below 1e-4 (orjson ``1e-5``/``0.00002`` vs Python ``1e-05``/``2e-05``), + which are re-printed with ``repr``. + """ + arr = np.ascontiguousarray(values, dtype=np.float64) + if arr.size == 0: + return [] + parts = orjson.dumps(arr, option=orjson.OPT_SERIALIZE_NUMPY)[1:-1].split(b",") + with np.errstate(invalid="ignore"): + small = np.flatnonzero((np.abs(arr) < 1e-4) & (arr != 0)) + for i in small: + parts[i] = repr(float(arr[i])).encode("ascii") + return parts + + +def float_list_bytes(values) -> bytes: + """``a, b, c`` (no brackets), as ``json.dumps`` separates list items.""" + return b", ".join(format_floats(values)) + + +def _py_values(values) -> list: + """Floats for ``json.dumps``; NaN -> None.""" + return [None if (v is None or (isinstance(v, float) and not math.isfinite(v))) else v for v in values] + + +def legacy_referencing(header) -> list: + """The ``coordinates`` of the reference system each legacy encoder wrote (preserved quirk).""" + domain, feature, role = header.domain_type, header.feature_type, header.time_axis + if domain == "Trajectory": + return list(_LATLON) if role == "hdate" else ["t", "x", "y", "z"] + if domain == "VerticalProfile": + return list(_LATLON) + if domain == "PointSeries": + if role == "step": + return list(_XYZ) + if role == "month" and feature == "position": + return list(_XYZ) + return list(_LATLON) + # MultiPoint + if role == "hdate": + return list(_LATLON) + if role == "date": + return list(_LATLON) if feature in ("boundingbox", "circle") else list(_XYZ) + if role == "month": + return list(_LATLON) if feature == "circle" else list(_XYZ) + return list(_XYZ) + + +def _parameter(p) -> dict: + return { + "type": "Parameter", + "description": {"en": p.description}, + "unit": {"symbol": p.unit}, + "observedProperty": {"id": p.shortname, "label": {"en": p.name}}, + } + + +def _range(name, values: list) -> dict: + return {"type": "NdArray", "dataType": "float", "shape": [len(values)], "axisNames": [name], "values": values} + + +def _series(groups) -> list: + """Consecutive buffered groups with equal metadata and levels form one coverage series.""" + series = [] + for g in groups: + if series and series[-1][0]["meta"] == g["meta"] and series[-1][0]["levels"] == g["levels"]: + series[-1].append(g) + else: + series.append([g]) + return series + + +def pointseries_coverages(groups, header, shortname) -> list: + """Legacy PointSeries layout (TimeSeries / Position encoders): one coverage per (point, series, level). + + A series is a run of field groups (time instants) sharing their coverage metadata; its ``t`` axis is + the groups' ``t`` values in order. Coverages are point-major (``for point: for series``) unless + ``header.extra["pointseries_order"] == "series_major"``. Position coverages from ``from_polytope`` + write the first level on their levelist axis whatever their level (preserved quirk). + """ + position_quirk = header.feature_type == "position" and header.time_axis == "date" + order = (header.extra or {}).get("pointseries_order", "point_major") + all_series = _series(groups) + + def coverage(series, i, level): + first = series[0] + levels = first["levels"] + if position_quirk and levels: + levelist = [levels[0]] + else: + levelist = [0 if level is None else level] + t = [v for g in series for v in g["t"]] + ranges = {} + for p in header.parameters: + if not any(p.id in g["values_params"] for g in series): + continue + vals = [] + for g in series: + arr = g["values"].get((p.id, level)) + vals.append(float(arr[i]) if arr is not None else None) + ranges[shortname[p.id]] = _range(shortname[p.id], _py_values(vals)) + return { + "mars:metadata": first["meta"], + "type": "Coverage", + "domain": { + "type": "Domain", + "axes": { + "latitude": {"values": [float(first["lat"][i])]}, + "longitude": {"values": [float(first["lon"][i])]}, + "levelist": {"values": levelist}, + "t": {"values": t}, + }, + }, + "ranges": ranges, + } + + out = [] + if order == "series_major": + for series in all_series: + for i in range(len(series[0]["lat"])): + for level in series[0]["levels"] or (None,): + out.append(coverage(series, i, level)) + return out + n_points = len(all_series[0][0]["lat"]) if all_series else 0 + for i in range(n_points): + for series in all_series: + for level in series[0]["levels"] or (None,): + out.append(coverage(series, i, level)) + return out + + +def verticalprofile_coverages(groups, header, shortname) -> list: + """Legacy VerticalProfile layout: one coverage per (point, group), all levels on the levelist axis.""" + out = [] + n_points = len(groups[0]["lat"]) if groups else 0 + for i in range(n_points): + for g in groups: + levels = list(g["levels"]) if g["levels"] else [0] + keys = list(g["levels"]) if g["levels"] else [None] + ranges = {} + for p in header.parameters: + if p.id not in g["values_params"]: + continue + vals = [float(g["values"][(p.id, k)][i]) for k in keys] + ranges[shortname[p.id]] = { + "type": "NdArray", + "dataType": "float", + "shape": [len(vals)], + "axisNames": ["levelist"], + "values": _py_values(vals), + } + out.append( + { + "mars:metadata": g["meta"], + "type": "Coverage", + "domain": { + "type": "Domain", + "axes": { + "latitude": {"values": [float(g["lat"][i])]}, + "longitude": {"values": [float(g["lon"][i])]}, + "levelist": {"values": levels}, + "t": {"values": list(g["t"])}, + }, + }, + "ranges": ranges, + } + ) + return out + + +def trajectory_coverages(groups, header, shortname, coordinates) -> list: + """Legacy Trajectory layout: one coverage per series; composite ``[t, lat, lon, level]``.""" + out = [] + for series in _series(groups): + composite = [] + values = {p.id: [] for p in header.parameters} + for g in series: + t = g["t"][0] if g["t"] else 0 + keys = list(g["levels"]) if g["levels"] else [None] + for k in keys: + level = 0 if k is None else k + composite.extend([t, float(la), float(lo), level] for la, lo in zip(g["lat"], g["lon"])) + for p in header.parameters: + arr = g["values"].get((p.id, k)) + values[p.id].extend(arr.tolist() if arr is not None else [None] * len(g["lat"])) + ranges = { + shortname[p.id]: _range(shortname[p.id], _py_values(values[p.id])) + for p in header.parameters + if any(p.id in g["values_params"] for g in series) + } + out.append( + { + "mars:metadata": series[0]["meta"], + "type": "Coverage", + "domain": { + "type": "Domain", + "axes": {"composite": {"dataType": "tuple", "coordinates": coordinates, "values": composite}}, + }, + "ranges": ranges, + } + ) + return out + + +class _OpenCoverage: + """State of the MultiPoint coverage being written.""" + + __slots__ = ("group", "n_tuples", "coords", "composite_closed", "param", "n_values") + + def __init__(self, group): + self.group = group + self.n_tuples = 0 + #: formatted (lat, lon) of every band, kept only when the coverage has more than one level + self.coords: list = [] + self.composite_closed = False + self.param = None + self.n_values = 0 + + +class CovjsonStreamEncoder: + """CoverageJSON encoder for the polytope-mars block stream (one instance per request).""" + + content_type = "application/prs.coverage+json" + file_extension = "covjson" + + def __init__(self, config=None): + self.config = dict(config or {}) + self._header = None + self.n_coverages = 0 + self._streaming = True + self._buffered: list = [] + self._current: _OpenCoverage | None = None + + @property + def header(self): + if self._header is None: + raise RuntimeError("CovjsonStreamEncoder.begin() must be called first") + return self._header + + def _open(self) -> _OpenCoverage: + if self._current is None: + raise RuntimeError("values block before the coordinates of its group") + return self._current + + # -- protocol ---------------------------------------------------------------------------------------- + + def begin(self, header) -> bytes: + self._header = header + self.n_coverages = 0 + self._shortname = {p.id: p.shortname for p in header.parameters} + self._coordinates = legacy_referencing(header) + self._streaming = header.domain_type == "MultiPoint" or ( + header.domain_type == "Trajectory" and header.time_axis == "hdate" + ) + self._with_t = header.domain_type == "MultiPoint" + self._buffered = [] + self._current = None + domain = "MultiPoint" if header.domain_type in ("shapefile",) else header.domain_type + return b'{"type": "CoverageCollection", "domainType": ' + _dumps(domain) + b', "coverages": [' + + def encode(self, block) -> bytes: + if hasattr(block, "lat"): + return self._coords(block) + if hasattr(block, "values") and hasattr(block, "param"): + return self._values(block) + return self._group_end(block) + + def end(self) -> bytes: + out = [] + if not self._streaming: + covs = self._buffered_coverages() + for i, cov in enumerate(covs): + out.append((b", " if self.n_coverages or i else b"") + _dumps(cov)) + self.n_coverages += len(covs) + self._buffered = [] + parameters = {p.shortname: _parameter(p) for p in self.header.parameters} + referencing = [{"coordinates": self._coordinates, "system": dict(_CRS)}] + out.append(b'], "referencing": ' + _dumps(referencing) + b', "parameters": ' + _dumps(parameters) + b"}") + return b"".join(out) + + # -- MultiPoint: stream ------------------------------------------------------------------------------- + + def _coords(self, block) -> bytes: + g = block.group + if not self._streaming: + cur = self._buffer_group(g) + cur["lat"].append(np.asarray(block.lat, dtype=np.float64)) + cur["lon"].append(np.asarray(block.lon, dtype=np.float64)) + return b"" + out = [] + if self._current is None or self._current.group is not g: + out.append(self._open_coverage(g)) + cur = self._open() + levels = list(g.levels) if g.levels else [0] + lat = format_floats(block.lat) + lon = format_floats(block.lon) + suffix = b", " + _dumps(levels[0]) + b"]" + if lat: + out.append(self._tuples(cur, lat, lon, suffix)) + if len(levels) > 1: + cur.coords.append((lat, lon)) + return b"".join(out) + + @staticmethod + def _tuples(cur: _OpenCoverage, lat: list, lon: list, suffix: bytes) -> bytes: + body = b", ".join([b"[" + a + b", " + o + suffix for a, o in zip(lat, lon)]) + sep = b", " if cur.n_tuples else b"" + cur.n_tuples += len(lat) + return sep + body + + def _open_coverage(self, g) -> bytes: + self._current = _OpenCoverage(g) + head = (b", " if self.n_coverages else b"") + b'{"mars:metadata": ' + _dumps(g.mars_metadata) + head += b', "type": "Coverage", "domain": {"type": "Domain", "axes": {' + if self._with_t: + head += b'"t": {"values": ' + _dumps(list(g.t)) + b"}, " + head += b'"composite": {"dataType": "tuple", "coordinates": ' + _dumps(self._coordinates) + b', "values": [' + self.n_coverages += 1 + return head + + def _close_composite(self) -> bytes: + cur = self._open() + if cur.composite_closed: + return b"" + cur.composite_closed = True + out = [] + levels = list(cur.group.levels) + for level in levels[1:]: + suffix = b", " + _dumps(level) + b"]" + for lat, lon in cur.coords: + if lat: + out.append(self._tuples(cur, lat, lon, suffix)) + cur.coords = [] + out.append(b']}}}, "ranges": {') + return b"".join(out) + + def _values(self, block) -> bytes: + g = block.group + if not self._streaming: + cur = self._buffer_group(g) + cur["values"].setdefault((block.param, block.level), []).append(np.asarray(block.values, dtype=np.float64)) + return b"" + cur = self._open() + out = [self._close_composite()] + if cur.param != block.param: + name = self._shortname.get(block.param, block.param) + n = g.n_points * max(1, len(g.levels)) + out.append(b"]}, " if cur.param is not None else b"") + out.append(_dumps(name) + b': {"type": "NdArray", "dataType": "float", "shape": [' + str(n).encode()) + out.append(b'], "axisNames": [' + _dumps(str(name)) + b'], "values": [') + cur.param = block.param + cur.n_values = 0 + body = float_list_bytes(block.values) + if body: + out.append((b", " if cur.n_values else b"") + body) + cur.n_values += len(block.values) + return b"".join(out) + + def _group_end(self, block) -> bytes: + if not self._streaming: + cur = self._buffered[-1] if self._buffered else None + if cur is not None and cur["group"] is block.group: + cur["lat"] = np.concatenate(cur["lat"]) if cur["lat"] else np.empty(0) + cur["lon"] = np.concatenate(cur["lon"]) if cur["lon"] else np.empty(0) + cur["values"] = {k: np.concatenate(v) for k, v in cur["values"].items()} + cur["done"] = True + return b"" + cur = self._current + if cur is None or cur.group is not block.group: + return b"" + out = [self._close_composite()] + out.append(b"]}}}" if cur.param is not None else b"}}") + self._current = None + return b"".join(out) + + # -- buffered domains --------------------------------------------------------------------------------- + + def _buffer_group(self, g) -> dict: + if not self._buffered or self._buffered[-1]["group"] is not g: + self._buffered.append( + { + "group": g, + "meta": g.mars_metadata, + "levels": tuple(g.levels), + "t": tuple(g.t), + "values_params": set(g.params), + "lat": [], + "lon": [], + "values": {}, + "done": False, + } + ) + return self._buffered[-1] + + def _buffered_coverages(self) -> list: + groups = [g for g in self._buffered if g["done"]] + if not groups: + return [] + domain = self.header.domain_type + if domain == "PointSeries": + return pointseries_coverages(groups, self.header, self._shortname) + if domain == "VerticalProfile": + return verticalprofile_coverages(groups, self.header, self._shortname) + if domain == "Trajectory": + return trajectory_coverages(groups, self.header, self._shortname, self._coordinates) + raise ValueError(f"Unsupported domain type {domain!r}") From 4b4e0daaf69f6c6e944921f181e554843dbdb14b Mon Sep 17 00:00:00 2001 From: james Date: Wed, 7 Oct 2026 19:49:37 +0100 Subject: [PATCH 02/12] Test the streaming encoder on synthetic blocks Checks legacy-identical bytes for MultiPoint, identical output for any band split, null for NaN, omitted ranges for absent params, float formatting against json.dumps, the referencing quirks, the point-feature layouts and that the module imports no polytope packages. --- tests/test_stream_encoder.py | 302 +++++++++++++++++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 tests/test_stream_encoder.py diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py new file mode 100644 index 0000000..1ea84e3 --- /dev/null +++ b/tests/test_stream_encoder.py @@ -0,0 +1,302 @@ +"""CovjsonStreamEncoder on synthetic blocks (no polytope needed). + +The blocks are local stand-ins with the attributes polytope-mars' block IR exposes; the encoder only +reads attributes. Expected bytes are ``json.dumps`` of the document the legacy encoders built. +""" + +import json +import subprocess +import sys +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +import pytest + +from covjsonkit.stream import CovjsonStreamEncoder, format_floats, legacy_referencing + + +@dataclass(frozen=True) +class Param: + id: str + shortname: str + name: str + unit: str + description: str = "" + + +@dataclass(frozen=True) +class Header: + feature_type: str + domain_type: str + time_axis: str + parameters: tuple + mars_metadata: dict = field(default_factory=dict) + extra: dict = field(default_factory=dict) + + +@dataclass(frozen=True) +class Group: + index: int + path: dict + t: tuple + params: tuple + levels: tuple + n_points: int + n_bands: int + mars_metadata: dict + + +@dataclass(frozen=True) +class Coords: + group: Any + band: int + offset: int + lat: np.ndarray + lon: np.ndarray + + +@dataclass(frozen=True) +class Values: + group: Any + param: str + level: Any + band: int + offset: int + values: np.ndarray + + +@dataclass(frozen=True) +class End: + group: Any + + +T2 = Param("167", "2t", "2 metre temperature", "K", "Température à 2 m") +U10 = Param("165", "10u", "10 metre U wind component", "m s**-1") +CRS = {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"} + + +def _parameters(params): + return { + p.shortname: { + "type": "Parameter", + "description": {"en": p.description}, + "unit": {"symbol": p.unit}, + "observedProperty": {"id": p.shortname, "label": {"en": p.name}}, + } + for p in params + } + + +def _nan_to_none(values): + return [None if v != v else v for v in np.asarray(values, dtype=np.float64).tolist()] + + +def multipoint_blocks(fields, lat, lon, levels=(), meta=None, t=("2024-01-01T00:00:00Z",), n_bands=1, index=0): + """Blocks of one MultiPoint group; ``fields`` = {(param, level): values}.""" + params = tuple(dict.fromkeys(p for p, _ in fields)) + bounds = np.array_split(np.arange(len(lat)), n_bands) + g = Group(index, {}, t, params, tuple(levels), len(lat), n_bands, meta or {"class": "od", "number": 0}) + starts = [idx[0].item() if len(idx) else 0 for idx in bounds] + blocks: list[Any] = [Coords(g, b, starts[b], lat[idx], lon[idx]) for b, idx in enumerate(bounds)] + for p in params: + for lev in levels or (None,): + for b, idx in enumerate(bounds): + blocks.append(Values(g, p, lev, b, starts[b], fields[(p, lev)][idx])) + blocks.append(End(g)) + return blocks + + +def encode(header, blocks) -> bytes: + enc = CovjsonStreamEncoder({"param_db": "ecmwf"}) + out = [enc.begin(header)] + out += [enc.encode(b) for b in blocks] + out.append(enc.end()) + return b"".join(out) + + +def legacy_multipoint(header, groups): + """What BoundingBox.from_polytope + json.dumps produced for the same data.""" + coords = legacy_referencing(header) + covs = [] + for fields, lat, lon, levels, meta, t in groups: + composite = [[float(a), float(o), lev] for lev in (levels or (0,)) for a, o in zip(lat, lon)] + ranges = {} + params = tuple(dict.fromkeys(p for p, _ in fields)) + short = {p.id: p.shortname for p in header.parameters} + for p in params: + vals = [v for lev in (levels or (None,)) for v in _nan_to_none(fields[(p, lev)])] + ranges[short[p]] = { + "type": "NdArray", + "dataType": "float", + "shape": [len(vals)], + "axisNames": [short[p]], + "values": vals, + } + covs.append( + { + "mars:metadata": meta, + "type": "Coverage", + "domain": { + "type": "Domain", + "axes": { + "t": {"values": list(t)}, + "composite": {"dataType": "tuple", "coordinates": coords, "values": composite}, + }, + }, + "ranges": ranges, + } + ) + doc = { + "type": "CoverageCollection", + "domainType": "MultiPoint", + "coverages": covs, + "referencing": [{"coordinates": coords, "system": CRS}], + "parameters": _parameters(header.parameters), + } + return json.dumps(doc).encode() + + +def _data(n, seed=0): + rng = np.random.default_rng(seed) + lat = np.round(rng.uniform(-90, 90, n), 12) + lon = np.round(rng.uniform(0, 360, n), 12) + return lat, lon + + +BBOX = Header("boundingbox", "MultiPoint", "date", (U10, T2)) + + +def test_multipoint_matches_legacy_json_dumps(): + lat, lon = _data(7) + fields = {("165", None): np.linspace(1, 7, 7) * 1.1e10 + 0.123, ("167", None): np.arange(7) * 0.5} + meta = {"class": "od", "Forecast date": "2024-01-01T00:00:00Z", "number": 0, "step": 0} + out = encode(BBOX, multipoint_blocks(fields, lat, lon, meta=meta)) + expected = legacy_multipoint(BBOX, [(fields, lat, lon, (), meta, ("2024-01-01T00:00:00Z",))]) + assert out == expected + json.loads(out) + + +@pytest.mark.parametrize("n_bands", [1, 2, 3, 13]) +def test_band_size_invariance(n_bands): + lat, lon = _data(13, 1) + levels = ("500", "850") + fields = {(p, lev): np.random.default_rng(2).normal(size=13) * 1e3 for p in ("165", "167") for lev in levels} + ref = encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels, n_bands=1)) + assert encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels, n_bands=n_bands)) == ref + group = (fields, lat, lon, levels, {"class": "od", "number": 0}, ("2024-01-01T00:00:00Z",)) + expected = legacy_multipoint(BBOX, [group]) + assert ref == expected + + +def test_nan_is_null_and_missing_param_is_omitted(): + lat, lon = _data(4, 3) + vals = np.array([1.0, np.nan, 3.0, np.nan]) + out = encode(BBOX, multipoint_blocks({("167", None): vals}, lat, lon)) + doc = json.loads(out) + assert b"NaN" not in out + ranges = doc["coverages"][0]["ranges"] + assert list(ranges) == ["2t"] # 10u not in group.params -> no range + assert ranges["2t"]["values"] == [1.0, None, 3.0, None] + assert list(doc["parameters"]) == ["10u", "2t"] # header parameters are always listed + + +def test_several_groups_and_empty_collection(): + lat, lon = _data(3, 4) + blocks = multipoint_blocks({("167", None): np.ones(3)}, lat, lon, index=0) + blocks += multipoint_blocks({("167", None): np.zeros(3)}, lat, lon, index=1, t=("2024-01-02T00:00:00Z",)) + doc = json.loads(encode(BBOX, blocks)) + assert [c["domain"]["axes"]["t"]["values"] for c in doc["coverages"]] == [ + ["2024-01-01T00:00:00Z"], + ["2024-01-02T00:00:00Z"], + ] + empty = json.loads(encode(BBOX, [])) + assert empty["coverages"] == [] and list(empty) == ["type", "domainType", "coverages", "referencing", "parameters"] + + +def test_float_formatting_matches_json_dumps(): + rng = np.random.default_rng(5) + vals = np.concatenate( + [ + rng.normal(size=5000) * 10.0 ** rng.integers(-12, 20, 5000), + [0.0, -0.0, 1.0, 1e16, 1e-4, 9.99e-5, 2e-5, 1e-7, -3e-9, 123456789.123, 5e-324, 1.7976931348623157e308], + ] + ) + assert b", ".join(format_floats(vals)) == json.dumps(vals.tolist()).encode()[1:-1] + assert format_floats(np.array([np.nan, np.inf])) == [b"null", b"null"] + + +@pytest.mark.parametrize( + "feature,role,coords", + [ + ("boundingbox", "date", ["latitude", "longitude", "levelist"]), + ("polygon", "date", ["x", "y", "z"]), + ("polygon", "step", ["x", "y", "z"]), + ("polygon", "hdate", ["latitude", "longitude", "levelist"]), + ("boundingbox", "month", ["x", "y", "z"]), + ("circle", "month", ["latitude", "longitude", "levelist"]), + ], +) +def test_multipoint_referencing_quirks(feature, role, coords): + referencing = legacy_referencing(Header(feature, "MultiPoint", role, ())) + assert referencing == coords + + +def _point_group(index, t, meta, values, levels=(), lat=(51.5, -33.9), lon=(0.1, 18.4)): + params = tuple(dict.fromkeys(p for p, _ in values)) + g = Group(index, {}, t, params, tuple(levels), len(lat), 1, meta) + blocks: list[Any] = [Coords(g, 0, 0, np.array(lat), np.array(lon))] + for (p, lev), v in values.items(): + blocks.append(Values(g, p, lev, 0, 0, np.asarray(v, dtype=float))) + return blocks + [End(g)] + + +def test_pointseries_layout_is_point_major_over_consecutive_groups(): + header = Header("timeseries", "PointSeries", "date", (T2,)) + meta = {"class": "od", "number": 0, "levelist": 0} + blocks = [] + for i, t in enumerate(["2024-01-01T00:00:00Z", "2024-01-01T01:00:00Z", "2024-01-01T02:00:00Z"]): + blocks += _point_group(i, (t,), meta, {("167", None): [i, 10 + i]}) + doc = json.loads(encode(header, blocks)) + assert len(doc["coverages"]) == 2 + first, second = doc["coverages"] + assert first["domain"]["axes"]["t"]["values"] == [ + "2024-01-01T00:00:00Z", + "2024-01-01T01:00:00Z", + "2024-01-01T02:00:00Z", + ] + assert first["ranges"]["2t"]["values"] == [0.0, 1.0, 2.0] + assert second["ranges"]["2t"]["values"] == [10.0, 11.0, 12.0] + assert second["domain"]["axes"]["latitude"]["values"] == [-33.9] + assert doc["referencing"][0]["coordinates"] == ["latitude", "longitude", "levelist"] + + +def test_verticalprofile_layout(): + header = Header("verticalprofile", "VerticalProfile", "date", (T2,)) + blocks = _point_group( + 0, ("2024-01-01T00:00:00Z",), {"number": 0}, {("167", 500): [1, 2], ("167", 850): [3, 4]}, levels=(500, 850) + ) + doc = json.loads(encode(header, blocks)) + assert [c["ranges"]["2t"]["values"] for c in doc["coverages"]] == [[1.0, 3.0], [2.0, 4.0]] + assert doc["coverages"][0]["domain"]["axes"]["levelist"]["values"] == [500, 850] + assert doc["coverages"][0]["ranges"]["2t"]["axisNames"] == ["levelist"] + + +def test_trajectory_layout(): + header = Header("trajectory", "Trajectory", "date", (T2,)) + blocks = _point_group(0, (0,), {"number": 0}, {("167", None): [1, 2]}) + doc = json.loads(encode(header, blocks)) + (cov,) = doc["coverages"] + assert cov["domain"]["axes"]["composite"]["values"] == [[0, 51.5, 0.1, 0], [0, -33.9, 18.4, 0]] + assert cov["domain"]["axes"]["composite"]["coordinates"] == ["t", "x", "y", "z"] + + +def test_no_polytope_imports(): + code = ( + "import sys, covjsonkit.stream; " + "roots = ('polytope_feature', 'polytope_mars', 'covjson_pydantic'); " + "bad = [m for m in sys.modules if m.split('.')[0] in roots]; " + "print(bad); sys.exit(1 if bad else 0)" + ) + proc = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True) + assert proc.returncode == 0, proc.stdout + proc.stderr From c7514352c922226fb034e4f59ac944fd7c969145 Mon Sep 17 00:00:00 2001 From: james Date: Wed, 7 Oct 2026 19:57:36 +0100 Subject: [PATCH 03/12] Format coordinates and values without a Python object per point Composite tuples and value lists are formatted by orjson over the whole array and adjusted to json.dumps separators with bytes.replace; arrays holding values below 1e-4 (printed differently by orjson) keep the per-value path. Halves peak memory of a 634k-point coverage. --- covjsonkit/stream.py | 49 ++++++++++++++++++++++++------------ tests/test_stream_encoder.py | 23 ++++++++++++++++- 2 files changed, 55 insertions(+), 17 deletions(-) diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py index a80081a..5d3084f 100644 --- a/covjsonkit/stream.py +++ b/covjsonkit/stream.py @@ -42,6 +42,12 @@ def _dumps(obj) -> bytes: return json.dumps(obj).encode("utf-8") +def _needs_repr(arr: np.ndarray) -> np.ndarray: + """Mask of the values orjson prints differently from ``json.dumps`` (0 < |x| < 1e-4).""" + with np.errstate(invalid="ignore"): + return (np.abs(arr) < 1e-4) & (arr != 0) + + def format_floats(values) -> list: """Each value as ``json.dumps`` writes a Python float, with NaN/inf as ``null``. @@ -53,16 +59,31 @@ def format_floats(values) -> list: if arr.size == 0: return [] parts = orjson.dumps(arr, option=orjson.OPT_SERIALIZE_NUMPY)[1:-1].split(b",") - with np.errstate(invalid="ignore"): - small = np.flatnonzero((np.abs(arr) < 1e-4) & (arr != 0)) - for i in small: + for i in np.flatnonzero(_needs_repr(arr)): parts[i] = repr(float(arr[i])).encode("ascii") return parts def float_list_bytes(values) -> bytes: """``a, b, c`` (no brackets), as ``json.dumps`` separates list items.""" - return b", ".join(format_floats(values)) + arr = np.ascontiguousarray(values, dtype=np.float64) + if not _needs_repr(arr).any(): + # fast path without one Python object per value + return orjson.dumps(arr, option=orjson.OPT_SERIALIZE_NUMPY)[1:-1].replace(b",", b", ") + return b", ".join(format_floats(arr)) + + +def composite_tuples(lat, lon, level) -> bytes: + """``[lat, lon, level], ...`` (no outer brackets) as ``json.dumps`` writes the composite values.""" + lat = np.ascontiguousarray(lat, dtype=np.float64) + lon = np.ascontiguousarray(lon, dtype=np.float64) + if lat.size == 0: + return b"" + suffix = b", " + _dumps(level) + b"]" + if not (_needs_repr(lat).any() or _needs_repr(lon).any()): + pairs = orjson.dumps(np.column_stack((lat, lon)), option=orjson.OPT_SERIALIZE_NUMPY)[1:-1] + return pairs.replace(b",", b", ").replace(b"]", suffix) + return b", ".join([b"[" + a + b", " + o + suffix for a, o in zip(format_floats(lat), format_floats(lon))]) def _py_values(values) -> list: @@ -257,7 +278,7 @@ class _OpenCoverage: def __init__(self, group): self.group = group self.n_tuples = 0 - #: formatted (lat, lon) of every band, kept only when the coverage has more than one level + #: (lat, lon) arrays of every band, kept only when the coverage has more than one level self.coords: list = [] self.composite_closed = False self.param = None @@ -339,18 +360,16 @@ def _coords(self, block) -> bytes: out.append(self._open_coverage(g)) cur = self._open() levels = list(g.levels) if g.levels else [0] - lat = format_floats(block.lat) - lon = format_floats(block.lon) - suffix = b", " + _dumps(levels[0]) + b"]" - if lat: - out.append(self._tuples(cur, lat, lon, suffix)) + out.append(self._tuples(cur, block.lat, block.lon, levels[0])) if len(levels) > 1: - cur.coords.append((lat, lon)) + cur.coords.append((np.array(block.lat, dtype=np.float64), np.array(block.lon, dtype=np.float64))) return b"".join(out) @staticmethod - def _tuples(cur: _OpenCoverage, lat: list, lon: list, suffix: bytes) -> bytes: - body = b", ".join([b"[" + a + b", " + o + suffix for a, o in zip(lat, lon)]) + def _tuples(cur: _OpenCoverage, lat, lon, level) -> bytes: + body = composite_tuples(lat, lon, level) + if not body: + return b"" sep = b", " if cur.n_tuples else b"" cur.n_tuples += len(lat) return sep + body @@ -373,10 +392,8 @@ def _close_composite(self) -> bytes: out = [] levels = list(cur.group.levels) for level in levels[1:]: - suffix = b", " + _dumps(level) + b"]" for lat, lon in cur.coords: - if lat: - out.append(self._tuples(cur, lat, lon, suffix)) + out.append(self._tuples(cur, lat, lon, level)) cur.coords = [] out.append(b']}}}, "ranges": {') return b"".join(out) diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py index 1ea84e3..dcccd73 100644 --- a/tests/test_stream_encoder.py +++ b/tests/test_stream_encoder.py @@ -13,7 +13,13 @@ import numpy as np import pytest -from covjsonkit.stream import CovjsonStreamEncoder, format_floats, legacy_referencing +from covjsonkit.stream import ( + CovjsonStreamEncoder, + composite_tuples, + float_list_bytes, + format_floats, + legacy_referencing, +) @dataclass(frozen=True) @@ -226,6 +232,21 @@ def test_float_formatting_matches_json_dumps(): assert format_floats(np.array([np.nan, np.inf])) == [b"null", b"null"] +@pytest.mark.parametrize("small", [False, True], ids=["fast-path", "repr-fallback"]) +@pytest.mark.parametrize("level", [0, "500", 850]) +def test_composite_and_value_bytes_match_json_dumps(small, level): + lat, lon = _data(2000, 6) + lon[3] = -0.0 + if small: + lon[7] = 3e-5 # forces the per-value path + pairs = [[a, o, level] for a, o in zip(lat.tolist(), lon.tolist())] + assert composite_tuples(lat, lon, level) == json.dumps(pairs).encode()[1:-1] + values = np.concatenate([lat * 1e9, [np.nan, 2e-7 if small else 2.0]]) + expected = json.dumps([None if v != v else v for v in values.tolist()]).encode()[1:-1] + assert float_list_bytes(values) == expected + assert composite_tuples(np.empty(0), np.empty(0), level) == b"" + + @pytest.mark.parametrize( "feature,role,coords", [ From e2bb09550cf94d08c2d9d70b355ca0335670e3eb Mon Sep 17 00:00:00 2001 From: james Date: Wed, 7 Oct 2026 20:08:04 +0100 Subject: [PATCH 04/12] Document the streaming encoder in CHANGES.md --- CHANGES.md | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 CHANGES.md diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 index 0000000..d79bf92 --- /dev/null +++ b/CHANGES.md @@ -0,0 +1,50 @@ +# Changes on `feat/streaming-encoder` + +Per-branch note for the PR description. The legacy encoder API (`Covjsonkit().encode(...)`, +`from_polytope*`, decoders, `param_db`) is unchanged; all existing tests pass. + +## New: `covjsonkit.stream.CovjsonStreamEncoder` + +A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars.blocks`, DESIGN §3): +`begin(header) -> bytes`, `encode(block) -> bytes`, `end() -> bytes`, `content_type = +"application/prs.coverage+json"`, `file_extension = "covjson"`. polytope-mars selects it for `format: covjson` +(lazy import, `polytope_mars.encoders`). + +- Blocks are read by attribute only; the module imports neither polytope-feature, polytope-mars nor + covjson-pydantic (tested). Parameter metadata comes resolved in the header (`ParamInfo`), so the encoder does + no `param_db` lookups. +- **MultiPoint** coverages are written as blocks arrive: coverage opening + `mars:metadata` + `t` on the first + coordinates block, `[lat, lon, level]` tuples (levels outer, points inner), each range opened on its first + values block and streamed band by band, closed on `GroupEnd`. Memory per coverage: one block, plus the + coordinates when the coverage has several levels (they are written once per level). +- **PointSeries, VerticalProfile, Trajectory** collections are buffered and written by `end()` (point features + are small and the legacy layout is point-major across field groups). The PointSeries layout is the single + function `pointseries_coverages` (to port covjsonkit `develop` #137 there); consecutive groups with equal + `mars_metadata` and `levels` form one series (one coverage per point and level). +- `referencing` and the composite `coordinates` reproduce what each legacy encoder method wrote + (`legacy_referencing`: `latitude/longitude/levelist`, `x/y/z` or `t/x/y/z` by feature and time-axis role). +- **Bytes:** identical to `json.dumps(legacy_covjson).encode()` (separators `", "`/`": "`, ASCII escapes, + key order `type, domainType, coverages, referencing, parameters`). Fragments with strings (metadata, + parameters, referencing) use `json.dumps`; numeric arrays use `orjson` with `OPT_SERIALIZE_NUMPY` over the + whole array, adjusted to `json.dumps` separators with `bytes.replace`. orjson and CPython agree on the + shortest round-trip digits; they differ only in how numbers with 0 < |x| < 1e-4 are spelled + (`1e-5` vs `1e-05`, `0.00002` vs `2e-05`), and arrays containing such values are formatted per value with + `repr`. `-0.0`, integer-valued floats and `1e+16`-style exponents already match (tested on 5,000 random + magnitudes from 1e-12 to 1e20). + +## Output differences against the legacy encoders + +- `NaN`/`inf` values are written as `null` (legacy `json.dumps` wrote the invalid JSON token `NaN`). +- A param absent from a field group (gribjump had no message) has no range in that coverage; a group with no + param at all produces no coverage (decided by polytope-mars, DESIGN §2.5). The collection's `parameters` + always lists every parameter of the header. +- Everything else (including the legacy quirks: space-separated datetimes on the `_step` path, the + `referencing` variants, int `realization`, efcl without `Forecast date`) is byte-identical; see the + polytope-mars golden corpus and its `CHANGES.md`. + +## Tests + +`tests/test_stream_encoder.py` (24 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, +band-size invariance (1, 2, 3 and 13 bands, with levels), `null` for NaN and omitted ranges, float and +composite formatting against `json.dumps` (fast and per-value paths), referencing quirks, PointSeries, +VerticalProfile and Trajectory layouts, no polytope imports. From c83888de96c6f6de7ed2a4c39bdb79feaed4ca09 Mon Sep 17 00:00:00 2001 From: james Date: Thu, 8 Oct 2026 00:43:29 +0100 Subject: [PATCH 05/12] Hand each block to the consumer in bounded fragments A whole-world coverage is one block of up to 26.3M points: its coordinates are ~800 MB of JSON text and its values ~470 MB, which the encoder used to return as a single bytes object (and the fe-worker keeps up to four of those in flight). encode_iter(block) now yields the block a slice at a time, each fragment at most max_fragment_bytes (constructor argument or config key, default 8 MiB): the (n, 3) composite array and the per-value float formatting are done per slice, so the text of a whole block never exists at once. encode(block) stays as b"".join(encode_iter(block)) for buffered callers, and the bytes are unchanged - separators, float spelling and the 1e-5-style per-value fallback are decided slice by slice but produce the same output as one pass over the array (tested at 64 B, 1 KiB and 8 MiB fragment limits against the single-fragment and legacy json.dumps bytes). Also reads numpy scalars with .item() instead of float() in the buffered point-feature layouts. --- covjsonkit/stream.py | 172 ++++++++++++++++++++++++----------- tests/test_stream_encoder.py | 74 +++++++++++++++ 2 files changed, 195 insertions(+), 51 deletions(-) diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py index 5d3084f..1fd58c0 100644 --- a/covjsonkit/stream.py +++ b/covjsonkit/stream.py @@ -1,10 +1,21 @@ """Streaming CoverageJSON encoder over a block stream. ``CovjsonStreamEncoder`` consumes the block IR produced by polytope-mars (``polytope_mars.blocks``) and -returns the CoverageJSON document in pieces: ``begin(header)`` -> collection opening, ``encode(block)`` -per block, ``end()`` -> collection closing. Concatenating the pieces gives exactly the bytes the legacy -path produced with ``json.dumps(encoder.from_polytope*(tree)).encode()``, except that missing values -(NaN) are written as ``null``. +returns the CoverageJSON document in pieces: ``begin(header)`` -> collection opening, +``encode_iter(block)`` -> the fragments of one block, ``end()`` -> collection closing. Concatenating the +pieces gives exactly the bytes the legacy path produced with +``json.dumps(encoder.from_polytope*(tree)).encode()``, except that missing values (NaN) are written as +``null``. + +``encode_iter`` yields one fragment per slice of the block's arrays, each at most +``max_fragment_bytes`` (constructor argument, default :data:`DEFAULT_MAX_FRAGMENT_BYTES` = 8 MiB), so +neither the encoder nor its consumer ever holds the text of a whole block: a whole-world O2560 +coordinates block is ~800 MB of JSON but is handed over 8 MiB at a time. The fragments of a block must +be consumed completely and in order before the next block is encoded (the encoder's state advances as +they are produced). ``encode(block)`` is the buffered equivalent, ``b"".join(encode_iter(block))``, kept +for callers that want one ``bytes`` per block. The structural fragments (a coverage's metadata and ``t`` +axis, the collection's ``referencing`` and ``parameters``) are written whole: they are a few hundred +bytes, and only point-feature collections, which are small by construction, are buffered to ``end()``. Blocks are read structurally (attributes only), so this module imports neither polytope-feature nor polytope-mars: @@ -26,11 +37,26 @@ import json import math +from typing import Iterator import numpy as np import orjson -__all__ = ["CovjsonStreamEncoder", "legacy_referencing", "pointseries_coverages"] +__all__ = [ + "CovjsonStreamEncoder", + "DEFAULT_MAX_FRAGMENT_BYTES", + "legacy_referencing", + "pointseries_coverages", +] + +#: Default upper bound on the size of one fragment of :meth:`CovjsonStreamEncoder.encode_iter`. +DEFAULT_MAX_FRAGMENT_BYTES = 8 * 1024 * 1024 + +#: Longest text one float64 can take (``-1.2345678901234567e-308``), used to size array slices. +_FLOAT_CHARS = 24 + +#: Text built in one serialiser call; a fragment is assembled from pieces of this size. +_WORK_BYTES = 256 * 1024 _CRS = {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"} _LATLON = ["latitude", "longitude", "levelist"] @@ -48,6 +74,33 @@ def _needs_repr(arr: np.ndarray) -> np.ndarray: return (np.abs(arr) < 1e-4) & (arr != 0) +def _fragment_limit(value) -> int: + """``value`` as a positive fragment size in bytes.""" + try: + limit = int(value) + except (TypeError, ValueError): + raise ValueError(f"max_fragment_bytes must be an integer number of bytes, got {value!r}") from None + if limit <= 0: + raise ValueError(f"max_fragment_bytes must be positive, got {limit}") + return limit + + +def _fragment_slices(n: int, item_bytes: int, max_fragment_bytes: int) -> Iterator[list]: + """Per fragment, the ``(start, stop)`` point ranges to serialise and join into it. + + ``item_bytes`` is an upper bound on the text of one point (its floats at their longest plus + separators), so a fragment stays within ``max_fragment_bytes``. A fragment always holds at least one + point: with a limit smaller than one point's text, the fragment is one point long and exceeds it. + A fragment is assembled from :data:`_WORK_BYTES`-sized ranges, which keeps the serialiser's + temporary buffers small whatever the fragment size. + """ + per_fragment = max(1, max_fragment_bytes // max(1, item_bytes)) + per_piece = min(per_fragment, max(1, _WORK_BYTES // max(1, item_bytes))) + for start in range(0, n, per_fragment): + stop = min(start + per_fragment, n) + yield [(a, min(a + per_piece, stop)) for a in range(start, stop, per_piece)] + + def format_floats(values) -> list: """Each value as ``json.dumps`` writes a Python float, with NaN/inf as ``null``. @@ -60,7 +113,7 @@ def format_floats(values) -> list: return [] parts = orjson.dumps(arr, option=orjson.OPT_SERIALIZE_NUMPY)[1:-1].split(b",") for i in np.flatnonzero(_needs_repr(arr)): - parts[i] = repr(float(arr[i])).encode("ascii") + parts[i] = repr(arr[i].item()).encode("ascii") return parts @@ -165,7 +218,7 @@ def coverage(series, i, level): vals = [] for g in series: arr = g["values"].get((p.id, level)) - vals.append(float(arr[i]) if arr is not None else None) + vals.append(arr[i].item() if arr is not None else None) ranges[shortname[p.id]] = _range(shortname[p.id], _py_values(vals)) return { "mars:metadata": first["meta"], @@ -173,8 +226,8 @@ def coverage(series, i, level): "domain": { "type": "Domain", "axes": { - "latitude": {"values": [float(first["lat"][i])]}, - "longitude": {"values": [float(first["lon"][i])]}, + "latitude": {"values": [first["lat"][i].item()]}, + "longitude": {"values": [first["lon"][i].item()]}, "levelist": {"values": levelist}, "t": {"values": t}, }, @@ -209,7 +262,7 @@ def verticalprofile_coverages(groups, header, shortname) -> list: for p in header.parameters: if p.id not in g["values_params"]: continue - vals = [float(g["values"][(p.id, k)][i]) for k in keys] + vals = [g["values"][(p.id, k)][i].item() for k in keys] ranges[shortname[p.id]] = { "type": "NdArray", "dataType": "float", @@ -224,8 +277,8 @@ def verticalprofile_coverages(groups, header, shortname) -> list: "domain": { "type": "Domain", "axes": { - "latitude": {"values": [float(g["lat"][i])]}, - "longitude": {"values": [float(g["lon"][i])]}, + "latitude": {"values": [g["lat"][i].item()]}, + "longitude": {"values": [g["lon"][i].item()]}, "levelist": {"values": levels}, "t": {"values": list(g["t"])}, }, @@ -247,7 +300,7 @@ def trajectory_coverages(groups, header, shortname, coordinates) -> list: keys = list(g["levels"]) if g["levels"] else [None] for k in keys: level = 0 if k is None else k - composite.extend([t, float(la), float(lo), level] for la, lo in zip(g["lat"], g["lon"])) + composite.extend([t, la.item(), lo.item(), level] for la, lo in zip(g["lat"], g["lon"])) for p in header.parameters: arr = g["values"].get((p.id, k)) values[p.id].extend(arr.tolist() if arr is not None else [None] * len(g["lat"])) @@ -291,8 +344,10 @@ class CovjsonStreamEncoder: content_type = "application/prs.coverage+json" file_extension = "covjson" - def __init__(self, config=None): + def __init__(self, config=None, max_fragment_bytes: int = DEFAULT_MAX_FRAGMENT_BYTES): self.config = dict(config or {}) + #: upper bound on one fragment of :meth:`encode_iter` (``config["max_fragment_bytes"]`` wins) + self.max_fragment_bytes = _fragment_limit(self.config.get("max_fragment_bytes") or max_fragment_bytes) self._header = None self.n_coverages = 0 self._streaming = True @@ -327,6 +382,15 @@ def begin(self, header) -> bytes: return b'{"type": "CoverageCollection", "domainType": ' + _dumps(domain) + b', "coverages": [' def encode(self, block) -> bytes: + """The whole block as one ``bytes`` (``b"".join(self.encode_iter(block))``).""" + return b"".join(self.encode_iter(block)) + + def encode_iter(self, block) -> Iterator[bytes]: + """The block as fragments of at most ``max_fragment_bytes``, each from a slice of its arrays. + + The fragments of a block must be consumed completely, in order, before the next block is + encoded: the encoder's state (separators, open ranges) advances as they are produced. + """ if hasattr(block, "lat"): return self._coords(block) if hasattr(block, "values") and hasattr(block, "param"): @@ -348,31 +412,36 @@ def end(self) -> bytes: # -- MultiPoint: stream ------------------------------------------------------------------------------- - def _coords(self, block) -> bytes: + def _coords(self, block) -> Iterator[bytes]: g = block.group if not self._streaming: cur = self._buffer_group(g) cur["lat"].append(np.asarray(block.lat, dtype=np.float64)) cur["lon"].append(np.asarray(block.lon, dtype=np.float64)) - return b"" - out = [] + return if self._current is None or self._current.group is not g: - out.append(self._open_coverage(g)) + yield self._open_coverage(g) cur = self._open() levels = list(g.levels) if g.levels else [0] - out.append(self._tuples(cur, block.lat, block.lon, levels[0])) + yield from self._tuples(cur, block.lat, block.lon, levels[0]) if len(levels) > 1: cur.coords.append((np.array(block.lat, dtype=np.float64), np.array(block.lon, dtype=np.float64))) - return b"".join(out) - @staticmethod - def _tuples(cur: _OpenCoverage, lat, lon, level) -> bytes: - body = composite_tuples(lat, lon, level) - if not body: - return b"" - sep = b", " if cur.n_tuples else b"" - cur.n_tuples += len(lat) - return sep + body + def _tuples(self, cur: _OpenCoverage, lat, lon, level) -> Iterator[bytes]: + """The composite tuples of one band at one level, sliced into fragments. + + The ``(n, 3)`` tuples are built per slice, so the text of a whole band never exists at once. + """ + lat = np.asarray(lat) + lon = np.asarray(lon) + item = 2 * _FLOAT_CHARS + 8 + len(_dumps(level)) # "[lat, lon, level], " + for ranges in _fragment_slices(lat.size, item, self.max_fragment_bytes): + body = b", ".join(composite_tuples(lat[a:b], lon[a:b], level) for a, b in ranges) + if not body: + continue + sep = b", " if cur.n_tuples else b"" + cur.n_tuples += ranges[-1][1] - ranges[0][0] + yield sep + body def _open_coverage(self, g) -> bytes: self._current = _OpenCoverage(g) @@ -384,43 +453,45 @@ def _open_coverage(self, g) -> bytes: self.n_coverages += 1 return head - def _close_composite(self) -> bytes: + def _close_composite(self) -> Iterator[bytes]: cur = self._open() if cur.composite_closed: - return b"" + return cur.composite_closed = True - out = [] levels = list(cur.group.levels) for level in levels[1:]: for lat, lon in cur.coords: - out.append(self._tuples(cur, lat, lon, level)) + yield from self._tuples(cur, lat, lon, level) cur.coords = [] - out.append(b']}}}, "ranges": {') - return b"".join(out) + yield b']}}}, "ranges": {' - def _values(self, block) -> bytes: + def _values(self, block) -> Iterator[bytes]: g = block.group if not self._streaming: cur = self._buffer_group(g) cur["values"].setdefault((block.param, block.level), []).append(np.asarray(block.values, dtype=np.float64)) - return b"" + return cur = self._open() - out = [self._close_composite()] + yield from self._close_composite() if cur.param != block.param: name = self._shortname.get(block.param, block.param) n = g.n_points * max(1, len(g.levels)) - out.append(b"]}, " if cur.param is not None else b"") - out.append(_dumps(name) + b': {"type": "NdArray", "dataType": "float", "shape": [' + str(n).encode()) - out.append(b'], "axisNames": [' + _dumps(str(name)) + b'], "values": [') + head = b"]}, " if cur.param is not None else b"" + head += _dumps(name) + b': {"type": "NdArray", "dataType": "float", "shape": [' + str(n).encode() + head += b'], "axisNames": [' + _dumps(str(name)) + b'], "values": [' cur.param = block.param cur.n_values = 0 - body = float_list_bytes(block.values) - if body: - out.append((b", " if cur.n_values else b"") + body) - cur.n_values += len(block.values) - return b"".join(out) + yield head + values = np.asarray(block.values) + for ranges in _fragment_slices(values.size, _FLOAT_CHARS + 2, self.max_fragment_bytes): + body = b", ".join(float_list_bytes(values[a:b]) for a, b in ranges) + if not body: + continue + sep = b", " if cur.n_values else b"" + cur.n_values += ranges[-1][1] - ranges[0][0] + yield sep + body - def _group_end(self, block) -> bytes: + def _group_end(self, block) -> Iterator[bytes]: if not self._streaming: cur = self._buffered[-1] if self._buffered else None if cur is not None and cur["group"] is block.group: @@ -428,14 +499,13 @@ def _group_end(self, block) -> bytes: cur["lon"] = np.concatenate(cur["lon"]) if cur["lon"] else np.empty(0) cur["values"] = {k: np.concatenate(v) for k, v in cur["values"].items()} cur["done"] = True - return b"" + return cur = self._current if cur is None or cur.group is not block.group: - return b"" - out = [self._close_composite()] - out.append(b"]}}}" if cur.param is not None else b"}}") + return + yield from self._close_composite() + yield b"]}}}" if cur.param is not None else b"}}" self._current = None - return b"".join(out) # -- buffered domains --------------------------------------------------------------------------------- diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py index dcccd73..d25ffcc 100644 --- a/tests/test_stream_encoder.py +++ b/tests/test_stream_encoder.py @@ -14,6 +14,7 @@ import pytest from covjsonkit.stream import ( + DEFAULT_MAX_FRAGMENT_BYTES, CovjsonStreamEncoder, composite_tuples, float_list_bytes, @@ -121,6 +122,16 @@ def encode(header, blocks) -> bytes: return b"".join(out) +def fragments(header, blocks, max_fragment_bytes) -> list: + """Every fragment of the document, block by block, with ``encode_iter``.""" + enc = CovjsonStreamEncoder({"param_db": "ecmwf"}, max_fragment_bytes=max_fragment_bytes) + out = [enc.begin(header)] + for block in blocks: + out += list(enc.encode_iter(block)) + out.append(enc.end()) + return out + + def legacy_multipoint(header, groups): """What BoundingBox.from_polytope + json.dumps produced for the same data.""" coords = legacy_referencing(header) @@ -195,6 +206,58 @@ def test_band_size_invariance(n_bands): assert ref == expected +def _mixed_fields(n, levels, seed=11): + """Values covering both float paths: a NaN, a ``1e-5``-style number, ``-0.0`` and plain floats.""" + rng = np.random.default_rng(seed) + fields = {(p, lev): np.round(rng.normal(size=n) * 1e3, 9) for p in ("165", "167") for lev in levels or (None,)} + first = fields[("165", levels[0] if levels else None)] + first[0], first[1], first[2] = np.nan, 1.5e-5, -0.0 + return fields + + +@pytest.mark.parametrize("max_fragment_bytes", [64, 1024, DEFAULT_MAX_FRAGMENT_BYTES], ids=["64B", "1KiB", "8MiB"]) +def test_fragments_join_to_the_single_fragment_bytes(max_fragment_bytes): + levels = ("500", "850") + lat, lon = _data(500, 8) + lon[4] = 3e-5 # forces the per-value path for the composite tuples as well + fields = _mixed_fields(500, levels) + blocks = multipoint_blocks(fields, lat, lon, levels=levels, n_bands=3) + joined = b"".join(fragments(BBOX, blocks, max_fragment_bytes)) + assert joined == encode(BBOX, blocks) + group = (fields, lat, lon, levels, {"class": "od", "number": 0}, ("2024-01-01T00:00:00Z",)) + assert joined == legacy_multipoint(BBOX, [group]) + + +@pytest.mark.parametrize("max_fragment_bytes", [1024, DEFAULT_MAX_FRAGMENT_BYTES], ids=["1KiB", "8MiB"]) +def test_no_fragment_exceeds_the_limit(max_fragment_bytes): + lat, lon = _data(5000, 9) + fields = {("167", None): np.linspace(200.0, 320.0, 5000)} + frags = fragments(BBOX, multipoint_blocks(fields, lat, lon, n_bands=2), max_fragment_bytes) + assert max(len(f) for f in frags) <= max_fragment_bytes + assert all(f for f in frags) # no empty fragment is handed out + assert len(frags) > 5 # the arrays really were sliced + + +def test_small_magnitude_values_are_formatted_per_value_in_every_slice(): + vals = np.array([1e-5, 2e-5, 0.00002, 9.99e-5, 3.0, np.nan]) + lat, lon = _data(6, 10) + blocks = multipoint_blocks({("167", None): vals}, lat, lon, n_bands=2) + out = b"".join(fragments(BBOX, blocks, 64)) + assert out == encode(BBOX, blocks) + assert b"1e-05" in out and b"2e-05" in out and b"1e-5" not in out + assert json.loads(out)["coverages"][0]["ranges"]["2t"]["values"] == [1e-5, 2e-5, 2e-5, 9.99e-5, 3.0, None] + + +def test_fragment_limit_is_configurable_and_validated(): + assert CovjsonStreamEncoder().max_fragment_bytes == DEFAULT_MAX_FRAGMENT_BYTES + assert CovjsonStreamEncoder({"max_fragment_bytes": 4096}).max_fragment_bytes == 4096 + assert CovjsonStreamEncoder(max_fragment_bytes=4096).max_fragment_bytes == 4096 + with pytest.raises(ValueError): + CovjsonStreamEncoder(max_fragment_bytes=-1) + with pytest.raises(ValueError): + CovjsonStreamEncoder({"max_fragment_bytes": "many"}) + + def test_nan_is_null_and_missing_param_is_omitted(): lat, lon = _data(4, 3) vals = np.array([1.0, np.nan, 3.0, np.nan]) @@ -292,6 +355,17 @@ def test_pointseries_layout_is_point_major_over_consecutive_groups(): assert doc["referencing"][0]["coordinates"] == ["latitude", "longitude", "levelist"] +def test_buffered_domains_produce_their_coverages_from_end(): + """Point features are buffered: their blocks yield no fragment, ``end()`` writes the coverages.""" + header = Header("timeseries", "PointSeries", "date", (T2,)) + meta = {"class": "od", "number": 0, "levelist": 0} + blocks = _point_group(0, ("2024-01-01T00:00:00Z",), meta, {("167", None): [1, 2]}) + enc = CovjsonStreamEncoder(max_fragment_bytes=64) + begin = enc.begin(header) + assert [f for block in blocks for f in enc.encode_iter(block)] == [] + assert begin + enc.end() == encode(header, blocks) + + def test_verticalprofile_layout(): header = Header("verticalprofile", "VerticalProfile", "date", (T2,)) blocks = _point_group( From d4400ea4cbd6567c8e1920374a5664502bbb38f0 Mon Sep 17 00:00:00 2001 From: james Date: Thu, 8 Oct 2026 00:43:36 +0100 Subject: [PATCH 06/12] Bound the peak RSS of encoding a 5M-point block to ~40 MB The measurement lives in tests/stream_memory_probe.py (run it directly for the numbers) and runs in a subprocess so the reported peak is only this encoding; it resets the process' peak RSS once the arrays are allocated, so a parent's footprint cannot mask the growth. 5M points -> 318 MB of CoverageJSON in 54 fragments of at most 8 MiB, for ~41 MB of peak RSS above the blocks themselves; 1M points cost ~30 MB, 20M points ~49 MB. The test asserts <100 MB and that growth does not follow the block size. --- tests/stream_memory_probe.py | 121 +++++++++++++++++++++++++++++++++++ tests/test_stream_memory.py | 40 ++++++++++++ 2 files changed, 161 insertions(+) create mode 100644 tests/stream_memory_probe.py create mode 100644 tests/test_stream_memory.py diff --git a/tests/stream_memory_probe.py b/tests/stream_memory_probe.py new file mode 100644 index 0000000..f8e88f2 --- /dev/null +++ b/tests/stream_memory_probe.py @@ -0,0 +1,121 @@ +"""Peak RSS of streaming one large MultiPoint coverage, as a subprocess probe. + +Run it directly to measure; it prints one JSON line:: + + python tests/stream_memory_probe.py [n_points] + +The blocks are duck-typed stand-ins for polytope-mars' ``CoordsBlock`` / ``ValuesBlock`` / ``GroupEnd``. +The process' peak RSS is reset once the arrays are allocated (``/proc/self/clear_refs``, so the peak a +parent process left behind does not mask the measurement), hence the reported growth is what encoding +costs on top of the data itself. ``test_stream_memory.py`` asserts the bounds. +""" + +from __future__ import annotations + +import argparse +import gc +import json +import resource +from types import SimpleNamespace + +import numpy as np + +from covjsonkit.stream import CovjsonStreamEncoder + + +def peak_rss_bytes() -> int: + """Peak resident set size of this process since the last :func:`reset_peak_rss`.""" + try: + with open("/proc/self/status", "r") as status: + for line in status: + if line.startswith("VmHWM:"): + return int(line.split()[1]) * 1024 + except OSError: + pass + return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss * 1024 + + +def reset_peak_rss() -> bool: + """Forget the peak RSS so far (Linux >= 4.0); ``False`` when the kernel does not allow it.""" + try: + with open("/proc/self/clear_refs", "w") as clear_refs: + clear_refs.write("5") + except OSError: + return False + return True + + +def one_coverage(n_points: int): + """Header and blocks of one coverage of ``n_points`` points with one parameter.""" + param = SimpleNamespace(id="167", shortname="2t", name="2 metre temperature", unit="K", description="") + header = SimpleNamespace( + feature_type="boundingbox", + domain_type="MultiPoint", + time_axis="date", + parameters=(param,), + mars_metadata={}, + extra={}, + ) + group = SimpleNamespace( + index=0, + path={}, + t=("2024-01-01T00:00:00Z",), + params=("167",), + levels=(), + n_points=n_points, + n_bands=1, + mars_metadata={"class": "od", "stream": "oper", "number": 0, "step": 0}, + ) + coords = SimpleNamespace( + group=group, + band=0, + offset=0, + lat=np.linspace(-90.0, 90.0, n_points), + lon=np.linspace(0.0, 359.9, n_points), + ) + values = SimpleNamespace( + group=group, + param="167", + level=None, + band=0, + offset=0, + values=np.linspace(200.0, 320.0, n_points), + ) + return header, [coords, values, SimpleNamespace(group=group)] + + +def measure(n_points: int) -> dict: + header, blocks = one_coverage(n_points) + encoder = CovjsonStreamEncoder() + gc.collect() + peak_reset = reset_peak_rss() + baseline = peak_rss_bytes() + + total = len(encoder.begin(header)) + n_fragments = 0 + largest = 0 + for block in blocks: + for fragment in encoder.encode_iter(block): + n_fragments += 1 + largest = max(largest, len(fragment)) + total += len(fragment) + total += len(encoder.end()) + + peak = peak_rss_bytes() + return { + "n_points": n_points, + "document_bytes": total, + "n_fragments": n_fragments, + "largest_fragment_bytes": largest, + "max_fragment_bytes": encoder.max_fragment_bytes, + "peak_was_reset": peak_reset, + "baseline_rss_bytes": baseline, + "peak_rss_bytes": peak, + "peak_rss_growth_bytes": peak - baseline, + } + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Measure the peak RSS of streaming one large coverage.") + parser.add_argument("n_points", nargs="?", type=int, default=5_000_000) + print(json.dumps(measure(parser.parse_args().n_points))) diff --git a/tests/test_stream_memory.py b/tests/test_stream_memory.py new file mode 100644 index 0000000..a1a6f0e --- /dev/null +++ b/tests/test_stream_memory.py @@ -0,0 +1,40 @@ +"""The streaming encoder's memory bound: peak RSS must not follow the size of a block. + +The measurement runs in a subprocess (``stream_memory_probe.py``) so the peak RSS it reports is only +this encoding; run that module directly to see the numbers. +""" + +import json +import subprocess +import sys +from pathlib import Path + +import pytest + +from covjsonkit.stream import DEFAULT_MAX_FRAGMENT_BYTES + +PROBE = str(Path(__file__).with_name("stream_memory_probe.py")) +MB = 1024 * 1024 + + +def _probe(n_points: int) -> dict: + run = subprocess.run([sys.executable, PROBE, str(n_points)], capture_output=True, text=True) + assert run.returncode == 0, run.stdout + run.stderr + last = run.stdout.strip().splitlines()[-1] if run.stdout.strip() else "" + try: + return json.loads(last) + except json.JSONDecodeError as exc: + raise AssertionError(f"probe printed no measurement: {last!r}\n{run.stderr}") from exc + + +def test_peak_rss_growth_is_bounded_and_independent_of_the_block_size(): + small = _probe(1_000_000) + big = _probe(5_000_000) + if not (small["peak_was_reset"] and big["peak_was_reset"]): + pytest.skip("this kernel does not support resetting the peak RSS via /proc/self/clear_refs") + assert big["document_bytes"] > 300 * MB # ~320 MB of JSON text out of 5M points + assert big["largest_fragment_bytes"] <= DEFAULT_MAX_FRAGMENT_BYTES + assert small["largest_fragment_bytes"] <= DEFAULT_MAX_FRAGMENT_BYTES + assert big["peak_rss_growth_bytes"] < 100 * MB, big + # five times the points for 5 times the text, but not for more memory + assert big["peak_rss_growth_bytes"] < small["peak_rss_growth_bytes"] + 32 * MB, f"{small} vs {big}" From 6f3b8db83c228171c66676da07319a57fee519d4 Mon Sep 17 00:00:00 2001 From: james Date: Thu, 8 Oct 2026 00:44:38 +0100 Subject: [PATCH 07/12] Document the fragment API and its memory bound in CHANGES.md --- CHANGES.md | 36 +++++++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index d79bf92..67d943e 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -6,16 +6,24 @@ Per-branch note for the PR description. The legacy encoder API (`Covjsonkit().en ## New: `covjsonkit.stream.CovjsonStreamEncoder` A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars.blocks`, DESIGN §3): -`begin(header) -> bytes`, `encode(block) -> bytes`, `end() -> bytes`, `content_type = -"application/prs.coverage+json"`, `file_extension = "covjson"`. polytope-mars selects it for `format: covjson` -(lazy import, `polytope_mars.encoders`). +`begin(header) -> bytes`, `encode_iter(block) -> Iterator[bytes]`, `encode(block) -> bytes`, +`end() -> bytes`, `content_type = "application/prs.coverage+json"`, `file_extension = "covjson"`. +polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.encoders`). - Blocks are read by attribute only; the module imports neither polytope-feature, polytope-mars nor covjson-pydantic (tested). Parameter metadata comes resolved in the header (`ParamInfo`), so the encoder does no `param_db` lookups. +- **Fragments:** `encode_iter(block)` yields the block one slice of its arrays at a time, each fragment at most + `max_fragment_bytes` (constructor argument or `encoders.covjson.max_fragment_bytes`, default 8 MiB), so the + text of a whole block never exists at once: a whole-world O2560 coverage is ~800 MB of coordinates and + ~470 MB of values, and encoding a 5M-point block costs ~40 MB of peak RSS whatever its size + (`tests/stream_memory_probe.py`). The fragments of a block must be consumed completely and in order before + the next block is encoded. `encode(block)` returns `b"".join(encode_iter(block))` for buffered callers; + both give the same bytes. Structural fragments (a coverage's `mars:metadata` and `t`, the collection's + `referencing` and `parameters`) are written whole - they are a few hundred bytes. - **MultiPoint** coverages are written as blocks arrive: coverage opening + `mars:metadata` + `t` on the first coordinates block, `[lat, lon, level]` tuples (levels outer, points inner), each range opened on its first - values block and streamed band by band, closed on `GroupEnd`. Memory per coverage: one block, plus the + values block and streamed band by band, closed on `GroupEnd`. Memory per coverage: one fragment, plus the coordinates when the coverage has several levels (they are written once per level). - **PointSeries, VerticalProfile, Trajectory** collections are buffered and written by `end()` (point features are small and the legacy layout is point-major across field groups). The PointSeries layout is the single @@ -24,9 +32,9 @@ A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars - `referencing` and the composite `coordinates` reproduce what each legacy encoder method wrote (`legacy_referencing`: `latitude/longitude/levelist`, `x/y/z` or `t/x/y/z` by feature and time-axis role). - **Bytes:** identical to `json.dumps(legacy_covjson).encode()` (separators `", "`/`": "`, ASCII escapes, - key order `type, domainType, coverages, referencing, parameters`). Fragments with strings (metadata, - parameters, referencing) use `json.dumps`; numeric arrays use `orjson` with `OPT_SERIALIZE_NUMPY` over the - whole array, adjusted to `json.dumps` separators with `bytes.replace`. orjson and CPython agree on the + key order `type, domainType, coverages, referencing, parameters`). Text with strings in it (metadata, + parameters, referencing) uses `json.dumps`; numeric arrays use `orjson` with `OPT_SERIALIZE_NUMPY` per + slice, adjusted to `json.dumps` separators with `bytes.replace`. orjson and CPython agree on the shortest round-trip digits; they differ only in how numbers with 0 < |x| < 1e-4 are spelled (`1e-5` vs `1e-05`, `0.00002` vs `2e-05`), and arrays containing such values are formatted per value with `repr`. `-0.0`, integer-valued floats and `1e+16`-style exponents already match (tested on 5,000 random @@ -44,7 +52,13 @@ A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars ## Tests -`tests/test_stream_encoder.py` (24 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, -band-size invariance (1, 2, 3 and 13 bands, with levels), `null` for NaN and omitted ranges, float and -composite formatting against `json.dumps` (fast and per-value paths), referencing quirks, PointSeries, -VerticalProfile and Trajectory layouts, no polytope imports. +`tests/test_stream_encoder.py` (32 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, +band-size invariance (1, 2, 3 and 13 bands, with levels), fragment joins equal to the single-fragment and +legacy bytes at 64 B / 1 KiB / 8 MiB fragment limits, the fragment limit itself, `null` for NaN and omitted +ranges, float and composite formatting against `json.dumps` (fast and per-value paths), referencing quirks, +PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. + +`tests/test_stream_memory.py` runs `tests/stream_memory_probe.py` in a subprocess: a 5M-point +`CoordsBlock` + `ValuesBlock` produce 318 MB of CoverageJSON in 54 fragments of at most 8 MiB for ~41 MB of +peak RSS above the blocks (1M points: ~30 MB; 20M points: ~49 MB), asserted below 100 MB and against any +growth with the block size. From 109b08117a4d0ddc9c69141125927572763daa1b Mon Sep 17 00:00:00 2001 From: james Date: Fri, 9 Oct 2026 00:12:03 +0100 Subject: [PATCH 08/12] Describe the streaming encoder by what it writes, not by a design section CHANGES.md names the block stream and the missing-param rule the encoder implements instead of citing an external design document, and opens with the compatibility statement a reader of the branch needs. --- CHANGES.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 67d943e..680c002 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,11 +1,11 @@ # Changes on `feat/streaming-encoder` -Per-branch note for the PR description. The legacy encoder API (`Covjsonkit().encode(...)`, -`from_polytope*`, decoders, `param_db`) is unchanged; all existing tests pass. +The legacy encoder API (`Covjsonkit().encode(...)`, `from_polytope*`, decoders, `param_db`) is +unchanged; all existing tests pass. ## New: `covjsonkit.stream.CovjsonStreamEncoder` -A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars.blocks`, DESIGN §3): +A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars.blocks`): `begin(header) -> bytes`, `encode_iter(block) -> Iterator[bytes]`, `encode(block) -> bytes`, `end() -> bytes`, `content_type = "application/prs.coverage+json"`, `file_extension = "covjson"`. polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.encoders`). @@ -44,7 +44,7 @@ polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.enco - `NaN`/`inf` values are written as `null` (legacy `json.dumps` wrote the invalid JSON token `NaN`). - A param absent from a field group (gribjump had no message) has no range in that coverage; a group with no - param at all produces no coverage (decided by polytope-mars, DESIGN §2.5). The collection's `parameters` + param at all produces no coverage, as polytope-mars reports it. The collection's `parameters` always lists every parameter of the header. - Everything else (including the legacy quirks: space-separated datetimes on the `_step` path, the `referencing` variants, int `realization`, efcl without `Forecast date`) is byte-identical; see the From 0a98f5d492283e4f8597101b34034a9e7041111c Mon Sep 17 00:00:00 2001 From: james Date: Fri, 9 Oct 2026 09:34:06 +0100 Subject: [PATCH 09/12] Read blocks by arrival order and bound fragments once --- CHANGES.md | 17 ++++++++++++-- covjsonkit/stream.py | 38 ++++++++++++++----------------- tests/stream_memory_probe.py | 5 ---- tests/test_stream_encoder.py | 44 +++++++++++++++++------------------- 4 files changed, 53 insertions(+), 51 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 680c002..0b174c2 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -52,8 +52,8 @@ polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.enco ## Tests -`tests/test_stream_encoder.py` (32 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, -band-size invariance (1, 2, 3 and 13 bands, with levels), fragment joins equal to the single-fragment and +`tests/test_stream_encoder.py` (29 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, +the same bytes for a group split over several blocks, fragment joins equal to the single-fragment and legacy bytes at 64 B / 1 KiB / 8 MiB fragment limits, the fragment limit itself, `null` for NaN and omitted ranges, float and composite formatting against `json.dumps` (fast and per-value paths), referencing quirks, PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. @@ -62,3 +62,16 @@ PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. `CoordsBlock` + `ValuesBlock` produce 318 MB of CoverageJSON in 54 fragments of at most 8 MiB for ~41 MB of peak RSS above the blocks (1M points: ~30 MB; 20M points: ~49 MB), asserted below 100 MB and against any growth with the block size. + +## One fragment bound, and blocks without band attributes + +- **The block IR has no `band` / `offset` / `n_bands`.** polytope-mars fetches a field whole, so those + attributes were always 0, 0 and 1; the encoder never read them (it places a block by arrival order and + its own value counters). A producer that does split a group's points over several consecutive blocks is + still supported and still writes the same bytes -- `_OpenCoverage.coords` keeps the coordinates of a + multi-level coverage per block -- and `tests/test_stream_encoder.py` pins it. +- **`encode_iter` slices a block once, not twice.** `_fragment_slices` yields one `(start, stop)` range per + fragment and the serialiser is called once per fragment. The 256 KiB sub-slicing inside a fragment did + not lower the peak: joining the pieces materialises the whole fragment while the pieces are still alive, + which is the same 2x a single `orjson.dumps` of the fragment costs. `max_fragment_bytes` is the bound, + and the ~41 MB peak for a 318 MB document (`tests/test_stream_memory.py`) is unchanged. diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py index 1fd58c0..0bf34b3 100644 --- a/covjsonkit/stream.py +++ b/covjsonkit/stream.py @@ -22,10 +22,12 @@ * header: ``feature_type``, ``domain_type``, ``time_axis``, ``parameters`` (each with ``id``, ``shortname``, ``name``, ``unit``, ``description``), ``extra``; -* group (``block.group``): ``index``, ``t``, ``params``, ``levels``, ``n_points``, ``n_bands``, - ``mars_metadata``; -* coordinates block: ``lat``, ``lon`` (float arrays), ``band``; values block: ``param``, ``level``, - ``band``, ``values`` (float array, NaN = missing); group end: only ``group``. +* group (``block.group``): ``index``, ``t``, ``params``, ``levels``, ``n_points``, ``mars_metadata``; +* coordinates block: ``lat``, ``lon`` (float arrays); values block: ``param``, ``level``, ``values`` + (float array, NaN = missing); group end: only ``group``. + +A group's points may arrive in one block or in several consecutive ones: blocks are placed by arrival +order, so a producer that splits a field writes the same bytes as one that does not. MultiPoint coverages are written as the blocks arrive (memory independent of the coverage size, apart from the coordinates of the current coverage when it has several levels). PointSeries, @@ -55,8 +57,6 @@ #: Longest text one float64 can take (``-1.2345678901234567e-308``), used to size array slices. _FLOAT_CHARS = 24 -#: Text built in one serialiser call; a fragment is assembled from pieces of this size. -_WORK_BYTES = 256 * 1024 _CRS = {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"} _LATLON = ["latitude", "longitude", "levelist"] @@ -85,20 +85,16 @@ def _fragment_limit(value) -> int: return limit -def _fragment_slices(n: int, item_bytes: int, max_fragment_bytes: int) -> Iterator[list]: - """Per fragment, the ``(start, stop)`` point ranges to serialise and join into it. +def _fragment_slices(n: int, item_bytes: int, max_fragment_bytes: int) -> Iterator[tuple]: + """Per fragment, the ``(start, stop)`` point range to serialise into it. ``item_bytes`` is an upper bound on the text of one point (its floats at their longest plus separators), so a fragment stays within ``max_fragment_bytes``. A fragment always holds at least one point: with a limit smaller than one point's text, the fragment is one point long and exceeds it. - A fragment is assembled from :data:`_WORK_BYTES`-sized ranges, which keeps the serialiser's - temporary buffers small whatever the fragment size. """ per_fragment = max(1, max_fragment_bytes // max(1, item_bytes)) - per_piece = min(per_fragment, max(1, _WORK_BYTES // max(1, item_bytes))) for start in range(0, n, per_fragment): - stop = min(start + per_fragment, n) - yield [(a, min(a + per_piece, stop)) for a in range(start, stop, per_piece)] + yield start, min(start + per_fragment, n) def format_floats(values) -> list: @@ -428,19 +424,19 @@ def _coords(self, block) -> Iterator[bytes]: cur.coords.append((np.array(block.lat, dtype=np.float64), np.array(block.lon, dtype=np.float64))) def _tuples(self, cur: _OpenCoverage, lat, lon, level) -> Iterator[bytes]: - """The composite tuples of one band at one level, sliced into fragments. + """The composite tuples of one block of points at one level, sliced into fragments. - The ``(n, 3)`` tuples are built per slice, so the text of a whole band never exists at once. + The ``(n, 3)`` tuples are built per slice, so the text of a whole block never exists at once. """ lat = np.asarray(lat) lon = np.asarray(lon) item = 2 * _FLOAT_CHARS + 8 + len(_dumps(level)) # "[lat, lon, level], " - for ranges in _fragment_slices(lat.size, item, self.max_fragment_bytes): - body = b", ".join(composite_tuples(lat[a:b], lon[a:b], level) for a, b in ranges) + for a, b in _fragment_slices(lat.size, item, self.max_fragment_bytes): + body = composite_tuples(lat[a:b], lon[a:b], level) if not body: continue sep = b", " if cur.n_tuples else b"" - cur.n_tuples += ranges[-1][1] - ranges[0][0] + cur.n_tuples += b - a yield sep + body def _open_coverage(self, g) -> bytes: @@ -483,12 +479,12 @@ def _values(self, block) -> Iterator[bytes]: cur.n_values = 0 yield head values = np.asarray(block.values) - for ranges in _fragment_slices(values.size, _FLOAT_CHARS + 2, self.max_fragment_bytes): - body = b", ".join(float_list_bytes(values[a:b]) for a, b in ranges) + for a, b in _fragment_slices(values.size, _FLOAT_CHARS + 2, self.max_fragment_bytes): + body = float_list_bytes(values[a:b]) if not body: continue sep = b", " if cur.n_values else b"" - cur.n_values += ranges[-1][1] - ranges[0][0] + cur.n_values += b - a yield sep + body def _group_end(self, block) -> Iterator[bytes]: diff --git a/tests/stream_memory_probe.py b/tests/stream_memory_probe.py index f8e88f2..d692f21 100644 --- a/tests/stream_memory_probe.py +++ b/tests/stream_memory_probe.py @@ -63,13 +63,10 @@ def one_coverage(n_points: int): params=("167",), levels=(), n_points=n_points, - n_bands=1, mars_metadata={"class": "od", "stream": "oper", "number": 0, "step": 0}, ) coords = SimpleNamespace( group=group, - band=0, - offset=0, lat=np.linspace(-90.0, 90.0, n_points), lon=np.linspace(0.0, 359.9, n_points), ) @@ -77,8 +74,6 @@ def one_coverage(n_points: int): group=group, param="167", level=None, - band=0, - offset=0, values=np.linspace(200.0, 320.0, n_points), ) return header, [coords, values, SimpleNamespace(group=group)] diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py index d25ffcc..00ea95e 100644 --- a/tests/test_stream_encoder.py +++ b/tests/test_stream_encoder.py @@ -50,15 +50,12 @@ class Group: params: tuple levels: tuple n_points: int - n_bands: int mars_metadata: dict @dataclass(frozen=True) class Coords: group: Any - band: int - offset: int lat: np.ndarray lon: np.ndarray @@ -68,8 +65,6 @@ class Values: group: Any param: str level: Any - band: int - offset: int values: np.ndarray @@ -99,17 +94,20 @@ def _nan_to_none(values): return [None if v != v else v for v in np.asarray(values, dtype=np.float64).tolist()] -def multipoint_blocks(fields, lat, lon, levels=(), meta=None, t=("2024-01-01T00:00:00Z",), n_bands=1, index=0): - """Blocks of one MultiPoint group; ``fields`` = {(param, level): values}.""" +def multipoint_blocks(fields, lat, lon, levels=(), meta=None, t=("2024-01-01T00:00:00Z",), n_blocks=1, index=0): + """Blocks of one MultiPoint group; ``fields`` = {(param, level): values}. + + ``n_blocks`` splits the group's points over that many consecutive blocks, which a producer that + cannot fetch a field whole would do and which must not change the bytes. + """ params = tuple(dict.fromkeys(p for p, _ in fields)) - bounds = np.array_split(np.arange(len(lat)), n_bands) - g = Group(index, {}, t, params, tuple(levels), len(lat), n_bands, meta or {"class": "od", "number": 0}) - starts = [idx[0].item() if len(idx) else 0 for idx in bounds] - blocks: list[Any] = [Coords(g, b, starts[b], lat[idx], lon[idx]) for b, idx in enumerate(bounds)] + bounds = np.array_split(np.arange(len(lat)), n_blocks) + g = Group(index, {}, t, params, tuple(levels), len(lat), meta or {"class": "od", "number": 0}) + blocks: list[Any] = [Coords(g, lat[idx], lon[idx]) for idx in bounds] for p in params: for lev in levels or (None,): - for b, idx in enumerate(bounds): - blocks.append(Values(g, p, lev, b, starts[b], fields[(p, lev)][idx])) + for idx in bounds: + blocks.append(Values(g, p, lev, fields[(p, lev)][idx])) blocks.append(End(g)) return blocks @@ -194,13 +192,13 @@ def test_multipoint_matches_legacy_json_dumps(): json.loads(out) -@pytest.mark.parametrize("n_bands", [1, 2, 3, 13]) -def test_band_size_invariance(n_bands): +def test_several_blocks_per_group_give_the_same_bytes(): + """Blocks are placed by arrival order, so a group split over several blocks encodes identically.""" lat, lon = _data(13, 1) levels = ("500", "850") fields = {(p, lev): np.random.default_rng(2).normal(size=13) * 1e3 for p in ("165", "167") for lev in levels} - ref = encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels, n_bands=1)) - assert encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels, n_bands=n_bands)) == ref + ref = encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels)) + assert encode(BBOX, multipoint_blocks(fields, lat, lon, levels=levels, n_blocks=3)) == ref group = (fields, lat, lon, levels, {"class": "od", "number": 0}, ("2024-01-01T00:00:00Z",)) expected = legacy_multipoint(BBOX, [group]) assert ref == expected @@ -221,7 +219,7 @@ def test_fragments_join_to_the_single_fragment_bytes(max_fragment_bytes): lat, lon = _data(500, 8) lon[4] = 3e-5 # forces the per-value path for the composite tuples as well fields = _mixed_fields(500, levels) - blocks = multipoint_blocks(fields, lat, lon, levels=levels, n_bands=3) + blocks = multipoint_blocks(fields, lat, lon, levels=levels, n_blocks=3) joined = b"".join(fragments(BBOX, blocks, max_fragment_bytes)) assert joined == encode(BBOX, blocks) group = (fields, lat, lon, levels, {"class": "od", "number": 0}, ("2024-01-01T00:00:00Z",)) @@ -232,7 +230,7 @@ def test_fragments_join_to_the_single_fragment_bytes(max_fragment_bytes): def test_no_fragment_exceeds_the_limit(max_fragment_bytes): lat, lon = _data(5000, 9) fields = {("167", None): np.linspace(200.0, 320.0, 5000)} - frags = fragments(BBOX, multipoint_blocks(fields, lat, lon, n_bands=2), max_fragment_bytes) + frags = fragments(BBOX, multipoint_blocks(fields, lat, lon, n_blocks=2), max_fragment_bytes) assert max(len(f) for f in frags) <= max_fragment_bytes assert all(f for f in frags) # no empty fragment is handed out assert len(frags) > 5 # the arrays really were sliced @@ -241,7 +239,7 @@ def test_no_fragment_exceeds_the_limit(max_fragment_bytes): def test_small_magnitude_values_are_formatted_per_value_in_every_slice(): vals = np.array([1e-5, 2e-5, 0.00002, 9.99e-5, 3.0, np.nan]) lat, lon = _data(6, 10) - blocks = multipoint_blocks({("167", None): vals}, lat, lon, n_bands=2) + blocks = multipoint_blocks({("167", None): vals}, lat, lon, n_blocks=2) out = b"".join(fragments(BBOX, blocks, 64)) assert out == encode(BBOX, blocks) assert b"1e-05" in out and b"2e-05" in out and b"1e-5" not in out @@ -328,10 +326,10 @@ def test_multipoint_referencing_quirks(feature, role, coords): def _point_group(index, t, meta, values, levels=(), lat=(51.5, -33.9), lon=(0.1, 18.4)): params = tuple(dict.fromkeys(p for p, _ in values)) - g = Group(index, {}, t, params, tuple(levels), len(lat), 1, meta) - blocks: list[Any] = [Coords(g, 0, 0, np.array(lat), np.array(lon))] + g = Group(index, {}, t, params, tuple(levels), len(lat), meta) + blocks: list[Any] = [Coords(g, np.array(lat), np.array(lon))] for (p, lev), v in values.items(): - blocks.append(Values(g, p, lev, 0, 0, np.asarray(v, dtype=float))) + blocks.append(Values(g, p, lev, np.asarray(v, dtype=float))) return blocks + [End(g)] From b6cf8135311c3530442758d3b367f40a75bf883e Mon Sep 17 00:00:00 2001 From: james Date: Fri, 9 Oct 2026 09:47:58 +0100 Subject: [PATCH 10/12] Write the referencing coordinates the header carries --- CHANGES.md | 17 ++++++++++++----- covjsonkit/stream.py | 30 ++---------------------------- tests/stream_memory_probe.py | 1 + tests/test_stream_encoder.py | 29 ++++++++++++----------------- 4 files changed, 27 insertions(+), 50 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 0b174c2..bb2cfac 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -29,8 +29,10 @@ polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.enco are small and the legacy layout is point-major across field groups). The PointSeries layout is the single function `pointseries_coverages` (to port covjsonkit `develop` #137 there); consecutive groups with equal `mars_metadata` and `levels` form one series (one coverage per point and level). -- `referencing` and the composite `coordinates` reproduce what each legacy encoder method wrote - (`legacy_referencing`: `latitude/longitude/levelist`, `x/y/z` or `t/x/y/z` by feature and time-axis role). +- `referencing` and every coverage's composite `coordinates` are the names the header carries + (`header.referencing_coordinates`), which is how the legacy variants + (`latitude/longitude/levelist`, `x/y/z`, `t/x/y/z`) are reproduced without the encoder knowing + anything about feature types. - **Bytes:** identical to `json.dumps(legacy_covjson).encode()` (separators `", "`/`": "`, ASCII escapes, key order `type, domainType, coverages, referencing, parameters`). Text with strings in it (metadata, parameters, referencing) uses `json.dumps`; numeric arrays use `orjson` with `OPT_SERIALIZE_NUMPY` per @@ -52,11 +54,11 @@ polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.enco ## Tests -`tests/test_stream_encoder.py` (29 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, +`tests/test_stream_encoder.py` (26 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, the same bytes for a group split over several blocks, fragment joins equal to the single-fragment and legacy bytes at 64 B / 1 KiB / 8 MiB fragment limits, the fragment limit itself, `null` for NaN and omitted -ranges, float and composite formatting against `json.dumps` (fast and per-value paths), referencing quirks, -PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. +ranges, float and composite formatting against `json.dumps` (fast and per-value paths), the header's +referencing coordinates, PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. `tests/test_stream_memory.py` runs `tests/stream_memory_probe.py` in a subprocess: a 5M-point `CoordsBlock` + `ValuesBlock` produce 318 MB of CoverageJSON in 54 fragments of at most 8 MiB for ~41 MB of @@ -75,3 +77,8 @@ growth with the block size. not lower the peak: joining the pieces materialises the whole fragment while the pieces are still alive, which is the same 2x a single `orjson.dumps` of the fragment costs. `max_fragment_bytes` is the bound, and the ~41 MB peak for a 318 MB document (`tests/test_stream_memory.py`) is unchanged. +- **The referencing table moved to polytope-mars** (`polytope_mars.legacy_format.referencing_coordinates`, + beside the other rules that reproduce the legacy output). The header carries the names as + `referencing_coordinates` and the encoder writes them, so `pointseries_coverages` is the only place left + that reasons about a feature type -- which is what makes the port of #137 onto `develop` a + one-function change. `covjsonkit.stream.legacy_referencing` is gone. diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py index 0bf34b3..4ff3bdd 100644 --- a/covjsonkit/stream.py +++ b/covjsonkit/stream.py @@ -21,7 +21,7 @@ polytope-mars: * header: ``feature_type``, ``domain_type``, ``time_axis``, ``parameters`` (each with ``id``, - ``shortname``, ``name``, ``unit``, ``description``), ``extra``; + ``shortname``, ``name``, ``unit``, ``description``), ``referencing_coordinates``, ``extra``; * group (``block.group``): ``index``, ``t``, ``params``, ``levels``, ``n_points``, ``mars_metadata``; * coordinates block: ``lat``, ``lon`` (float arrays); values block: ``param``, ``level``, ``values`` (float array, NaN = missing); group end: only ``group``. @@ -47,7 +47,6 @@ __all__ = [ "CovjsonStreamEncoder", "DEFAULT_MAX_FRAGMENT_BYTES", - "legacy_referencing", "pointseries_coverages", ] @@ -59,8 +58,6 @@ _CRS = {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"} -_LATLON = ["latitude", "longitude", "levelist"] -_XYZ = ["x", "y", "z"] def _dumps(obj) -> bytes: @@ -140,29 +137,6 @@ def _py_values(values) -> list: return [None if (v is None or (isinstance(v, float) and not math.isfinite(v))) else v for v in values] -def legacy_referencing(header) -> list: - """The ``coordinates`` of the reference system each legacy encoder wrote (preserved quirk).""" - domain, feature, role = header.domain_type, header.feature_type, header.time_axis - if domain == "Trajectory": - return list(_LATLON) if role == "hdate" else ["t", "x", "y", "z"] - if domain == "VerticalProfile": - return list(_LATLON) - if domain == "PointSeries": - if role == "step": - return list(_XYZ) - if role == "month" and feature == "position": - return list(_XYZ) - return list(_LATLON) - # MultiPoint - if role == "hdate": - return list(_LATLON) - if role == "date": - return list(_LATLON) if feature in ("boundingbox", "circle") else list(_XYZ) - if role == "month": - return list(_LATLON) if feature == "circle" else list(_XYZ) - return list(_XYZ) - - def _parameter(p) -> dict: return { "type": "Parameter", @@ -367,7 +341,7 @@ def begin(self, header) -> bytes: self._header = header self.n_coverages = 0 self._shortname = {p.id: p.shortname for p in header.parameters} - self._coordinates = legacy_referencing(header) + self._coordinates = list(header.referencing_coordinates) self._streaming = header.domain_type == "MultiPoint" or ( header.domain_type == "Trajectory" and header.time_axis == "hdate" ) diff --git a/tests/stream_memory_probe.py b/tests/stream_memory_probe.py index d692f21..dabb899 100644 --- a/tests/stream_memory_probe.py +++ b/tests/stream_memory_probe.py @@ -54,6 +54,7 @@ def one_coverage(n_points: int): time_axis="date", parameters=(param,), mars_metadata={}, + referencing_coordinates=("latitude", "longitude", "levelist"), extra={}, ) group = SimpleNamespace( diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py index 00ea95e..b1efbfb 100644 --- a/tests/test_stream_encoder.py +++ b/tests/test_stream_encoder.py @@ -19,7 +19,6 @@ composite_tuples, float_list_bytes, format_floats, - legacy_referencing, ) @@ -38,6 +37,8 @@ class Header: domain_type: str time_axis: str parameters: tuple + #: polytope-mars decides these (``legacy_format.referencing_coordinates``); the encoder writes them + referencing_coordinates: tuple = ("latitude", "longitude", "levelist") mars_metadata: dict = field(default_factory=dict) extra: dict = field(default_factory=dict) @@ -132,7 +133,7 @@ def fragments(header, blocks, max_fragment_bytes) -> list: def legacy_multipoint(header, groups): """What BoundingBox.from_polytope + json.dumps produced for the same data.""" - coords = legacy_referencing(header) + coords = list(header.referencing_coordinates) covs = [] for fields, lat, lon, levels, meta, t in groups: composite = [[float(a), float(o), lev] for lev in (levels or (0,)) for a, o in zip(lat, lon)] @@ -308,20 +309,14 @@ def test_composite_and_value_bytes_match_json_dumps(small, level): assert composite_tuples(np.empty(0), np.empty(0), level) == b"" -@pytest.mark.parametrize( - "feature,role,coords", - [ - ("boundingbox", "date", ["latitude", "longitude", "levelist"]), - ("polygon", "date", ["x", "y", "z"]), - ("polygon", "step", ["x", "y", "z"]), - ("polygon", "hdate", ["latitude", "longitude", "levelist"]), - ("boundingbox", "month", ["x", "y", "z"]), - ("circle", "month", ["latitude", "longitude", "levelist"]), - ], -) -def test_multipoint_referencing_quirks(feature, role, coords): - referencing = legacy_referencing(Header(feature, "MultiPoint", role, ())) - assert referencing == coords +@pytest.mark.parametrize("coords", [("latitude", "longitude", "levelist"), ("x", "y", "z")]) +def test_the_header_decides_the_referencing_coordinates(coords): + """Which names a collection declares is polytope-mars' rule; the encoder writes what it is given.""" + header = Header("boundingbox", "MultiPoint", "date", (T2,), referencing_coordinates=coords) + lat, lon = _data(3, 7) + doc = json.loads(encode(header, multipoint_blocks({("167", None): np.ones(3)}, lat, lon))) + assert doc["referencing"][0]["coordinates"] == list(coords) + assert doc["coverages"][0]["domain"]["axes"]["composite"]["coordinates"] == list(coords) def _point_group(index, t, meta, values, levels=(), lat=(51.5, -33.9), lon=(0.1, 18.4)): @@ -376,7 +371,7 @@ def test_verticalprofile_layout(): def test_trajectory_layout(): - header = Header("trajectory", "Trajectory", "date", (T2,)) + header = Header("trajectory", "Trajectory", "date", (T2,), referencing_coordinates=("t", "x", "y", "z")) blocks = _point_group(0, (0,), {"number": 0}, {("167", None): [1, 2]}) doc = json.loads(encode(header, blocks)) (cov,) = doc["coverages"] From a4c5d5986c76bbb7daeb5ff9a546c660153c7116 Mon Sep 17 00:00:00 2001 From: james Date: Fri, 9 Oct 2026 12:09:18 +0100 Subject: [PATCH 11/12] Remove the branch notes file; the pull request describes the change --- CHANGES.md | 84 ------------------------------------------------------ 1 file changed, 84 deletions(-) delete mode 100644 CHANGES.md diff --git a/CHANGES.md b/CHANGES.md deleted file mode 100644 index bb2cfac..0000000 --- a/CHANGES.md +++ /dev/null @@ -1,84 +0,0 @@ -# Changes on `feat/streaming-encoder` - -The legacy encoder API (`Covjsonkit().encode(...)`, `from_polytope*`, decoders, `param_db`) is -unchanged; all existing tests pass. - -## New: `covjsonkit.stream.CovjsonStreamEncoder` - -A CoverageJSON encoder over the block stream polytope-mars emits (`polytope_mars.blocks`): -`begin(header) -> bytes`, `encode_iter(block) -> Iterator[bytes]`, `encode(block) -> bytes`, -`end() -> bytes`, `content_type = "application/prs.coverage+json"`, `file_extension = "covjson"`. -polytope-mars selects it for `format: covjson` (lazy import, `polytope_mars.encoders`). - -- Blocks are read by attribute only; the module imports neither polytope-feature, polytope-mars nor - covjson-pydantic (tested). Parameter metadata comes resolved in the header (`ParamInfo`), so the encoder does - no `param_db` lookups. -- **Fragments:** `encode_iter(block)` yields the block one slice of its arrays at a time, each fragment at most - `max_fragment_bytes` (constructor argument or `encoders.covjson.max_fragment_bytes`, default 8 MiB), so the - text of a whole block never exists at once: a whole-world O2560 coverage is ~800 MB of coordinates and - ~470 MB of values, and encoding a 5M-point block costs ~40 MB of peak RSS whatever its size - (`tests/stream_memory_probe.py`). The fragments of a block must be consumed completely and in order before - the next block is encoded. `encode(block)` returns `b"".join(encode_iter(block))` for buffered callers; - both give the same bytes. Structural fragments (a coverage's `mars:metadata` and `t`, the collection's - `referencing` and `parameters`) are written whole - they are a few hundred bytes. -- **MultiPoint** coverages are written as blocks arrive: coverage opening + `mars:metadata` + `t` on the first - coordinates block, `[lat, lon, level]` tuples (levels outer, points inner), each range opened on its first - values block and streamed band by band, closed on `GroupEnd`. Memory per coverage: one fragment, plus the - coordinates when the coverage has several levels (they are written once per level). -- **PointSeries, VerticalProfile, Trajectory** collections are buffered and written by `end()` (point features - are small and the legacy layout is point-major across field groups). The PointSeries layout is the single - function `pointseries_coverages` (to port covjsonkit `develop` #137 there); consecutive groups with equal - `mars_metadata` and `levels` form one series (one coverage per point and level). -- `referencing` and every coverage's composite `coordinates` are the names the header carries - (`header.referencing_coordinates`), which is how the legacy variants - (`latitude/longitude/levelist`, `x/y/z`, `t/x/y/z`) are reproduced without the encoder knowing - anything about feature types. -- **Bytes:** identical to `json.dumps(legacy_covjson).encode()` (separators `", "`/`": "`, ASCII escapes, - key order `type, domainType, coverages, referencing, parameters`). Text with strings in it (metadata, - parameters, referencing) uses `json.dumps`; numeric arrays use `orjson` with `OPT_SERIALIZE_NUMPY` per - slice, adjusted to `json.dumps` separators with `bytes.replace`. orjson and CPython agree on the - shortest round-trip digits; they differ only in how numbers with 0 < |x| < 1e-4 are spelled - (`1e-5` vs `1e-05`, `0.00002` vs `2e-05`), and arrays containing such values are formatted per value with - `repr`. `-0.0`, integer-valued floats and `1e+16`-style exponents already match (tested on 5,000 random - magnitudes from 1e-12 to 1e20). - -## Output differences against the legacy encoders - -- `NaN`/`inf` values are written as `null` (legacy `json.dumps` wrote the invalid JSON token `NaN`). -- A param absent from a field group (gribjump had no message) has no range in that coverage; a group with no - param at all produces no coverage, as polytope-mars reports it. The collection's `parameters` - always lists every parameter of the header. -- Everything else (including the legacy quirks: space-separated datetimes on the `_step` path, the - `referencing` variants, int `realization`, efcl without `Forecast date`) is byte-identical; see the - polytope-mars golden corpus and its `CHANGES.md`. - -## Tests - -`tests/test_stream_encoder.py` (26 tests, synthetic blocks, no polytope): legacy-identical MultiPoint bytes, -the same bytes for a group split over several blocks, fragment joins equal to the single-fragment and -legacy bytes at 64 B / 1 KiB / 8 MiB fragment limits, the fragment limit itself, `null` for NaN and omitted -ranges, float and composite formatting against `json.dumps` (fast and per-value paths), the header's -referencing coordinates, PointSeries, VerticalProfile and Trajectory layouts, no polytope imports. - -`tests/test_stream_memory.py` runs `tests/stream_memory_probe.py` in a subprocess: a 5M-point -`CoordsBlock` + `ValuesBlock` produce 318 MB of CoverageJSON in 54 fragments of at most 8 MiB for ~41 MB of -peak RSS above the blocks (1M points: ~30 MB; 20M points: ~49 MB), asserted below 100 MB and against any -growth with the block size. - -## One fragment bound, and blocks without band attributes - -- **The block IR has no `band` / `offset` / `n_bands`.** polytope-mars fetches a field whole, so those - attributes were always 0, 0 and 1; the encoder never read them (it places a block by arrival order and - its own value counters). A producer that does split a group's points over several consecutive blocks is - still supported and still writes the same bytes -- `_OpenCoverage.coords` keeps the coordinates of a - multi-level coverage per block -- and `tests/test_stream_encoder.py` pins it. -- **`encode_iter` slices a block once, not twice.** `_fragment_slices` yields one `(start, stop)` range per - fragment and the serialiser is called once per fragment. The 256 KiB sub-slicing inside a fragment did - not lower the peak: joining the pieces materialises the whole fragment while the pieces are still alive, - which is the same 2x a single `orjson.dumps` of the fragment costs. `max_fragment_bytes` is the bound, - and the ~41 MB peak for a 318 MB document (`tests/test_stream_memory.py`) is unchanged. -- **The referencing table moved to polytope-mars** (`polytope_mars.legacy_format.referencing_coordinates`, - beside the other rules that reproduce the legacy output). The header carries the names as - `referencing_coordinates` and the encoder writes them, so `pointseries_coverages` is the only place left - that reasons about a feature type -- which is what makes the port of #137 onto `develop` a - one-function change. `covjsonkit.stream.legacy_referencing` is gone. From 0cee1e4922df8fd4c8847339a1948e8b7df77d5c Mon Sep 17 00:00:00 2001 From: james Date: Fri, 9 Oct 2026 16:11:16 +0100 Subject: [PATCH 12/12] List in the collection's parameters only the parameters its coverages hold A parameter gribjump has no message for in any field (mx2t6 at hourly steps, a parameter a stream does not produce) has no range in any coverage. The collection's parameters object is written after the coverages, so it can list what the result holds instead of what was requested, and a reader may take it as the list of present parameters. An empty collection has {}. --- covjsonkit/stream.py | 6 +++++- tests/test_stream_encoder.py | 21 ++++++++++++++++++--- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/covjsonkit/stream.py b/covjsonkit/stream.py index 4ff3bdd..ab490a5 100644 --- a/covjsonkit/stream.py +++ b/covjsonkit/stream.py @@ -342,6 +342,9 @@ def begin(self, header) -> bytes: self.n_coverages = 0 self._shortname = {p.id: p.shortname for p in header.parameters} self._coordinates = list(header.referencing_coordinates) + #: ids of the parameters some coverage has a range for; the collection's ``parameters`` lists + #: these, so a reader can take it as the list of what the result holds + self._present: set = set() self._streaming = header.domain_type == "MultiPoint" or ( header.domain_type == "Trajectory" and header.time_axis == "hdate" ) @@ -375,7 +378,7 @@ def end(self) -> bytes: out.append((b", " if self.n_coverages or i else b"") + _dumps(cov)) self.n_coverages += len(covs) self._buffered = [] - parameters = {p.shortname: _parameter(p) for p in self.header.parameters} + parameters = {p.shortname: _parameter(p) for p in self.header.parameters if p.id in self._present} referencing = [{"coordinates": self._coordinates, "system": dict(_CRS)}] out.append(b'], "referencing": ' + _dumps(referencing) + b', "parameters": ' + _dumps(parameters) + b"}") return b"".join(out) @@ -437,6 +440,7 @@ def _close_composite(self) -> Iterator[bytes]: def _values(self, block) -> Iterator[bytes]: g = block.group + self._present.add(block.param) if not self._streaming: cur = self._buffer_group(g) cur["values"].setdefault((block.param, block.level), []).append(np.asarray(block.values, dtype=np.float64)) diff --git a/tests/test_stream_encoder.py b/tests/test_stream_encoder.py index b1efbfb..6f388c8 100644 --- a/tests/test_stream_encoder.py +++ b/tests/test_stream_encoder.py @@ -132,13 +132,18 @@ def fragments(header, blocks, max_fragment_bytes) -> list: def legacy_multipoint(header, groups): - """What BoundingBox.from_polytope + json.dumps produced for the same data.""" + """What BoundingBox.from_polytope + json.dumps produced for the same data. + + The collection's ``parameters`` lists the parameters some coverage has a range for, in header order. + """ coords = list(header.referencing_coordinates) covs = [] + present = set() for fields, lat, lon, levels, meta, t in groups: composite = [[float(a), float(o), lev] for lev in (levels or (0,)) for a, o in zip(lat, lon)] ranges = {} params = tuple(dict.fromkeys(p for p, _ in fields)) + present.update(params) short = {p.id: p.shortname for p in header.parameters} for p in params: vals = [v for lev in (levels or (None,)) for v in _nan_to_none(fields[(p, lev)])] @@ -168,7 +173,7 @@ def legacy_multipoint(header, groups): "domainType": "MultiPoint", "coverages": covs, "referencing": [{"coordinates": coords, "system": CRS}], - "parameters": _parameters(header.parameters), + "parameters": _parameters([p for p in header.parameters if p.id in present]), } return json.dumps(doc).encode() @@ -266,7 +271,7 @@ def test_nan_is_null_and_missing_param_is_omitted(): ranges = doc["coverages"][0]["ranges"] assert list(ranges) == ["2t"] # 10u not in group.params -> no range assert ranges["2t"]["values"] == [1.0, None, 3.0, None] - assert list(doc["parameters"]) == ["10u", "2t"] # header parameters are always listed + assert list(doc["parameters"]) == ["2t"] # a parameter with no range in any coverage is not listed def test_several_groups_and_empty_collection(): @@ -280,6 +285,16 @@ def test_several_groups_and_empty_collection(): ] empty = json.loads(encode(BBOX, [])) assert empty["coverages"] == [] and list(empty) == ["type", "domainType", "coverages", "referencing", "parameters"] + assert empty["parameters"] == {} + + +def test_parameters_lists_what_the_coverages_hold_in_header_order(): + lat, lon = _data(3, 4) + blocks = multipoint_blocks({("167", None): np.ones(3)}, lat, lon, index=0) + blocks += multipoint_blocks({("165", None): np.zeros(3)}, lat, lon, index=1, t=("2024-01-02T00:00:00Z",)) + doc = json.loads(encode(BBOX, blocks)) + assert [list(c["ranges"]) for c in doc["coverages"]] == [["2t"], ["10u"]] + assert list(doc["parameters"]) == ["10u", "2t"] def test_float_formatting_matches_json_dumps():