From 792d3b351980cb76e7dc0bb97e54e0b7f4f2f019 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 2 Jul 2026 10:59:48 +0100 Subject: [PATCH 1/4] Add initial prototype of to_grib functionality, currently only works for bounding box and oper data, need to add grid info from other source as currently hard coded --- covjsonkit/decoder/BoundingBox.py | 194 +++++++++ covjsonkit/decoder/Circle.py | 3 + covjsonkit/decoder/Frame.py | 3 + covjsonkit/decoder/Grid.py | 3 + covjsonkit/decoder/Path.py | 3 + covjsonkit/decoder/Position.py | 3 + covjsonkit/decoder/Shapefile.py | 3 + covjsonkit/decoder/TimeSeries.py | 3 + covjsonkit/decoder/VerticalProfile.py | 3 + covjsonkit/decoder/Wkt.py | 3 + covjsonkit/decoder/decoder.py | 4 + covjsonkit/decoder/grib_backends/__init__.py | 47 +++ covjsonkit/decoder/grib_backends/base.py | 31 ++ .../decoder/grib_backends/eccodes_backend.py | 245 +++++++++++ .../grib_backends/mars2grib_backend.py | 42 ++ covjsonkit/encoder/BoundingBox.py | 4 +- pyproject.toml | 3 + tests/data/test_oper_multipoint_coverage.json | 95 +++++ tests/test_grib_area_comparison.py | 391 ++++++++++++++++++ tests/test_grib_bounding_box.py | 359 ++++++++++++++++ 20 files changed, 1441 insertions(+), 1 deletion(-) create mode 100644 covjsonkit/decoder/grib_backends/__init__.py create mode 100644 covjsonkit/decoder/grib_backends/base.py create mode 100644 covjsonkit/decoder/grib_backends/eccodes_backend.py create mode 100644 covjsonkit/decoder/grib_backends/mars2grib_backend.py create mode 100644 tests/data/test_oper_multipoint_coverage.json create mode 100644 tests/test_grib_area_comparison.py create mode 100644 tests/test_grib_bounding_box.py diff --git a/covjsonkit/decoder/BoundingBox.py b/covjsonkit/decoder/BoundingBox.py index 34e0deb..3ed90e2 100644 --- a/covjsonkit/decoder/BoundingBox.py +++ b/covjsonkit/decoder/BoundingBox.py @@ -1,3 +1,5 @@ +import logging + import numpy as np try: @@ -11,6 +13,16 @@ from ..encoder.encoder import sort_step_values from .decoder import Decoder +logger = logging.getLogger(__name__) + +# Default grid metadata for ECMWF oper data (TCo1279 → O1280 reduced Gaussian). +# Used when mars:grid is absent from the CoverageJSON. Will be replaced once +# polytope exposes real grid information. +_OPER_GRID_DEFAULTS = { + "gridType": "reduced_gg", + "N": 1280, +} + class BoundingBox(Decoder): def __init__(self, covjson): @@ -216,3 +228,185 @@ def to_xarray(self): ds.attrs["date"] = self.get_coordinates()["t"]["values"][0] return ds + + # ------------------------------------------------------------------ + # GRIB export + # ------------------------------------------------------------------ + + def to_grib(self, output_path="output.grib", backend="auto"): + """Convert the CoverageJSON to a multi-message GRIB file. + + Produces one GRIB message per field — i.e. for each unique + combination of (parameter, step, number, date/time) found in + the coverage collection. This mirrors the output of a standard + MARS ``retrieve`` with an ``area`` keyword. + + Args: + output_path: Filesystem path for the output GRIB file. + backend: GRIB encoding backend to use. One of ``"auto"`` + (try pymars2grib first, fall back to eccodes), + ``"mars2grib"``, or ``"eccodes"``. + + Returns: + The *output_path* that was written. + + Raises: + ImportError: If no suitable GRIB backend is available. + """ + from .grib_backends import get_backend + + grib_backend = get_backend(backend) + + messages = [] + + for coverage in self.coverages: + mars_metadata = coverage.get("mars:metadata", {}) + grid_metadata = coverage.get("mars:grid", {}) + + mars_dict = self._build_mars_dict(mars_metadata, coverage) + misc_dict = self._build_misc_dict(grid_metadata, coverage) + + # Compute sort order for N→S, W→E point ordering (MARS convention) + coords = coverage["domain"]["axes"]["composite"]["values"] + sort_idx = self._nswe_sort_indices(coords) + + # One GRIB message per parameter + for param_shortname in self.parameters: + values = coverage["ranges"][param_shortname]["values"] + + # Reorder values to N→S, W→E + sorted_values = [values[i] for i in sort_idx] + + field_mars = {**mars_dict, "param": self._shortname_to_param_id(param_shortname)} + + msg_bytes = grib_backend.encode_message(sorted_values, field_mars, misc_dict) + messages.append(msg_bytes) + + with open(output_path, "wb") as fh: + for msg in messages: + fh.write(msg) + + logger.info("Wrote %d GRIB message(s) to %s", len(messages), output_path) + return output_path + + # ------------------------------------------------------------------ + # Private helpers for to_grib + # ------------------------------------------------------------------ + + @staticmethod + def _build_mars_dict(mars_metadata, coverage): + """Normalise ``mars:metadata`` into the dict expected by GRIB backends. + + Handles the ``Forecast date`` ISO-8601 string that polytope-mars + puts on each coverage, splitting it into separate ``date`` and + ``time`` keys. + """ + mars = {} + + # Direct MARS keys + for key in ("class", "stream", "type", "expver", "levtype", "domain"): + if key in mars_metadata: + mars[key] = mars_metadata[key] + + # Date / time + forecast_date = mars_metadata.get("Forecast date", "") + if forecast_date: + # "2025-06-23T00:00:00Z" → date=20250623, time=0000 + dt_str = str(forecast_date).replace("Z", "") + if "T" in dt_str: + date_part, time_part = dt_str.split("T", 1) + else: + date_part = dt_str + time_part = "0000" + mars["date"] = date_part.replace("-", "") + mars["time"] = time_part.replace(":", "")[:4] + elif "date" in mars_metadata: + mars["date"] = str(mars_metadata["date"]) + if "time" in mars_metadata: + mars["time"] = str(mars_metadata["time"]) + + # Step + if "step" in mars_metadata: + mars["step"] = str(mars_metadata["step"]) + + # Ensemble number + if "number" in mars_metadata: + mars["number"] = str(mars_metadata["number"]) + + # Level — use the z-coordinate from the first composite point + coords = coverage["domain"]["axes"]["composite"]["values"] + if coords and len(coords[0]) > 2: + level = coords[0][2] + if level != 0: + mars["levelist"] = str(int(level)) + + return mars + + @staticmethod + def _build_misc_dict(grid_metadata, coverage): + """Build the ``misc`` dict with grid geometry for the GRIB backend. + + When ``mars:grid`` is not present on the coverage (i.e. polytope + does not yet expose grid info), oper defaults are applied so that + the pipeline can be tested end-to-end. + """ + misc = {} + + if grid_metadata: + misc.update(grid_metadata) + + # Apply oper defaults when grid info is missing + if "gridType" not in misc: + logger.warning( + "No gridType in mars:grid metadata — assuming ECMWF oper defaults " + "(reduced_gg N1280). This will be replaced when polytope provides " + "real grid information." + ) + for key, default in _OPER_GRID_DEFAULTS.items(): + misc.setdefault(key, default) + + # Compute area from coordinates if not already present + if "area" not in misc: + coords = coverage["domain"]["axes"]["composite"]["values"] + if coords: + lats = [c[0] for c in coords] + lons = [c[1] for c in coords] + misc["area"] = [max(lats), min(lons), min(lats), max(lons)] # N/W/S/E + + # Compute Nj (number of latitude rows) and pl (points per row) + # from the coordinates if not already provided. + if "Nj" not in misc or "pl" not in misc: + coords = coverage["domain"]["axes"]["composite"]["values"] + if coords: + from collections import Counter + + lat_counts = Counter(round(c[0], 9) for c in coords) + # Sort latitudes N→S + sorted_lats = sorted(lat_counts.keys(), reverse=True) + misc.setdefault("Nj", len(sorted_lats)) + misc.setdefault("pl", [lat_counts[lat] for lat in sorted_lats]) + + return misc + + @staticmethod + def _nswe_sort_indices(coords): + """Return indices that sort composite coords into N→S, W→E order. + + MARS GRIB convention: first grid point is the north-west corner, + scanning west→east within each latitude row, rows ordered north→south. + """ + # coords is a list of [lat, lon, level] tuples + # Sort by latitude descending (N→S), then longitude ascending (W→E) + indexed = list(enumerate(coords)) + indexed.sort(key=lambda item: (-item[1][0], item[1][1])) + return [i for i, _ in indexed] + + def _shortname_to_param_id(self, shortname): + """Map a parameter shortname (e.g. ``'2t'``) to its numeric param ID.""" + from covjsonkit.param_db import get_param_id_from_db + + try: + return str(get_param_id_from_db(shortname)) + except (KeyError, Exception): + # If the shortname is not in the DB, pass it through as-is + return shortname diff --git a/covjsonkit/decoder/Circle.py b/covjsonkit/decoder/Circle.py index a3a1ae9..cf61f91 100644 --- a/covjsonkit/decoder/Circle.py +++ b/covjsonkit/decoder/Circle.py @@ -158,3 +158,6 @@ def to_xarray(self): ds.attrs["date"] = self.get_coordinates()["t"]["values"][0] return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Frame.py b/covjsonkit/decoder/Frame.py index 67a15b6..a60c35c 100644 --- a/covjsonkit/decoder/Frame.py +++ b/covjsonkit/decoder/Frame.py @@ -121,3 +121,6 @@ def to_xarray(self): ds.attrs[key] = val return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index 1abf68f..b9cd5b0 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -228,3 +228,6 @@ def to_xarray(self): ds[pname].attrs[k] = v return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Path.py b/covjsonkit/decoder/Path.py index 5001a24..4def7f8 100644 --- a/covjsonkit/decoder/Path.py +++ b/covjsonkit/decoder/Path.py @@ -154,3 +154,6 @@ def to_xarray(self): ds.attrs[mars_metadata] = self.mars_metadata[0][mars_metadata] return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Position.py b/covjsonkit/decoder/Position.py index 4331582..3b6d9c0 100644 --- a/covjsonkit/decoder/Position.py +++ b/covjsonkit/decoder/Position.py @@ -213,3 +213,6 @@ def _covers_domain(self, coverage, num, date, x, y, z): and coverage["domain"]["axes"][self.y_name]["values"] == y and coverage["domain"]["axes"][self.z_name]["values"] == z ) + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Shapefile.py b/covjsonkit/decoder/Shapefile.py index 3342e04..730d85b 100644 --- a/covjsonkit/decoder/Shapefile.py +++ b/covjsonkit/decoder/Shapefile.py @@ -124,3 +124,6 @@ def to_xarray(self): ds.attrs[key] = val return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/TimeSeries.py b/covjsonkit/decoder/TimeSeries.py index 28f5702..f3f5526 100644 --- a/covjsonkit/decoder/TimeSeries.py +++ b/covjsonkit/decoder/TimeSeries.py @@ -226,6 +226,9 @@ def _covers_domain(self, coverage, num, date, x, y, z): and coverage["domain"]["axes"][self.z_name]["values"] == z ) + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") + def _to_xarray_no_forecast_date(self): """Convert monthly-means CovJSON (no 'Forecast date' in metadata) to xarray. diff --git a/covjsonkit/decoder/VerticalProfile.py b/covjsonkit/decoder/VerticalProfile.py index 23dc1f4..34a9bf0 100644 --- a/covjsonkit/decoder/VerticalProfile.py +++ b/covjsonkit/decoder/VerticalProfile.py @@ -222,3 +222,6 @@ def to_xarray(self): return ds[0] return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/Wkt.py b/covjsonkit/decoder/Wkt.py index 8f16707..3c7b554 100644 --- a/covjsonkit/decoder/Wkt.py +++ b/covjsonkit/decoder/Wkt.py @@ -121,3 +121,6 @@ def to_xarray(self): ds.attrs[key] = val return ds + + def to_grib(self, output_path="output.grib", backend="auto"): + raise NotImplementedError("to_grib() is only supported for BoundingBox domains.") diff --git a/covjsonkit/decoder/decoder.py b/covjsonkit/decoder/decoder.py index 6922e46..0eb9e8b 100644 --- a/covjsonkit/decoder/decoder.py +++ b/covjsonkit/decoder/decoder.py @@ -122,3 +122,7 @@ def to_geojson(self): @abstractmethod def to_geotiff(self): pass + + @abstractmethod + def to_grib(self, output_path="output.grib", backend="auto"): + pass diff --git a/covjsonkit/decoder/grib_backends/__init__.py b/covjsonkit/decoder/grib_backends/__init__.py new file mode 100644 index 0000000..90c37e1 --- /dev/null +++ b/covjsonkit/decoder/grib_backends/__init__.py @@ -0,0 +1,47 @@ +"""Swappable GRIB encoding backends. + +Use :func:`get_backend` to obtain the best available backend at runtime. +The factory tries ``pymars2grib`` first (preferred) and falls back to +``eccodes`` when it is not installed. +""" + +from .base import GribBackend # noqa: F401 + + +def get_backend(preferred: str = "auto") -> GribBackend: + """Return a ready-to-use GRIB encoding backend. + + Args: + preferred: One of ``"auto"``, ``"mars2grib"``, or ``"eccodes"``. + * ``"auto"`` – try mars2grib first, fall back to eccodes. + * ``"mars2grib"`` – require mars2grib (raises ImportError if + unavailable). + * ``"eccodes"`` – use eccodes directly. + + Returns: + A :class:`GribBackend` instance. + + Raises: + ImportError: If the explicitly requested backend is not installed. + """ + if preferred in ("mars2grib", "auto"): + try: + from .mars2grib_backend import Mars2GribBackend + + return Mars2GribBackend() + except ImportError: + if preferred == "mars2grib": + raise ImportError( + "pymars2grib is not installed. " + "Build metkit from source with pybind11 support, or use backend='eccodes'." + ) + + try: + from .eccodes_backend import EccodesBackend + + return EccodesBackend() + except ImportError: + raise ImportError( + "No GRIB backend available. Install eccodes: pip install eccodes, " + "or build pymars2grib from metkit source." + ) diff --git a/covjsonkit/decoder/grib_backends/base.py b/covjsonkit/decoder/grib_backends/base.py new file mode 100644 index 0000000..8b728c5 --- /dev/null +++ b/covjsonkit/decoder/grib_backends/base.py @@ -0,0 +1,31 @@ +from abc import ABC, abstractmethod + + +class GribBackend(ABC): + """Abstract backend interface for encoding GRIB messages from MARS metadata + values. + + Implementations translate (values, mars_dict, misc_dict) into a single + encoded GRIB message returned as ``bytes``. Two concrete backends are + provided: + + * :class:`~.eccodes_backend.EccodesBackend` – uses the eccodes Python API + (available via ``pip install eccodes``). + * :class:`~.mars2grib_backend.Mars2GribBackend` – wraps ECMWF's + ``pymars2grib`` pybind11 module (requires metkit built from source). + """ + + @abstractmethod + def encode_message(self, values: list, mars: dict, misc: dict) -> bytes: + """Encode a single GRIB message. + + Args: + values: Field data values (one per grid point). + mars: MARS keys describing the field (class, stream, type, date, + time, step, param, levtype, levelist, …). + misc: Non-MARS metadata such as grid geometry (gridType, N, area, + …) and packing options. + + Returns: + The fully encoded GRIB message as raw bytes. + """ + pass diff --git a/covjsonkit/decoder/grib_backends/eccodes_backend.py b/covjsonkit/decoder/grib_backends/eccodes_backend.py new file mode 100644 index 0000000..3a37a95 --- /dev/null +++ b/covjsonkit/decoder/grib_backends/eccodes_backend.py @@ -0,0 +1,245 @@ +"""GRIB encoding backend using the eccodes Python API. + +This is the fallback backend when pymars2grib is not available. +It builds a GRIB2 message from scratch using eccodes sample-based +creation, setting keys derived from MARS metadata and grid geometry. +""" + +import logging + +import eccodes + +from .base import GribBackend + +logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# MARS levtype → eccodes typeOfLevel mapping +# --------------------------------------------------------------------------- +_LEVTYPE_MAP = { + "sfc": "surface", + "pl": "isobaricInhPa", + "ml": "hybrid", + "pv": "potentialVorticity", + "pt": "theta", + "dp": "depthBelowSea", +} + +# --------------------------------------------------------------------------- +# MARS stream → GRIB2 code-table values +# --------------------------------------------------------------------------- +_STREAM_MAP = { + "oper": "oper", + "enfo": "enfo", + "efov": "efov", + "scda": "scda", + "scwv": "scwv", + "wave": "wave", + "waef": "waef", + "moda": "moda", +} + +# --------------------------------------------------------------------------- +# MARS type → GRIB2 code-table values +# --------------------------------------------------------------------------- +_TYPE_MAP = { + "an": "an", + "fc": "fc", + "pf": "pf", + "cf": "cf", + "em": "em", + "es": "es", +} + + +class EccodesBackend(GribBackend): + """Encode GRIB2 messages using the eccodes Python API.""" + + def __init__(self): + # Verify eccodes is functional + try: + eccodes.codes_get_api_version() + except Exception as exc: + raise ImportError(f"eccodes is installed but not functional: {exc}") from exc + + def encode_message(self, values: list, mars: dict, misc: dict) -> bytes: + """Encode a single GRIB2 message. + + Args: + values: Field data values. + mars: MARS keys (class, stream, type, date, time, step, param, + levtype, levelist, number, …). + misc: Grid geometry and encoding options (gridType, N, area, …). + + Returns: + Encoded GRIB2 message as bytes. + """ + sample_id = eccodes.codes_grib_new_from_samples("GRIB2") + try: + # Enable the ECMWF local definition section so that MARS + # namespace keys (marsClass, marsStream, marsType, …) are + # available for reading and writing. + eccodes.codes_set(sample_id, "centre", "ecmf") + eccodes.codes_set_long(sample_id, "setLocalDefinition", 1) + + self._set_identification(sample_id, mars) + self._set_temporal(sample_id, mars) + self._set_product(sample_id, mars) + self._set_grid(sample_id, mars, misc, len(values)) + self._set_ensemble(sample_id, mars) + self._set_values(sample_id, values) + + return eccodes.codes_get_message(sample_id) + finally: + eccodes.codes_release(sample_id) + + # ------------------------------------------------------------------ + # Private helpers + # ------------------------------------------------------------------ + + @staticmethod + def _set_identification(gid: int, mars: dict) -> None: + """Set GRIB identification keys from MARS dict.""" + if "class" in mars: + eccodes.codes_set(gid, "marsClass", mars["class"]) + if "stream" in mars: + stream = _STREAM_MAP.get(mars["stream"], mars["stream"]) + eccodes.codes_set(gid, "marsStream", stream) + if "type" in mars: + mars_type = _TYPE_MAP.get(mars["type"], mars["type"]) + eccodes.codes_set(gid, "marsType", mars_type) + if "expver" in mars: + eccodes.codes_set(gid, "experimentVersionNumber", mars["expver"]) + + @staticmethod + def _set_temporal(gid: int, mars: dict) -> None: + """Set date, time, and step keys.""" + if "date" in mars: + date_str = str(mars["date"]).replace("-", "") + eccodes.codes_set_long(gid, "dataDate", int(date_str)) + if "time" in mars: + time_str = str(mars["time"]).replace(":", "") + # Pad to 4 digits (e.g. "0" → "0000", "12" → "1200") + time_str = time_str.ljust(4, "0") + eccodes.codes_set_long(gid, "dataTime", int(time_str)) + if "step" in mars: + eccodes.codes_set(gid, "stepRange", str(mars["step"])) + + @staticmethod + def _set_product(gid: int, mars: dict) -> None: + """Set parameter and level keys. + + Important: ``typeOfLevel`` and ``level`` must be set BEFORE + ``paramId`` because eccodes re-resolves paramId when the type + of fixed surface changes (e.g. 10v ↔ v depending on surface + vs. upper-air). + """ + levtype = mars.get("levtype", "sfc") + type_of_level = _LEVTYPE_MAP.get(levtype, levtype) + eccodes.codes_set(gid, "typeOfLevel", type_of_level) + + if "levelist" in mars: + eccodes.codes_set_long(gid, "level", int(mars["levelist"])) + elif levtype == "sfc": + eccodes.codes_set_long(gid, "level", 0) + + if "param" in mars: + param = mars["param"] + try: + eccodes.codes_set_long(gid, "paramId", int(param)) + except (ValueError, TypeError): + # If param is a shortname string, set it that way + eccodes.codes_set(gid, "shortName", str(param)) + + @staticmethod + def _set_ensemble(gid: int, mars: dict) -> None: + """Set ensemble-specific keys if present. + + ``productDefinitionTemplateNumber`` must be set to 1 first so + that the ``perturbationNumber`` key becomes available in GRIB2. + """ + if "number" in mars: + number = int(mars["number"]) + if number > 0: + # Switch to ensemble product definition template FIRST + eccodes.codes_set_long(gid, "productDefinitionTemplateNumber", 1) + eccodes.codes_set_long(gid, "perturbationNumber", number) + + def _set_grid(self, gid: int, mars: dict, misc: dict, num_values: int) -> None: + """Set grid definition from misc dict. + + Supports reduced_gg (reduced Gaussian) and regular_ll (regular lat/lon). + Falls back to a simple points-based approach if grid type is unknown. + """ + grid_type = misc.get("gridType", "reduced_gg") + + if grid_type == "reduced_gg": + self._set_reduced_gaussian_grid(gid, misc, num_values) + elif grid_type == "regular_ll": + self._set_regular_ll_grid(gid, misc, num_values) + else: + # Fallback: set as unstructured grid with explicit coordinates + logger.warning("Unknown gridType '%s', falling back to unstructured grid", grid_type) + self._set_unstructured_grid(gid, misc, num_values) + + @staticmethod + def _set_reduced_gaussian_grid(gid: int, misc: dict, num_values: int) -> None: + """Configure a reduced Gaussian grid (sub-area).""" + eccodes.codes_set(gid, "gridType", "reduced_gg") + + n = misc.get("N", 1280) + eccodes.codes_set_long(gid, "N", n) + + # Area bounds (N/W/S/E) in degrees + if "area" in misc: + area = misc["area"] + eccodes.codes_set_double(gid, "latitudeOfFirstGridPointInDegrees", float(area[0])) + eccodes.codes_set_double(gid, "longitudeOfFirstGridPointInDegrees", float(area[1])) + eccodes.codes_set_double(gid, "latitudeOfLastGridPointInDegrees", float(area[2])) + eccodes.codes_set_double(gid, "longitudeOfLastGridPointInDegrees", float(area[3])) + + # Nj (number of latitude rows in the sub-area) + if "Nj" in misc: + eccodes.codes_set_long(gid, "Nj", int(misc["Nj"])) + + # pl array (number of points per latitude row) — required for sub-area + if "pl" in misc: + eccodes.codes_set_array(gid, "pl", misc["pl"]) + + eccodes.codes_set_long(gid, "numberOfDataPoints", num_values) + + @staticmethod + def _set_regular_ll_grid(gid: int, misc: dict, num_values: int) -> None: + """Configure a regular lat/lon grid (sub-area).""" + eccodes.codes_set(gid, "gridType", "regular_ll") + + if "area" in misc: + area = misc["area"] + eccodes.codes_set_double(gid, "latitudeOfFirstGridPointInDegrees", float(area[0])) + eccodes.codes_set_double(gid, "longitudeOfFirstGridPointInDegrees", float(area[1])) + eccodes.codes_set_double(gid, "latitudeOfLastGridPointInDegrees", float(area[2])) + eccodes.codes_set_double(gid, "longitudeOfLastGridPointInDegrees", float(area[3])) + + if "Dx" in misc: + eccodes.codes_set_double(gid, "iDirectionIncrementInDegrees", float(misc["Dx"])) + if "Dy" in misc: + eccodes.codes_set_double(gid, "jDirectionIncrementInDegrees", float(misc["Dy"])) + if "Ni" in misc: + eccodes.codes_set_long(gid, "Ni", int(misc["Ni"])) + if "Nj" in misc: + eccodes.codes_set_long(gid, "Nj", int(misc["Nj"])) + + eccodes.codes_set_long(gid, "numberOfDataPoints", num_values) + + @staticmethod + def _set_unstructured_grid(gid: int, misc: dict, num_values: int) -> None: + """Fallback: configure as an unstructured grid.""" + eccodes.codes_set(gid, "gridType", "unstructured_grid") + eccodes.codes_set_long(gid, "numberOfDataPoints", num_values) + + @staticmethod + def _set_values(gid: int, values: list) -> None: + """Pack field values into the GRIB data section.""" + # Use 16 bits per value (MARS default for oper data) + eccodes.codes_set_long(gid, "bitsPerValue", 16) + eccodes.codes_set_values(gid, values) diff --git a/covjsonkit/decoder/grib_backends/mars2grib_backend.py b/covjsonkit/decoder/grib_backends/mars2grib_backend.py new file mode 100644 index 0000000..727bdfe --- /dev/null +++ b/covjsonkit/decoder/grib_backends/mars2grib_backend.py @@ -0,0 +1,42 @@ +"""GRIB encoding backend using pymars2grib (metkit's Python binding). + +This is the preferred backend when available. pymars2grib handles all +MARS-to-GRIB key resolution internally, so this wrapper is intentionally +thin. + +pymars2grib is not yet available on PyPI — it must be built from metkit +source with pybind11 support. When it is not installed the factory in +``__init__.py`` falls back to the eccodes backend. +""" + +from .base import GribBackend + + +class Mars2GribBackend(GribBackend): + """Encode GRIB2 messages via pymars2grib.""" + + def __init__(self): + try: + from pymars2grib import Mars2Grib + except ImportError: + raise ImportError( + "pymars2grib is not installed. " + "Build metkit from source with pybind11 support to use this backend." + ) + self._encoder = Mars2Grib() + + def encode_message(self, values: list, mars: dict, misc: dict) -> bytes: + """Encode a single GRIB2 message. + + Delegates entirely to pymars2grib which resolves GRIB header + layout from the MARS dictionary and injects the values. + + Args: + values: Field data values. + mars: MARS keys describing the field. + misc: Auxiliary metadata (grid geometry, packing hints, …). + + Returns: + Encoded GRIB2 message as bytes. + """ + return self._encoder.encode(values, mars, misc) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index d80d96f..7285733 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -12,13 +12,15 @@ def __init__(self, type, domaintype): self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - def add_coverage(self, mars_metadata, coords, values): + def add_coverage(self, mars_metadata, coords, values, grid_metadata=None): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) + if grid_metadata: + new_coverage["mars:grid"] = grid_metadata self.add_domain(new_coverage, coords) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) diff --git a/pyproject.toml b/pyproject.toml index a1e954a..4334334 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,9 @@ geo = [ "rasterio", "shapely", ] +grib = [ + "eccodes", +] tests = [ "pytest", "polytope-python", diff --git a/tests/data/test_oper_multipoint_coverage.json b/tests/data/test_oper_multipoint_coverage.json new file mode 100644 index 0000000..3d7fa25 --- /dev/null +++ b/tests/data/test_oper_multipoint_coverage.json @@ -0,0 +1,95 @@ +{ + "type": "CoverageCollection", + "domainType": "MultiPoint", + "coverages": [ + { + "mars:metadata": { + "class": "od", + "Forecast date": "2025-06-23T00:00:00Z", + "domain": "g", + "expver": "0001", + "levtype": "sfc", + "step": 0, + "stream": "oper", + "type": "an" + }, + "type": "Coverage", + "domain": { + "type": "Domain", + "axes": { + "t": { + "values": ["2025-06-23T00:00:00Z"] + }, + "composite": { + "dataType": "tuple", + "coordinates": ["latitude", "longitude", "levelist"], + "values": [ + [0.035149384216, 0.0, 0], + [0.035149384216, 0.070093457944, 0], + [0.035149384216, 0.140186915888, 0], + [0.105448152647, 0.0, 0], + [0.105448152647, 0.070148090413, 0], + [0.105448152647, 0.140296180826, 0], + [0.175746921078, 0.0, 0], + [0.175746921078, 0.070202808112, 0], + [0.175746921078, 0.140405616225, 0] + ] + } + } + }, + "ranges": { + "10v": { + "type": "NdArray", + "dataType": "float", + "shape": [9], + "axisNames": ["10v"], + "values": [ + 6.2300872802734375, 6.3160247802734375, 6.2593841552734375, + 6.1285247802734375, 6.1968841552734375, 6.1265716552734375, + 6.1255950927734375, 6.1304779052734375, 5.9986419677734375 + ] + }, + "2t": { + "type": "NdArray", + "dataType": "float", + "shape": [9], + "axisNames": ["2t"], + "values": [ + 299.01934814453125, 299.06427001953125, 299.14044189453125, + 298.97247314453125, 299.01739501953125, 299.07403564453125, + 298.93731689453125, 298.99786376953125, 299.09356689453125 + ] + } + } + } + ], + "referencing": [ + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84" + } + } + ], + "parameters": { + "10v": { + "type": "Parameter", + "description": {"en": "10 metre V wind component"}, + "unit": {"symbol": "m s**-1"}, + "observedProperty": { + "id": "10v", + "label": {"en": "10 metre V wind component"} + } + }, + "2t": { + "type": "Parameter", + "description": {"en": "2 metre temperature"}, + "unit": {"symbol": "K"}, + "observedProperty": { + "id": "2t", + "label": {"en": "2 metre temperature"} + } + } + } +} diff --git a/tests/test_grib_area_comparison.py b/tests/test_grib_area_comparison.py new file mode 100644 index 0000000..a2997d9 --- /dev/null +++ b/tests/test_grib_area_comparison.py @@ -0,0 +1,391 @@ +"""Integration test: Compare MARS area output vs BoundingBox to_grib() output. + +This test retrieves the same data via two polytope paths: + 1. Polytope mars-od datasource with `area` keyword → GRIB directly + 2. Polytope boundingbox feature request → CoverageJSON → to_grib() + +Then compares the two GRIB files field-by-field to verify our +to_grib() produces output matching what MARS returns natively. + +Requirements: + - Valid polytope credentials (~/.polytopeapirc) + - Network access to polytope.ecmwf.int + - eccodes Python package + +Run: + pytest tests/test_grib_area_comparison.py -v -m integration + # or directly: + python tests/test_grib_area_comparison.py +""" + +import json +import os +import sys +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required") + +# Skip by default unless explicitly requested via marker or env var +pytestmark = pytest.mark.skipif( + os.environ.get("RUN_INTEGRATION_TESTS", "0") != "1", + reason="Integration test — set RUN_INTEGRATION_TESTS=1 to run", +) + +# --------------------------------------------------------------------------- +# Configuration — adjust these for your test case +# --------------------------------------------------------------------------- + +# A small bounding box near the equator (matches our test fixtures) +AREA = [0.18, 0.0, 0.03, 0.15] # N/W/S/E in degrees + +# MARS request keys (common to both retrieval paths) +MARS_REQUEST = { + "class": "od", + "stream": "oper", + "type": "fc", + "levtype": "sfc", + "param": "2t/10v", + "step": "0", + "date": "-1", # yesterday (most likely available) + "time": "0000", +} + +# Polytope collection for standard MARS retrieves +COLLECTION = "ecmwf-mars" + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _read_all_grib_messages(filepath): + """Read all messages from a GRIB file, returning list of dicts with keys+values.""" + messages = [] + with open(filepath, "rb") as f: + while True: + gid = eccodes.codes_grib_new_from_file(f) + if gid is None: + break + try: + msg = { + "paramId": eccodes.codes_get(gid, "paramId"), + "shortName": eccodes.codes_get(gid, "shortName"), + "dataDate": eccodes.codes_get(gid, "dataDate"), + "dataTime": eccodes.codes_get(gid, "dataTime"), + "stepRange": eccodes.codes_get(gid, "stepRange"), + "gridType": eccodes.codes_get(gid, "gridType"), + "numberOfDataPoints": eccodes.codes_get(gid, "numberOfDataPoints"), + "values": eccodes.codes_get_values(gid).tolist(), + } + try: + msg["N"] = eccodes.codes_get(gid, "N") + except Exception: + pass + try: + msg["latitudeOfFirstGridPointInDegrees"] = eccodes.codes_get( + gid, "latitudeOfFirstGridPointInDegrees" + ) + msg["longitudeOfFirstGridPointInDegrees"] = eccodes.codes_get( + gid, "longitudeOfFirstGridPointInDegrees" + ) + msg["latitudeOfLastGridPointInDegrees"] = eccodes.codes_get( + gid, "latitudeOfLastGridPointInDegrees" + ) + msg["longitudeOfLastGridPointInDegrees"] = eccodes.codes_get( + gid, "longitudeOfLastGridPointInDegrees" + ) + except Exception: + pass + messages.append(msg) + finally: + eccodes.codes_release(gid) + return messages + + +def _retrieve_mars_area(output_path): + """Retrieve data via polytope mars-od datasource with area keyword → GRIB.""" + from polytope.api import Client + + client = Client() + + # Dict-based request with area key routes to mars-od datasource + request = {**MARS_REQUEST, "area": f"{AREA[0]}/{AREA[1]}/{AREA[2]}/{AREA[3]}"} + + print(f"\n{'='*60}") + print("PATH 1: Polytope mars-od retrieve with area keyword") + print(f"{'='*60}") + print(f"Collection: {COLLECTION}") + print(f"Request: {request}") + print(f"Output: {output_path}") + + client.retrieve( + COLLECTION, + request, + output_file=output_path, + pointer=False, + ) + + file_size = os.path.getsize(output_path) + print(f"Retrieved {file_size} bytes") + return output_path + + +def _retrieve_boundingbox_covjson(): + """Retrieve data via polytope BoundingBox feature → CoverageJSON dict.""" + from polytope.api import Client + + client = Client() + + # Build dict-based request (polytope requires dict format for feature routing) + request = {**MARS_REQUEST} + request["feature"] = { + "type": "boundingbox", + "points": [[AREA[0], AREA[1]], [AREA[2], AREA[3]]], + } + + print(f"\n{'='*60}") + print("PATH 2: Polytope BoundingBox feature → CoverageJSON") + print(f"{'='*60}") + print(f"Collection: {COLLECTION}") + print(f"Request: {request}") + + # Feature requests return CoverageJSON (write to temp file, then parse) + with tempfile.NamedTemporaryFile(suffix=".covjson", delete=False, mode="w") as tmp: + covjson_path = tmp.name + + try: + client.retrieve( + COLLECTION, + request, + output_file=covjson_path, + ) + + file_size = os.path.getsize(covjson_path) + print(f"Retrieved {file_size} bytes of CoverageJSON") + + with open(covjson_path) as f: + covjson = json.load(f) + + print(f"Type: {covjson.get('type')}") + print(f"Domain type: {covjson.get('domainType')}") + print(f"Coverages: {len(covjson.get('coverages', []))}") + + return covjson + finally: + if os.path.exists(covjson_path): + os.unlink(covjson_path) + + +def _convert_covjson_to_grib(covjson, output_path): + """Convert CoverageJSON to GRIB using our to_grib() method.""" + from covjsonkit.api import Covjsonkit + + print(f"\n{'='*60}") + print("Converting CoverageJSON → GRIB via to_grib()") + print(f"{'='*60}") + + decoder = Covjsonkit().decode(covjson) + result = decoder.to_grib(output_path, backend="eccodes") + + file_size = os.path.getsize(output_path) + print(f"Wrote {file_size} bytes to {output_path}") + return result + + +def _compare_grib_files(mars_grib_path, covjson_grib_path): + """Compare two GRIB files and report differences. + + Returns (match: bool, report: str) + """ + mars_msgs = _read_all_grib_messages(mars_grib_path) + covjson_msgs = _read_all_grib_messages(covjson_grib_path) + + report_lines = [] + report_lines.append(f"\n{'='*60}") + report_lines.append("COMPARISON RESULTS") + report_lines.append(f"{'='*60}") + report_lines.append(f"MARS area messages: {len(mars_msgs)}") + report_lines.append(f"CovJSON→GRIB messages: {len(covjson_msgs)}") + report_lines.append("") + + all_match = True + + # Compare metadata + mars_params = sorted(set(m["shortName"] for m in mars_msgs)) + covjson_params = sorted(set(m["shortName"] for m in covjson_msgs)) + report_lines.append(f"MARS params: {mars_params}") + report_lines.append(f"CovJSON params: {covjson_params}") + + if mars_params != covjson_params: + report_lines.append("⚠️ Parameter mismatch!") + all_match = False + + # Match messages by (shortName, stepRange) and compare + for mars_msg in mars_msgs: + key = (mars_msg["shortName"], mars_msg["stepRange"]) + matching = [m for m in covjson_msgs if (m["shortName"], m["stepRange"]) == key] + + if not matching: + report_lines.append(f"\n❌ No CovJSON match for {key}") + all_match = False + continue + + covjson_msg = matching[0] + report_lines.append(f"\n--- Field: {key[0]}, step={key[1]} ---") + + # Compare number of data points + if mars_msg["numberOfDataPoints"] != covjson_msg["numberOfDataPoints"]: + report_lines.append( + f" ⚠️ Point count differs: MARS={mars_msg['numberOfDataPoints']} " + f"vs CovJSON={covjson_msg['numberOfDataPoints']}" + ) + all_match = False + else: + report_lines.append(f" ✓ Point count: {mars_msg['numberOfDataPoints']}") + + # Compare grid type + report_lines.append(f" Grid: MARS={mars_msg['gridType']} vs CovJSON={covjson_msg['gridType']}") + if mars_msg["gridType"] != covjson_msg["gridType"]: + report_lines.append(" ⚠️ Grid type differs (expected — see notes)") + + # Compare values + mars_vals = mars_msg["values"] + covjson_vals = covjson_msg["values"] + + if len(mars_vals) != len(covjson_vals): + report_lines.append(f" ⚠️ Value count differs: {len(mars_vals)} vs {len(covjson_vals)}") + all_match = False + else: + # Statistics on differences + diffs = [abs(a - b) for a, b in zip(mars_vals, covjson_vals)] + max_diff = max(diffs) if diffs else 0 + mean_diff = sum(diffs) / len(diffs) if diffs else 0 + n_exact = sum(1 for d in diffs if d == 0.0) + + report_lines.append(f" Values ({len(mars_vals)} points):") + report_lines.append(f" Max abs diff: {max_diff:.10e}") + report_lines.append(f" Mean abs diff: {mean_diff:.10e}") + report_lines.append(f" Exact matches: {n_exact}/{len(mars_vals)}") + + # Check if values are "close enough" (within GRIB packing tolerance) + tolerance = 0.001 # 1e-3 — generous for GRIB2 simple packing + if max_diff > tolerance: + report_lines.append(f" ⚠️ Max diff {max_diff} exceeds tolerance {tolerance}") + all_match = False + else: + report_lines.append(f" ✓ All values within tolerance ({tolerance})") + + # Compare area bounds if available + for bound_key in [ + "latitudeOfFirstGridPointInDegrees", + "longitudeOfFirstGridPointInDegrees", + "latitudeOfLastGridPointInDegrees", + "longitudeOfLastGridPointInDegrees", + ]: + if bound_key in mars_msg and bound_key in covjson_msg: + mars_val = mars_msg[bound_key] + covjson_val = covjson_msg[bound_key] + diff = abs(mars_val - covjson_val) + status = "✓" if diff < 0.001 else "⚠️" + report_lines.append(f" {status} {bound_key}: MARS={mars_val:.6f} CovJSON={covjson_val:.6f}") + + # Summary + report_lines.append(f"\n{'='*60}") + if all_match: + report_lines.append("✅ OVERALL: GRIB files match within tolerance") + else: + report_lines.append("❌ OVERALL: Differences found (see details above)") + report_lines.append(f"{'='*60}") + + report = "\n".join(report_lines) + return all_match, report + + +# --------------------------------------------------------------------------- +# Pytest test +# --------------------------------------------------------------------------- + + +class TestGribAreaComparison: + """Compare MARS area retrieve vs BoundingBox→to_grib(). + + Run with: RUN_INTEGRATION_TESTS=1 pytest tests/test_grib_area_comparison.py -v -s + """ + + def test_area_vs_boundingbox_to_grib(self): + """Full roundtrip comparison: MARS area vs CoverageJSON→GRIB.""" + mars_grib = tempfile.NamedTemporaryFile(suffix="_mars_area.grib", delete=False) + covjson_grib = tempfile.NamedTemporaryFile(suffix="_covjson.grib", delete=False) + mars_grib.close() + covjson_grib.close() + + try: + # Path 1: MARS area → GRIB + _retrieve_mars_area(mars_grib.name) + + # Path 2: BoundingBox feature → CoverageJSON → GRIB + covjson = _retrieve_boundingbox_covjson() + _convert_covjson_to_grib(covjson, covjson_grib.name) + + # Compare + match, report = _compare_grib_files(mars_grib.name, covjson_grib.name) + print(report) + + # Don't assert match yet — first run is diagnostic + # Uncomment below once you've validated the output: + # assert match, report + + finally: + for path in [mars_grib.name, covjson_grib.name]: + if os.path.exists(path): + os.unlink(path) + + +# --------------------------------------------------------------------------- +# Standalone execution +# --------------------------------------------------------------------------- + +if __name__ == "__main__": + """Run directly without pytest for quick manual testing.""" + print("GRIB Area Comparison Test") + print(f"Area: N={AREA[0]}, W={AREA[1]}, S={AREA[2]}, E={AREA[3]}") + print(f"Request: {MARS_REQUEST}") + + mars_grib_path = "/tmp/mars_area_output.grib" + covjson_grib_path = "/tmp/covjson_to_grib_output.grib" + + try: + # Path 1: Standard MARS area retrieve + _retrieve_mars_area(mars_grib_path) + + # Path 2: BoundingBox feature → CoverageJSON → to_grib + covjson = _retrieve_boundingbox_covjson() + + # Save CoverageJSON for inspection + covjson_path = "/tmp/boundingbox_covjson.json" + with open(covjson_path, "w") as f: + json.dump(covjson, f, indent=2) + print(f"\nCoverageJSON saved to {covjson_path}") + + _convert_covjson_to_grib(covjson, covjson_grib_path) + + # Compare + match, report = _compare_grib_files(mars_grib_path, covjson_grib_path) + print(report) + + print(f"\nFiles preserved for inspection:") + print(f" MARS area GRIB: {mars_grib_path}") + print(f" CovJSON→GRIB: {covjson_grib_path}") + print(f" CoverageJSON: {covjson_path}") + + sys.exit(0 if match else 1) + + except Exception as e: + print(f"\n❌ Error: {e}") + import traceback + + traceback.print_exc() + sys.exit(2) diff --git a/tests/test_grib_bounding_box.py b/tests/test_grib_bounding_box.py new file mode 100644 index 0000000..bda8303 --- /dev/null +++ b/tests/test_grib_bounding_box.py @@ -0,0 +1,359 @@ +"""Tests for BoundingBox.to_grib() — CoverageJSON → GRIB conversion.""" + +import copy +import json +import os +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required for GRIB tests") + +from covjsonkit.api import Covjsonkit # noqa: E402 + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def oper_covjson(): + """Load the oper (no-ensemble) multipoint test fixture.""" + path = os.path.join(os.path.dirname(__file__), "data", "test_oper_multipoint_coverage.json") + with open(path) as f: + return json.load(f) + + +@pytest.fixture +def enfo_covjson(): + """Load the enfo (ensemble) multipoint test fixture.""" + path = os.path.join(os.path.dirname(__file__), "data", "test_multipoint_coverage.json") + with open(path) as f: + return json.load(f) + + +def _read_grib_messages(filepath): + """Read all GRIB messages from a file and return a list of (keys, values) dicts.""" + messages = [] + with open(filepath, "rb") as f: + while True: + gid = eccodes.codes_grib_new_from_file(f) + if gid is None: + break + try: + msg = { + "paramId": eccodes.codes_get(gid, "paramId"), + "shortName": eccodes.codes_get(gid, "shortName"), + "dataDate": eccodes.codes_get(gid, "dataDate"), + "dataTime": eccodes.codes_get(gid, "dataTime"), + "gridType": eccodes.codes_get(gid, "gridType"), + "numberOfDataPoints": eccodes.codes_get(gid, "numberOfDataPoints"), + "values": eccodes.codes_get_values(gid).tolist(), + } + # Optional keys that may not be in every message + try: + msg["marsClass"] = eccodes.codes_get(gid, "marsClass") + except Exception: + pass + try: + msg["marsStream"] = eccodes.codes_get(gid, "marsStream") + except Exception: + pass + try: + msg["marsType"] = eccodes.codes_get(gid, "marsType") + except Exception: + pass + try: + msg["typeOfLevel"] = eccodes.codes_get(gid, "typeOfLevel") + except Exception: + pass + try: + msg["perturbationNumber"] = eccodes.codes_get(gid, "perturbationNumber") + except Exception: + msg["perturbationNumber"] = 0 + try: + msg["stepRange"] = eccodes.codes_get(gid, "stepRange") + except Exception: + pass + messages.append(msg) + finally: + eccodes.codes_release(gid) + return messages + + +# --------------------------------------------------------------------------- +# Step 6 stub: other decoders raise NotImplementedError +# --------------------------------------------------------------------------- + +class TestToGribStubs: + """Verify that non-BoundingBox decoders raise NotImplementedError.""" + + def test_timeseries_raises(self): + from covjsonkit.decoder.TimeSeries import TimeSeries + with pytest.raises(NotImplementedError): + # We just need to call to_grib on the class — but we can't + # instantiate without valid covjson, so test the method directly. + TimeSeries.to_grib(None) + + def test_vertical_profile_raises(self): + from covjsonkit.decoder.VerticalProfile import VerticalProfile + with pytest.raises(NotImplementedError): + VerticalProfile.to_grib(None) + + +# --------------------------------------------------------------------------- +# Core tests +# --------------------------------------------------------------------------- + +class TestToGribOper: + """Tests using the oper (no-ensemble) fixture.""" + + def test_message_count(self, oper_covjson): + """1 coverage × 2 params = 2 GRIB messages.""" + decoder = Covjsonkit().decode(oper_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + assert len(messages) == 2 + finally: + os.unlink(output) + + def test_roundtrip_values(self, oper_covjson): + """Values in GRIB should match the CoverageJSON input (within packing tolerance).""" + decoder = Covjsonkit().decode(oper_covjson) + original_2t = oper_covjson["coverages"][0]["ranges"]["2t"]["values"] + original_10v = oper_covjson["coverages"][0]["ranges"]["10v"]["values"] + + # Compute expected N→S/W→E sort order (same as to_grib applies) + coords = oper_covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] + from covjsonkit.decoder.BoundingBox import BoundingBox + + sort_idx = BoundingBox._nswe_sort_indices(coords) + expected_2t = [original_2t[i] for i in sort_idx] + expected_10v = [original_10v[i] for i in sort_idx] + + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + + # Find messages by paramId + msg_2t = next(m for m in messages if m["paramId"] == 167) + msg_10v = next(m for m in messages if m["paramId"] == 166) + + assert len(msg_2t["values"]) == len(expected_2t) + assert len(msg_10v["values"]) == len(expected_10v) + + # GRIB packing introduces small rounding — allow tolerance + for orig, decoded in zip(expected_2t, msg_2t["values"]): + assert abs(orig - decoded) < 0.01, f"2t mismatch: {orig} vs {decoded}" + for orig, decoded in zip(expected_10v, msg_10v["values"]): + assert abs(orig - decoded) < 0.01, f"10v mismatch: {orig} vs {decoded}" + finally: + os.unlink(output) + + def test_mars_keys(self, oper_covjson): + """GRIB keys should match the mars:metadata from CoverageJSON.""" + decoder = Covjsonkit().decode(oper_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + + for msg in messages: + assert msg["dataDate"] == 20250623 + assert msg["dataTime"] == 0 + # MARS levtype=sfc maps to different eccodes typeOfLevel + # depending on the parameter (e.g. "surface" for sp, + # "heightAboveGround" for 2t/10v) + assert msg["typeOfLevel"] in ("surface", "heightAboveGround") + assert msg["numberOfDataPoints"] == 9 + finally: + os.unlink(output) + + def test_grid_type_defaults(self, oper_covjson): + """Without mars:grid, the decoder should apply reduced_gg defaults.""" + decoder = Covjsonkit().decode(oper_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + for msg in messages: + assert msg["gridType"] == "reduced_gg" + finally: + os.unlink(output) + + def test_output_path_returned(self, oper_covjson): + """to_grib() should return the output path.""" + decoder = Covjsonkit().decode(oper_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + result = decoder.to_grib(output, backend="eccodes") + assert result == output + finally: + os.unlink(output) + + +class TestToGribEnfo: + """Tests using the enfo (ensemble) fixture.""" + + def test_message_count(self, enfo_covjson): + """4 coverages (2 numbers × 2 steps) × 2 params = 8 GRIB messages.""" + decoder = Covjsonkit().decode(enfo_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + assert len(messages) == 8 + finally: + os.unlink(output) + + def test_ensemble_numbers_present(self, enfo_covjson): + """Each ensemble member should have perturbationNumber set.""" + decoder = Covjsonkit().decode(enfo_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + perturb_nums = {m["perturbationNumber"] for m in messages} + assert 1 in perturb_nums + assert 2 in perturb_nums + finally: + os.unlink(output) + + def test_per_field_uniqueness(self, enfo_covjson): + """Each GRIB message should represent a unique (param, step, number) combo.""" + decoder = Covjsonkit().decode(enfo_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + combos = set() + for msg in messages: + key = (msg["paramId"], msg.get("stepRange", "0"), msg["perturbationNumber"]) + combos.add(key) + # Should have as many unique combos as messages + assert len(combos) == len(messages) + finally: + os.unlink(output) + + +class TestToGribExplicitGrid: + """Test with explicit mars:grid metadata injected into fixture.""" + + def test_explicit_grid_metadata(self, oper_covjson): + """When mars:grid is present, its values should be used instead of defaults.""" + covjson = copy.deepcopy(oper_covjson) + # Inject explicit grid metadata + covjson["coverages"][0]["mars:grid"] = { + "gridType": "regular_ll", + "Ni": 3, + "Nj": 3, + "Dx": 0.07, + "Dy": 0.07, + "area": [0.18, 0.0, 0.03, 0.15], + } + + decoder = Covjsonkit().decode(covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + for msg in messages: + assert msg["gridType"] == "regular_ll" + finally: + os.unlink(output) + + +class TestBackendSelection: + """Test backend factory behaviour.""" + + def test_auto_selects_eccodes(self): + """With pymars2grib unavailable, auto should fall back to eccodes.""" + from covjsonkit.decoder.grib_backends import get_backend + from covjsonkit.decoder.grib_backends.eccodes_backend import EccodesBackend + + backend = get_backend("auto") + assert isinstance(backend, EccodesBackend) + + def test_explicit_eccodes(self): + """Explicitly requesting eccodes should work.""" + from covjsonkit.decoder.grib_backends import get_backend + from covjsonkit.decoder.grib_backends.eccodes_backend import EccodesBackend + + backend = get_backend("eccodes") + assert isinstance(backend, EccodesBackend) + + def test_mars2grib_raises_when_unavailable(self): + """Explicitly requesting mars2grib should raise if not installed.""" + from covjsonkit.decoder.grib_backends import get_backend + + with pytest.raises(ImportError, match="pymars2grib"): + get_backend("mars2grib") + + +class TestBuildHelpers: + """Unit tests for the private helper methods on BoundingBox.""" + + def test_build_mars_dict_date_parsing(self, oper_covjson): + """_build_mars_dict should parse ISO-8601 Forecast date correctly.""" + from covjsonkit.decoder.BoundingBox import BoundingBox + + coverage = oper_covjson["coverages"][0] + mars_metadata = coverage["mars:metadata"] + mars = BoundingBox._build_mars_dict(mars_metadata, coverage) + + assert mars["date"] == "20250623" + assert mars["time"] == "0000" + assert mars["class"] == "od" + assert mars["stream"] == "oper" + assert mars["type"] == "an" + assert mars["step"] == "0" + + def test_build_misc_dict_defaults(self, oper_covjson): + """_build_misc_dict should apply oper defaults when mars:grid is absent.""" + from covjsonkit.decoder.BoundingBox import BoundingBox + + coverage = oper_covjson["coverages"][0] + misc = BoundingBox._build_misc_dict({}, coverage) + + assert misc["gridType"] == "reduced_gg" + assert misc["N"] == 1280 + assert "area" in misc + # Area should be [max_lat, min_lon, min_lat, max_lon] + assert misc["area"][0] > misc["area"][2] # N > S + + def test_build_misc_dict_explicit_grid(self, oper_covjson): + """_build_misc_dict should use explicit grid metadata when provided.""" + from covjsonkit.decoder.BoundingBox import BoundingBox + + coverage = oper_covjson["coverages"][0] + grid_meta = {"gridType": "regular_ll", "Ni": 10, "Nj": 10} + misc = BoundingBox._build_misc_dict(grid_meta, coverage) + + assert misc["gridType"] == "regular_ll" + assert misc["Ni"] == 10 + # Should NOT have oper defaults + assert "N" not in misc + + def test_shortname_to_param_id(self, oper_covjson): + """_shortname_to_param_id should resolve known shortnames.""" + decoder = Covjsonkit().decode(oper_covjson) + assert decoder._shortname_to_param_id("2t") == "167" + assert decoder._shortname_to_param_id("10v") == "166" + + def test_shortname_to_param_id_unknown(self, oper_covjson): + """Unknown shortnames should pass through as-is.""" + decoder = Covjsonkit().decode(oper_covjson) + result = decoder._shortname_to_param_id("totally_unknown_param_xyz") + assert result == "totally_unknown_param_xyz" From e31c2339519c588b55de3846d7881f9cb169685a Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 2 Jul 2026 14:21:58 +0100 Subject: [PATCH 2/4] Now parse grid info that is passed from polytope tree and add it as grid metadata in covjson --- covjsonkit/encoder/BoundingBox.py | 14 ++++++-- covjsonkit/encoder/Circle.py | 2 +- covjsonkit/encoder/Frame.py | 2 +- covjsonkit/encoder/Grid.py | 2 +- covjsonkit/encoder/Path.py | 2 +- covjsonkit/encoder/Position.py | 2 +- covjsonkit/encoder/Shapefile.py | 2 +- covjsonkit/encoder/TimeSeries.py | 2 +- covjsonkit/encoder/VerticalProfile.py | 2 +- covjsonkit/encoder/Wkt.py | 2 +- covjsonkit/encoder/encoder.py | 2 +- tests/test_grib_bounding_box.py | 49 ++++++++++++++++++++++++--- 12 files changed, 66 insertions(+), 17 deletions(-) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index 7285733..e6c4d5b 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -119,8 +119,16 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: - """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (BoundingBox) CoverageJSON collection.""" + def from_polytope(self, result, date_key: str = "date", grid_metadata=None) -> dict: + """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (BoundingBox) CoverageJSON collection. + + Args: + result: The polytope TensorIndexTree result. + date_key: The name of the date axis in the tree. + grid_metadata: Optional dict with grid geometry from polytope config + (e.g. ``{"gridType": "reduced_gg", "N": 1280}``). When provided, + stored as ``mars:grid`` on each coverage. + """ coords = {} mars_metadata = {} range_dict = {} @@ -200,7 +208,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm["number"] = num mm["step"] = normalize_step_value(step) mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step]) + self.add_coverage(mm, coords[date], val_dict[step], grid_metadata=grid_metadata) return self.covjson diff --git a/covjsonkit/encoder/Circle.py b/covjsonkit/encoder/Circle.py index 7c2beb9..4e9ea3d 100644 --- a/covjsonkit/encoder/Circle.py +++ b/covjsonkit/encoder/Circle.py @@ -114,7 +114,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (Circle) CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Frame.py b/covjsonkit/encoder/Frame.py index 22bf4c6..e6f88eb 100644 --- a/covjsonkit/encoder/Frame.py +++ b/covjsonkit/encoder/Frame.py @@ -114,7 +114,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (Frame) CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index 5f3d56d..d9214f2 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -121,7 +121,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a Grid CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Path.py b/covjsonkit/encoder/Path.py index 1e9ff73..2d7d1ba 100644 --- a/covjsonkit/encoder/Path.py +++ b/covjsonkit/encoder/Path.py @@ -114,7 +114,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a Trajectory (Path) CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index 26fef03..f2bb92b 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -123,7 +123,7 @@ def from_xarray(self, datasets): return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a PointSeries (Position) CoverageJSON collection. Args: diff --git a/covjsonkit/encoder/Shapefile.py b/covjsonkit/encoder/Shapefile.py index 897b7fb..77010a2 100644 --- a/covjsonkit/encoder/Shapefile.py +++ b/covjsonkit/encoder/Shapefile.py @@ -114,7 +114,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (Shapefile) CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/TimeSeries.py b/covjsonkit/encoder/TimeSeries.py index 310f630..d77f34c 100644 --- a/covjsonkit/encoder/TimeSeries.py +++ b/covjsonkit/encoder/TimeSeries.py @@ -123,7 +123,7 @@ def from_xarray(self, datasets): return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a PointSeries CoverageJSON collection. Args: diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 0781f89..c942ba1 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -121,7 +121,7 @@ def from_xarray(self, datasets): return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a VerticalProfile CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Wkt.py b/covjsonkit/encoder/Wkt.py index 931e6e3..6c65d23 100644 --- a/covjsonkit/encoder/Wkt.py +++ b/covjsonkit/encoder/Wkt.py @@ -118,7 +118,7 @@ def from_xarray(self, dataset): # Return the generated CoverageJSON return self.covjson - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: """Encode a polytope ``TensorIndexTree`` result into a MultiPoint (Wkt/Polygon) CoverageJSON collection.""" coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/encoder.py b/covjsonkit/encoder/encoder.py index 0332de2..48a3940 100644 --- a/covjsonkit/encoder/encoder.py +++ b/covjsonkit/encoder/encoder.py @@ -768,7 +768,7 @@ def from_xarray(self, dataset): pass @abstractmethod - def from_polytope(self, result, date_key: str = "date") -> dict: + def from_polytope(self, result, date_key: str = "date", **kwargs) -> dict: pass def from_polytope_reforecast(self, result) -> dict: diff --git a/tests/test_grib_bounding_box.py b/tests/test_grib_bounding_box.py index bda8303..4c9df67 100644 --- a/tests/test_grib_bounding_box.py +++ b/tests/test_grib_bounding_box.py @@ -248,12 +248,36 @@ def test_per_field_uniqueness(self, enfo_covjson): class TestToGribExplicitGrid: - """Test with explicit mars:grid metadata injected into fixture.""" + """Test with explicit mars:grid metadata (as polytope-mars would provide from config).""" - def test_explicit_grid_metadata(self, oper_covjson): - """When mars:grid is present, its values should be used instead of defaults.""" + def test_octahedral_grid_from_config(self, oper_covjson): + """mars:grid from octahedral mapper config (type=octahedral, resolution=1280).""" covjson = copy.deepcopy(oper_covjson) - # Inject explicit grid metadata + # This is what polytope-mars _get_grid_metadata() produces for: + # axis_name: values, transformations: [{name: mapper, type: octahedral, resolution: 1280}] + covjson["coverages"][0]["mars:grid"] = { + "gridType": "reduced_gg", + "N": 1280, + } + + decoder = Covjsonkit().decode(covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + for msg in messages: + assert msg["gridType"] == "reduced_gg" + assert msg["numberOfDataPoints"] == 9 + finally: + os.unlink(output) + + def test_regular_ll_grid_from_config(self, oper_covjson): + """mars:grid from local_regular mapper config (e.g. EFAS).""" + covjson = copy.deepcopy(oper_covjson) + # This is what polytope-mars _get_grid_metadata() produces for: + # axis_name: values, transformations: [{name: mapper, type: local_regular, + # resolution: [2969, 4529], local: [22.76, 72.24, -25.24, 50.24]}] covjson["coverages"][0]["mars:grid"] = { "gridType": "regular_ll", "Ni": 3, @@ -274,6 +298,23 @@ def test_explicit_grid_metadata(self, oper_covjson): finally: os.unlink(output) + def test_no_grid_metadata_uses_defaults(self, oper_covjson): + """Without mars:grid, decoder falls back to reduced_gg N=1280 defaults.""" + # oper_covjson has no mars:grid key — verify it uses defaults + assert "mars:grid" not in oper_covjson["coverages"][0] + + decoder = Covjsonkit().decode(oper_covjson) + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + try: + decoder.to_grib(output, backend="eccodes") + messages = _read_grib_messages(output) + for msg in messages: + # Defaults applied + assert msg["gridType"] == "reduced_gg" + finally: + os.unlink(output) + class TestBackendSelection: """Test backend factory behaviour.""" From 44454b0a9e5f94d2b36aed9433b275000478e645 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 2 Jul 2026 14:31:19 +0100 Subject: [PATCH 3/4] Formatting --- covjsonkit/decoder/grib_backends/mars2grib_backend.py | 3 +-- covjsonkit/encoder/Position.py | 1 - covjsonkit/encoder/TimeSeries.py | 1 - covjsonkit/encoder/VerticalProfile.py | 1 - tests/test_grib_area_comparison.py | 6 ++---- tests/test_grib_bounding_box.py | 6 +++++- 6 files changed, 8 insertions(+), 10 deletions(-) diff --git a/covjsonkit/decoder/grib_backends/mars2grib_backend.py b/covjsonkit/decoder/grib_backends/mars2grib_backend.py index 727bdfe..6438292 100644 --- a/covjsonkit/decoder/grib_backends/mars2grib_backend.py +++ b/covjsonkit/decoder/grib_backends/mars2grib_backend.py @@ -20,8 +20,7 @@ def __init__(self): from pymars2grib import Mars2Grib except ImportError: raise ImportError( - "pymars2grib is not installed. " - "Build metkit from source with pybind11 support to use this backend." + "pymars2grib is not installed. " "Build metkit from source with pybind11 support to use this backend." ) self._encoder = Mars2Grib() diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index f2bb92b..7b1e8a1 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -98,7 +98,6 @@ def from_xarray(self, datasets): self.add_parameter(data_var) for dataset in datasets: - # Process each "number" in the dataset for num in dataset["number"].values: dv_dict = {} diff --git a/covjsonkit/encoder/TimeSeries.py b/covjsonkit/encoder/TimeSeries.py index d77f34c..8fb13d1 100644 --- a/covjsonkit/encoder/TimeSeries.py +++ b/covjsonkit/encoder/TimeSeries.py @@ -98,7 +98,6 @@ def from_xarray(self, datasets): self.add_parameter(data_var) for dataset in datasets: - # Process each "number" in the dataset for num in dataset["number"].values: dv_dict = {} diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index c942ba1..3bf58d3 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -94,7 +94,6 @@ def from_xarray(self, datasets): self.add_parameter(data_var) for dataset in datasets: - # Process each "number" in the dataset for num in dataset["number"].values: for step in dataset["time"].values: diff --git a/tests/test_grib_area_comparison.py b/tests/test_grib_area_comparison.py index a2997d9..475f31d 100644 --- a/tests/test_grib_area_comparison.py +++ b/tests/test_grib_area_comparison.py @@ -91,9 +91,7 @@ def _read_all_grib_messages(filepath): msg["longitudeOfFirstGridPointInDegrees"] = eccodes.codes_get( gid, "longitudeOfFirstGridPointInDegrees" ) - msg["latitudeOfLastGridPointInDegrees"] = eccodes.codes_get( - gid, "latitudeOfLastGridPointInDegrees" - ) + msg["latitudeOfLastGridPointInDegrees"] = eccodes.codes_get(gid, "latitudeOfLastGridPointInDegrees") msg["longitudeOfLastGridPointInDegrees"] = eccodes.codes_get( gid, "longitudeOfLastGridPointInDegrees" ) @@ -376,7 +374,7 @@ def test_area_vs_boundingbox_to_grib(self): match, report = _compare_grib_files(mars_grib_path, covjson_grib_path) print(report) - print(f"\nFiles preserved for inspection:") + print("\nFiles preserved for inspection:") print(f" MARS area GRIB: {mars_grib_path}") print(f" CovJSON→GRIB: {covjson_grib_path}") print(f" CoverageJSON: {covjson_path}") diff --git a/tests/test_grib_bounding_box.py b/tests/test_grib_bounding_box.py index 4c9df67..4b3c623 100644 --- a/tests/test_grib_bounding_box.py +++ b/tests/test_grib_bounding_box.py @@ -11,11 +11,11 @@ from covjsonkit.api import Covjsonkit # noqa: E402 - # --------------------------------------------------------------------------- # Fixtures # --------------------------------------------------------------------------- + @pytest.fixture def oper_covjson(): """Load the oper (no-ensemble) multipoint test fixture.""" @@ -85,11 +85,13 @@ def _read_grib_messages(filepath): # Step 6 stub: other decoders raise NotImplementedError # --------------------------------------------------------------------------- + class TestToGribStubs: """Verify that non-BoundingBox decoders raise NotImplementedError.""" def test_timeseries_raises(self): from covjsonkit.decoder.TimeSeries import TimeSeries + with pytest.raises(NotImplementedError): # We just need to call to_grib on the class — but we can't # instantiate without valid covjson, so test the method directly. @@ -97,6 +99,7 @@ def test_timeseries_raises(self): def test_vertical_profile_raises(self): from covjsonkit.decoder.VerticalProfile import VerticalProfile + with pytest.raises(NotImplementedError): VerticalProfile.to_grib(None) @@ -105,6 +108,7 @@ def test_vertical_profile_raises(self): # Core tests # --------------------------------------------------------------------------- + class TestToGribOper: """Tests using the oper (no-ensemble) fixture.""" From 928c16b477f158d239d992628614b37a7b8f545e Mon Sep 17 00:00:00 2001 From: awarde96 Date: Wed, 22 Jul 2026 15:53:06 +0200 Subject: [PATCH 4/4] Add prototype covjson2grib code using the mars2grib functionality from metkit --- covjsonkit/decoder/grib_backends/__init__.py | 22 +- .../grib_backends/mars2grib_backend.py | 197 ++++++- .../grib_backends/mars2grib_native_backend.py | 163 ++++++ tests/test_grib_bounding_box.py | 30 +- tests/test_grib_config_driven.py | 512 ++++++++++++++++++ tests/test_grib_mars_vs_bbox_backends.py | 302 +++++++++++ tests/test_grib_snap_out_comparison.py | 340 ++++++++++++ tests/test_grib_snap_out_polygon.py | 190 +++++++ tests/test_pymars2grib_smoke.py | 299 ++++++++++ 9 files changed, 2026 insertions(+), 29 deletions(-) create mode 100644 covjsonkit/decoder/grib_backends/mars2grib_native_backend.py create mode 100644 tests/test_grib_config_driven.py create mode 100644 tests/test_grib_mars_vs_bbox_backends.py create mode 100644 tests/test_grib_snap_out_comparison.py create mode 100644 tests/test_grib_snap_out_polygon.py create mode 100644 tests/test_pymars2grib_smoke.py diff --git a/covjsonkit/decoder/grib_backends/__init__.py b/covjsonkit/decoder/grib_backends/__init__.py index 90c37e1..2fd5125 100644 --- a/covjsonkit/decoder/grib_backends/__init__.py +++ b/covjsonkit/decoder/grib_backends/__init__.py @@ -1,8 +1,11 @@ """Swappable GRIB encoding backends. Use :func:`get_backend` to obtain the best available backend at runtime. -The factory tries ``pymars2grib`` first (preferred) and falls back to -``eccodes`` when it is not installed. +Available backends (in preference order for ``"auto"``): + + * ``mars2grib`` – hybrid: mars2grib for header, eccodes for section 3 + * ``mars2grib_native`` – pure mars2grib via compound gridSpec + skipSection3 + * ``eccodes`` – hand-rolled with eccodes only """ from .base import GribBackend # noqa: F401 @@ -12,10 +15,12 @@ def get_backend(preferred: str = "auto") -> GribBackend: """Return a ready-to-use GRIB encoding backend. Args: - preferred: One of ``"auto"``, ``"mars2grib"``, or ``"eccodes"``. - * ``"auto"`` – try mars2grib first, fall back to eccodes. - * ``"mars2grib"`` – require mars2grib (raises ImportError if - unavailable). + preferred: One of ``"auto"``, ``"mars2grib"``, ``"mars2grib_native"``, + or ``"eccodes"``. + * ``"auto"`` – try mars2grib (hybrid), fall back to eccodes. + * ``"mars2grib"`` – hybrid mars2grib + eccodes section 3. + * ``"mars2grib_native"`` – pure mars2grib via compound gridSpec + + ``skipSection3``. Requires ``ECCODES_ECKIT_GEO`` env var. * ``"eccodes"`` – use eccodes directly. Returns: @@ -24,6 +29,11 @@ def get_backend(preferred: str = "auto") -> GribBackend: Raises: ImportError: If the explicitly requested backend is not installed. """ + if preferred == "mars2grib_native": + from .mars2grib_native_backend import Mars2GribNativeBackend + + return Mars2GribNativeBackend() + if preferred in ("mars2grib", "auto"): try: from .mars2grib_backend import Mars2GribBackend diff --git a/covjsonkit/decoder/grib_backends/mars2grib_backend.py b/covjsonkit/decoder/grib_backends/mars2grib_backend.py index 6438292..42347c7 100644 --- a/covjsonkit/decoder/grib_backends/mars2grib_backend.py +++ b/covjsonkit/decoder/grib_backends/mars2grib_backend.py @@ -1,41 +1,204 @@ """GRIB encoding backend using pymars2grib (metkit's Python binding). -This is the preferred backend when available. pymars2grib handles all -MARS-to-GRIB key resolution internally, so this wrapper is intentionally -thin. +Hybrid approach: pymars2grib generates a full GRIB message with correct +header layout (sections 1, 4, 5) from MARS keys. Then eccodes overwrites +section 3 (grid geometry) with sub-area information and updates the data +values. + +NOTE: We do NOT use ``skipSection3``. In the current metkit+eccodes +build, ``skipSection3=True`` triggers the ``gridSpec`` string setter path, +and eccodes reports ``GridSpec::pack_string not available``. Instead we +let mars2grib emit a normal template and overwrite section 3 via eccodes +afterwards. pymars2grib is not yet available on PyPI — it must be built from metkit -source with pybind11 support. When it is not installed the factory in -``__init__.py`` falls back to the eccodes backend. +source with pybind11 support. When not installed, the factory in +``__init__.py`` falls back to the eccodes-only backend. + +Build instructions: see ``mars2grib-bundle/build.sh``. """ +import logging + from .base import GribBackend +logger = logging.getLogger(__name__) + +# MARS keys that pymars2grib expects as integers (long). +# The C++ layer uses strict type checking via eckit::LocalConfiguration. +_INT_MARS_KEYS = frozenset({"param", "step", "number", "levelist", "date", "time"}) + +# Ensemble stream/type combinations that require numberOfForecastsInEnsemble +# in the misc dict — mars2grib has no default deduction for this key. +_ENSEMBLE_TYPES = frozenset({"pf", "cf", "em", "es", "fcmean", "fcstdev"}) +_ENSEMBLE_STREAMS = frozenset({"enfo", "elda", "waef", "eefo", "efov", "efho"}) + + +def _coerce_mars_types(mars: dict) -> dict: + """Coerce MARS dict values to the types pymars2grib expects.""" + result = {} + for key, value in mars.items(): + if key in _INT_MARS_KEYS: + try: + result[key] = int(value) + except (ValueError, TypeError): + result[key] = value + else: + result[key] = value + return result + + +def _build_mars2grib_mars_dict(mars: dict, misc: dict) -> dict: + """Build the mars dict that pymars2grib expects. + + Adds required keys (origin, grid, packing) that mars2grib needs for + header layout resolution but our _build_mars_dict doesn't provide. + Also strips ``number`` for non-ensemble types — polytope-mars sets + ``number: 0`` on oper coverages, which incorrectly routes ``fc`` into + mars2grib's ensemble path and aborts on ``typeOfEnsembleForecast``. + """ + result = _coerce_mars_types(mars) + + # Drop `number` for non-ensemble types. polytope-mars includes + # `number: 0` on oper coverages; mars2grib treats its presence as a + # signal to encode as an ensemble member. + stream = str(result.get("stream", "")).lower() + typ = str(result.get("type", "")).lower() + is_ensemble = stream in _ENSEMBLE_STREAMS or typ in _ENSEMBLE_TYPES + if not is_ensemble and "number" in result: + result.pop("number") + + # origin — required by mars2grib + result.setdefault("origin", "ecmf") + + # packing — CCSDS matches operational MARS data + result.setdefault("packing", "ccsds") + + # grid — required by mars2grib; derive from misc gridType/N. + # This drives the template's section 3, which we overwrite anyway, + # but mars2grib needs *some* valid grid to route the encoder. + if "grid" not in result: + grid_type = misc.get("gridType", "reduced_gg") + n = misc.get("N", 1280) + if grid_type == "reduced_gg": + result["grid"] = f"O{n}" + elif grid_type == "regular_ll": + dx = misc.get("iDirectionIncrementInDegrees", 0.1) + dy = misc.get("jDirectionIncrementInDegrees", 0.1) + result["grid"] = f"{dy}/{dx}" + else: + result["grid"] = f"O{n}" + logger.warning("Unknown gridType '%s', defaulting to O%d", grid_type, n) + + return result + + +def _build_mars2grib_misc_dict(mars: dict, misc: dict) -> dict: + """Build the misc dict that pymars2grib expects. + + mars2grib silently ignores geometry keys (gridType, Nj, pl, area, …) + — those are only used by our eccodes overwrite step. But some + concepts (notably ensembles) need explicit deduction hints here. + """ + result = {} + + # Ensemble deduction: mars2grib has no default for numberOfForecastsInEnsemble + stream = str(mars.get("stream", "")).lower() + typ = str(mars.get("type", "")).lower() + if stream in _ENSEMBLE_STREAMS or typ in _ENSEMBLE_TYPES: + # Prefer explicit misc override; otherwise fall back to ECMWF ENS size + result["numberOfForecastsInEnsemble"] = int(misc.get("numberOfForecastsInEnsemble", 51)) + + return result + + +def _apply_section3_with_eccodes(template_bytes: bytes, misc: dict, values: list) -> bytes: + """Overwrite GRIB section 3 (grid definition) and data values via eccodes.""" + import eccodes + + gid = eccodes.codes_new_from_message(template_bytes) + + try: + grid_type = misc.get("gridType", "reduced_gg") + eccodes.codes_set(gid, "gridType", grid_type) + + if grid_type == "reduced_gg": + n = misc.get("N", 1280) + nj = misc.get("Nj", 1) + pl = misc.get("pl", [len(values)]) + + eccodes.codes_set_long(gid, "numberOfParallelsBetweenAPoleAndTheEquator", int(n)) + eccodes.codes_set_long(gid, "Nj", int(nj)) + eccodes.codes_set_long_array(gid, "pl", [int(p) for p in pl]) + + elif grid_type == "regular_ll": + for key in ("Ni", "Nj"): + if key in misc: + eccodes.codes_set_long(gid, key, int(misc[key])) + for key in ("iDirectionIncrementInDegrees", "jDirectionIncrementInDegrees"): + if key in misc: + eccodes.codes_set_double(gid, key, float(misc[key])) + + # Area bounds + area = misc.get("area") + if area and len(area) == 4: + eccodes.codes_set_double(gid, "latitudeOfFirstGridPointInDegrees", float(area[0])) + eccodes.codes_set_double(gid, "longitudeOfFirstGridPointInDegrees", float(area[1])) + eccodes.codes_set_double(gid, "latitudeOfLastGridPointInDegrees", float(area[2])) + eccodes.codes_set_double(gid, "longitudeOfLastGridPointInDegrees", float(area[3])) + + # Packing precision + bits_per_value = misc.get("bitsPerValue", 16) + eccodes.codes_set_long(gid, "bitsPerValue", int(bits_per_value)) + + # Set the actual data values + eccodes.codes_set_values(gid, [float(v) for v in values]) + + return eccodes.codes_get_message(gid) + + finally: + eccodes.codes_release(gid) + class Mars2GribBackend(GribBackend): - """Encode GRIB2 messages via pymars2grib.""" + """Encode GRIB2 messages via pymars2grib + eccodes hybrid. + + mars2grib resolves the GRIB header layout (sections 1, 4, 5) from + MARS keys. eccodes then overwrites section 3 (grid geometry) with + sub-area info and updates the data values. + """ def __init__(self): try: from pymars2grib import Mars2Grib except ImportError: raise ImportError( - "pymars2grib is not installed. " "Build metkit from source with pybind11 support to use this backend." + "pymars2grib is not installed. " + "Build metkit from source with pybind11 support to use this backend. " + "See mars2grib-bundle/build.sh." ) + # No skipSection3 — that path is broken in current builds. + # We let mars2grib produce a normal template and overwrite section 3. self._encoder = Mars2Grib() def encode_message(self, values: list, mars: dict, misc: dict) -> bytes: - """Encode a single GRIB2 message. + """Encode a single GRIB2 message using the hybrid approach. - Delegates entirely to pymars2grib which resolves GRIB header - layout from the MARS dictionary and injects the values. + 1. mars2grib generates a full GRIB message from MARS keys + 2. eccodes overwrites section 3 (grid geometry) with sub-area info + 3. eccodes updates the data values Args: - values: Field data values. - mars: MARS keys describing the field. - misc: Auxiliary metadata (grid geometry, packing hints, …). - - Returns: - Encoded GRIB2 message as bytes. + values: Field data values (sub-area, N→S/W→E ordered). + mars: MARS keys describing the field (class, stream, type, param, …). + misc: Grid geometry, encoding hints, ensemble deductions, … """ - return self._encoder.encode(values, mars, misc) + typed_mars = _build_mars2grib_mars_dict(mars, misc) + typed_misc = _build_mars2grib_misc_dict(mars, misc) + + # mars2grib requires values matching its template's grid size, but + # since we overwrite section 3 the template's grid size is thrown + # away. Just pass the real values — eccodes will re-set them. + template_bytes = self._encoder.encode(values, typed_mars, typed_misc) + + return _apply_section3_with_eccodes(template_bytes, misc, values) diff --git a/covjsonkit/decoder/grib_backends/mars2grib_native_backend.py b/covjsonkit/decoder/grib_backends/mars2grib_native_backend.py new file mode 100644 index 0000000..f70885c --- /dev/null +++ b/covjsonkit/decoder/grib_backends/mars2grib_native_backend.py @@ -0,0 +1,163 @@ +"""GRIB encoding backend using pymars2grib with **native sub-area support** +via compound gridSpec + ``skipSection3``. + +This is the "clean" mars2grib path recommended by the metkit developer and +documented in ``Mars2GribEncoding.md``. Unlike the hybrid backend +(``mars2grib_backend.py``), it does NOT touch section 3 with eccodes: + + 1. Build a compound gridSpec string: + mars["grid"] = "{grid: O1280, area: [N, W, S, E]}" + 2. Encode with ``Mars2Grib(options={"skipSection3": True})`` + 3. metkit forwards the string to eccodes' ``gridSpec`` accessor + (``GridSpec::pack_string``), eckit-geo builds the reduced-gg + sub-area natively, section 3 is written by eccodes with the + correct corners, ``numberOfDataPoints``, and full ``pl``. + +Runtime requirement: + ``ECCODES_ECKIT_GEO=2`` must be exported (or ``=1``). Without it, + ``GridSpec::pack_string`` returns ``GRIB_NOT_IMPLEMENTED`` — see + ``eccodes/geo/eckit.h`` for the ``EckitGeoLevel`` enum. Our + ``mars2grib-bundle/env.sh`` sets this correctly. + +Ensemble note: + Same as the hybrid backend — mars2grib has no default deduction for + ``numberOfForecastsInEnsemble``, so we inject it into ``misc`` when + the stream/type is ensemble. +""" + +import logging +import os + +from .base import GribBackend +from .mars2grib_backend import _ENSEMBLE_STREAMS, _ENSEMBLE_TYPES, _coerce_mars_types + +logger = logging.getLogger(__name__) + + +def _build_compound_grid_spec(misc: dict, values_len: int) -> str: + """Build the compound gridSpec string that eckit-geo can decode. + + Format:: + + {grid: , area: [N, W, S, E]} + + for reduced-gg grids, or:: + + {grid: [dx, dy], area: [N, W, S, E]} + + for regular_ll. The area is taken from ``misc["area"]`` if present, + otherwise from the coverage bounds computed by + ``BoundingBox._build_misc_dict``. + """ + grid_type = misc.get("gridType", "reduced_gg") + + if grid_type == "reduced_gg": + n = misc.get("N", 1280) + grid_str = f"O{int(n)}" + elif grid_type == "regular_ll": + dx = misc.get("iDirectionIncrementInDegrees", 0.1) + dy = misc.get("jDirectionIncrementInDegrees", 0.1) + grid_str = f"[{dx}, {dy}]" + else: + n = misc.get("N", 1280) + grid_str = f"O{int(n)}" + logger.warning("Unknown gridType '%s', defaulting to O%d", grid_type, n) + + area = misc.get("area") + if area and len(area) == 4: + return f"{{grid: {grid_str}, area: [{area[0]}, {area[1]}, {area[2]}, {area[3]}]}}" + + # No area — full-globe (uncommon for our use case) + logger.warning("No area in misc dict; encoding as full globe (values_len=%d)", values_len) + return f"{{grid: {grid_str}}}" + + +def _build_mars_dict(mars: dict, misc: dict) -> dict: + """Build the mars dict for the compound-gridSpec path. + + Same as the hybrid backend, but ``grid`` is a compound spec, and we + also strip ``number`` for non-ensemble types (polytope-mars sets + ``number: 0`` on oper coverages, which otherwise routes ``fc`` into + ensemble encoding and aborts). + """ + result = _coerce_mars_types(mars) + + # Strip spurious `number: 0` on oper coverages + stream = str(result.get("stream", "")).lower() + typ = str(result.get("type", "")).lower() + is_ensemble = stream in _ENSEMBLE_STREAMS or typ in _ENSEMBLE_TYPES + if not is_ensemble and "number" in result: + result.pop("number") + + result.setdefault("origin", "ecmf") + result.setdefault("packing", "ccsds") + + # Compound gridSpec — verbatim to eccodes' gridSpec accessor + result["grid"] = _build_compound_grid_spec(misc, values_len=0) + return result + + +def _build_misc_dict(mars: dict, misc: dict, values_len: int) -> dict: + """Build the misc dict — mars2grib only cares about a few keys here. + + Geometry keys (gridType, Nj, pl, area…) are ignored by mars2grib — + they go into the compound gridSpec string on ``mars["grid"]`` instead. + """ + result = {} + + # Ensemble deduction + stream = str(mars.get("stream", "")).lower() + typ = str(mars.get("type", "")).lower() + if stream in _ENSEMBLE_STREAMS or typ in _ENSEMBLE_TYPES: + result["numberOfForecastsInEnsemble"] = int(misc.get("numberOfForecastsInEnsemble", 51)) + + # Packing precision + if "bitsPerValue" in misc: + result["bitsPerValue"] = int(misc["bitsPerValue"]) + + return result + + +class Mars2GribNativeBackend(GribBackend): + """Encode GRIB2 messages natively via ``Mars2Grib(skipSection3=True)``. + + Uses compound gridSpec on ``mars["grid"]`` so eckit-geo builds the + sub-area natively and eccodes writes section 3 correctly. No manual + section-3 overwrite required. + + Trade-off vs the hybrid backend: + * ``+`` No hand-rolled ``Nj``/``pl`` computation — eckit-geo owns it. + * ``+`` Corner lat/lons match MARS's exact snap-out convention. + * ``-`` Requires ``ECCODES_ECKIT_GEO`` env var; without it, aborts. + """ + + def __init__(self): + try: + from pymars2grib import Mars2Grib + except ImportError: + raise ImportError( + "pymars2grib is not installed. " + "Build metkit from source with pybind11 support to use this backend. " + "See mars2grib-bundle/build.sh." + ) + + # Runtime env-var check — GridSpec::pack_string is gated on this. + if os.environ.get("ECCODES_ECKIT_GEO", "0") == "0": + logger.warning( + "ECCODES_ECKIT_GEO is not set — mars2grib_native will abort " + "when encoding. Run `source mars2grib-bundle/env.sh` first." + ) + + self._encoder = Mars2Grib(options={"skipSection3": True}) + + def encode_message(self, values: list, mars: dict, misc: dict) -> bytes: + """Encode a single GRIB2 message. + + Args: + values: Field data values (sub-area, N→S/W→E ordered). + mars: MARS keys. + misc: Grid geometry (used to build the compound gridSpec). + """ + typed_mars = _build_mars_dict(mars, misc) + typed_misc = _build_misc_dict(mars, misc, values_len=len(values)) + return self._encoder.encode([float(v) for v in values], typed_mars, typed_misc) diff --git a/tests/test_grib_bounding_box.py b/tests/test_grib_bounding_box.py index 4b3c623..2b64334 100644 --- a/tests/test_grib_bounding_box.py +++ b/tests/test_grib_bounding_box.py @@ -323,13 +323,26 @@ def test_no_grid_metadata_uses_defaults(self, oper_covjson): class TestBackendSelection: """Test backend factory behaviour.""" - def test_auto_selects_eccodes(self): - """With pymars2grib unavailable, auto should fall back to eccodes.""" + def test_auto_selects_available_backend(self): + """Auto picks mars2grib if importable, otherwise eccodes.""" from covjsonkit.decoder.grib_backends import get_backend from covjsonkit.decoder.grib_backends.eccodes_backend import EccodesBackend backend = get_backend("auto") - assert isinstance(backend, EccodesBackend) + try: + from covjsonkit.decoder.grib_backends.mars2grib_backend import ( + Mars2GribBackend, + ) + + # pymars2grib available → auto should give mars2grib + try: + import pymars2grib # noqa: F401 + + assert isinstance(backend, Mars2GribBackend) + except ImportError: + assert isinstance(backend, EccodesBackend) + except ImportError: + assert isinstance(backend, EccodesBackend) def test_explicit_eccodes(self): """Explicitly requesting eccodes should work.""" @@ -340,11 +353,16 @@ def test_explicit_eccodes(self): assert isinstance(backend, EccodesBackend) def test_mars2grib_raises_when_unavailable(self): - """Explicitly requesting mars2grib should raise if not installed.""" + """Explicitly requesting mars2grib should raise if pymars2grib absent.""" from covjsonkit.decoder.grib_backends import get_backend - with pytest.raises(ImportError, match="pymars2grib"): - get_backend("mars2grib") + try: + import pymars2grib # noqa: F401 + + pytest.skip("pymars2grib is installed — this test only runs without it") + except ImportError: + with pytest.raises(ImportError, match="pymars2grib"): + get_backend("mars2grib") class TestBuildHelpers: diff --git a/tests/test_grib_config_driven.py b/tests/test_grib_config_driven.py new file mode 100644 index 0000000..7230222 --- /dev/null +++ b/tests/test_grib_config_driven.py @@ -0,0 +1,512 @@ +"""Integration test: Config-driven grid metadata → to_grib() vs MARS area. + +This test simulates the full polytope-mars pipeline: + 1. Retrieve CoverageJSON via polytope BoundingBox feature request + 2. Simulate polytope-mars _get_grid_metadata() by extracting grid info + from a config (as polytope-server would provide) + 3. Inject mars:grid into each coverage (as the encoder would do) + 4. Convert enriched CoverageJSON → GRIB via to_grib() + 5. Compare against MARS area retrieve (reference GRIB) + +This validates the full chain: + polytope-server config → _get_grid_metadata() → encoder mars:grid + → decoder reads mars:grid → to_grib() produces correct output + +Requirements: + - Valid polytope credentials (~/.polytopeapirc) + - Network access to polytope.ecmwf.int + - eccodes Python package + +Run: + RUN_INTEGRATION_TESTS=1 pytest tests/test_grib_config_driven.py -v -s +""" + +import json +import os +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required") + +_skip_integration = pytest.mark.skipif( + os.environ.get("RUN_INTEGRATION_TESTS", "0") != "1", + reason="Integration test — set RUN_INTEGRATION_TESTS=1 to run", +) + +# --------------------------------------------------------------------------- +# Configuration +# --------------------------------------------------------------------------- + +AREA = [0.18, 0.0, 0.03, 0.15] # N/W/S/E + +MARS_REQUEST = { + "class": "od", + "stream": "oper", + "type": "fc", + "levtype": "sfc", + "param": "2t/10v", + "step": "0", + "date": "-1", + "time": "0000", +} + +COLLECTION = "ecmwf-mars" + +# Simulated polytope-server axis_config (matches polytope-od production config) +AXIS_CONFIG_OCTAHEDRAL = { + "axis_name": "values", + "transformations": [ + { + "name": "mapper", + "type": "octahedral", + "resolution": 1280, + "axes": ["latitude", "longitude"], + "local": None, + } + ], +} + +# EFAS-style regular_ll config for testing alternate grid types +AXIS_CONFIG_REGULAR_LL = { + "axis_name": "values", + "transformations": [ + { + "name": "mapper", + "type": "local_regular", + "resolution": [2969, 4529], + "axes": ["latitude", "longitude"], + "local": [22.76, 72.24, -25.24, 50.24], + } + ], +} + + +# --------------------------------------------------------------------------- +# Simulate polytope-mars _get_grid_metadata() +# --------------------------------------------------------------------------- + + +def _get_grid_metadata_from_config(axis_config: dict) -> dict: + """Replicate polytope-mars _get_grid_metadata() logic from a config dict. + + This is a standalone version of PolytopeMars._get_grid_metadata() that + operates on a raw axis_config dict rather than requiring a full + PolytopeMarsConfig instance. + """ + for transform in axis_config.get("transformations", []): + if transform.get("name") != "mapper": + continue + + mapper_type = transform.get("type") + resolution = transform.get("resolution") + + if mapper_type == "octahedral": + return {"gridType": "reduced_gg", "N": int(resolution)} + elif mapper_type == "local_regular": + meta = {"gridType": "regular_ll"} + if isinstance(resolution, list) and len(resolution) == 2: + meta["Nj"] = int(resolution[0]) + meta["Ni"] = int(resolution[1]) + local = transform.get("local") + if local and isinstance(local, list) and len(local) == 4: + meta["area"] = [float(v) for v in local] + return meta + else: + return {"mapperType": mapper_type} + + return {} + + +def _inject_grid_metadata(covjson: dict, grid_metadata: dict) -> dict: + """Inject mars:grid into each coverage (as the encoder would do).""" + if not grid_metadata: + return covjson + + for coverage in covjson.get("coverages", []): + coverage["mars:grid"] = grid_metadata + + return covjson + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _read_all_grib_messages(filepath): + """Read all messages from a GRIB file.""" + messages = [] + with open(filepath, "rb") as f: + while True: + gid = eccodes.codes_grib_new_from_file(f) + if gid is None: + break + try: + msg = { + "paramId": eccodes.codes_get(gid, "paramId"), + "shortName": eccodes.codes_get(gid, "shortName"), + "dataDate": eccodes.codes_get(gid, "dataDate"), + "dataTime": eccodes.codes_get(gid, "dataTime"), + "stepRange": eccodes.codes_get(gid, "stepRange"), + "gridType": eccodes.codes_get(gid, "gridType"), + "numberOfDataPoints": eccodes.codes_get(gid, "numberOfDataPoints"), + "bitsPerValue": eccodes.codes_get(gid, "bitsPerValue"), + "values": eccodes.codes_get_values(gid).tolist(), + } + try: + msg["N"] = eccodes.codes_get(gid, "N") + except Exception: + pass + try: + msg["Nj"] = eccodes.codes_get(gid, "Nj") + except Exception: + pass + try: + msg["Ni"] = eccodes.codes_get(gid, "Ni") + except Exception: + pass + try: + msg["latitudeOfFirstGridPointInDegrees"] = eccodes.codes_get( + gid, "latitudeOfFirstGridPointInDegrees" + ) + msg["longitudeOfFirstGridPointInDegrees"] = eccodes.codes_get( + gid, "longitudeOfFirstGridPointInDegrees" + ) + msg["latitudeOfLastGridPointInDegrees"] = eccodes.codes_get(gid, "latitudeOfLastGridPointInDegrees") + msg["longitudeOfLastGridPointInDegrees"] = eccodes.codes_get( + gid, "longitudeOfLastGridPointInDegrees" + ) + except Exception: + pass + messages.append(msg) + finally: + eccodes.codes_release(gid) + return messages + + +def _retrieve_mars_area(output_path): + """Retrieve reference GRIB via polytope mars-od with area keyword.""" + from polytope.api import Client + + client = Client() + request = {**MARS_REQUEST, "area": f"{AREA[0]}/{AREA[1]}/{AREA[2]}/{AREA[3]}"} + + print(f"\n{'='*60}") + print("REFERENCE: Polytope mars-od retrieve with area keyword") + print(f"{'='*60}") + print(f"Request: {request}") + + client.retrieve(COLLECTION, request, output_file=output_path, pointer=False) + + file_size = os.path.getsize(output_path) + print(f"Retrieved {file_size} bytes") + return output_path + + +def _retrieve_boundingbox_covjson(): + """Retrieve CoverageJSON via polytope BoundingBox feature.""" + from polytope.api import Client + + client = Client() + request = {**MARS_REQUEST} + request["feature"] = { + "type": "boundingbox", + "points": [[AREA[0], AREA[1]], [AREA[2], AREA[3]]], + } + + print(f"\n{'='*60}") + print("FEATURE: Polytope BoundingBox → CoverageJSON") + print(f"{'='*60}") + print(f"Request: {request}") + + with tempfile.NamedTemporaryFile(suffix=".covjson", delete=False, mode="w") as tmp: + covjson_path = tmp.name + + try: + client.retrieve(COLLECTION, request, output_file=covjson_path) + with open(covjson_path) as f: + covjson = json.load(f) + + print(f"Coverages: {len(covjson.get('coverages', []))}") + return covjson + finally: + if os.path.exists(covjson_path): + os.unlink(covjson_path) + + +def _convert_to_grib(covjson, output_path): + """Convert CoverageJSON to GRIB via to_grib().""" + from covjsonkit.api import Covjsonkit + + decoder = Covjsonkit().decode(covjson) + decoder.to_grib(output_path, backend="eccodes") + + file_size = os.path.getsize(output_path) + print(f"Wrote {file_size} bytes GRIB") + return output_path + + +def _compare_messages(mars_msgs, covjson_msgs, tolerance=0.001): + """Compare two sets of GRIB messages. Returns (all_match, report_lines).""" + report = [] + all_match = True + + report.append(f"MARS messages: {len(mars_msgs)}, CovJSON messages: {len(covjson_msgs)}") + + if len(mars_msgs) != len(covjson_msgs): + report.append(" ⚠️ Message count mismatch") + all_match = False + + for mars_msg in mars_msgs: + key = (mars_msg["shortName"], mars_msg["stepRange"]) + matching = [m for m in covjson_msgs if (m["shortName"], m["stepRange"]) == key] + + if not matching: + report.append(f" ❌ No match for {key}") + all_match = False + continue + + covjson_msg = matching[0] + report.append(f"\n Field: {key[0]} step={key[1]}") + + # Point count + if mars_msg["numberOfDataPoints"] != covjson_msg["numberOfDataPoints"]: + report.append( + f" ⚠️ Points: MARS={mars_msg['numberOfDataPoints']} " f"CovJSON={covjson_msg['numberOfDataPoints']}" + ) + all_match = False + else: + report.append(f" ✓ Points: {mars_msg['numberOfDataPoints']}") + + # Grid type + if mars_msg["gridType"] != covjson_msg["gridType"]: + report.append(f" ⚠️ Grid: MARS={mars_msg['gridType']} CovJSON={covjson_msg['gridType']}") + all_match = False + else: + report.append(f" ✓ Grid: {mars_msg['gridType']}") + + # N value (for reduced_gg) + if "N" in mars_msg and "N" in covjson_msg: + if mars_msg["N"] != covjson_msg["N"]: + report.append(f" ⚠️ N: MARS={mars_msg['N']} CovJSON={covjson_msg['N']}") + all_match = False + else: + report.append(f" ✓ N: {mars_msg['N']}") + + # Values + mars_vals = mars_msg["values"] + covjson_vals = covjson_msg["values"] + + if len(mars_vals) != len(covjson_vals): + report.append(f" ⚠️ Value count: {len(mars_vals)} vs {len(covjson_vals)}") + all_match = False + else: + diffs = [abs(a - b) for a, b in zip(mars_vals, covjson_vals)] + max_diff = max(diffs) if diffs else 0 + mean_diff = sum(diffs) / len(diffs) if diffs else 0 + n_exact = sum(1 for d in diffs if d == 0.0) + + report.append( + f" Values: max_diff={max_diff:.2e}, mean_diff={mean_diff:.2e}, exact={n_exact}/{len(mars_vals)}" + ) + + if max_diff > tolerance: + report.append(f" ⚠️ Exceeds tolerance {tolerance}") + all_match = False + else: + report.append(" ✓ Within tolerance") + + # Area bounds + for bound_key in [ + "latitudeOfFirstGridPointInDegrees", + "latitudeOfLastGridPointInDegrees", + "longitudeOfFirstGridPointInDegrees", + "longitudeOfLastGridPointInDegrees", + ]: + if bound_key in mars_msg and bound_key in covjson_msg: + diff = abs(mars_msg[bound_key] - covjson_msg[bound_key]) + status = "✓" if diff < 0.001 else "⚠️" + short_key = bound_key.replace("InDegrees", "").replace("OfGridPoint", "") + report.append(f" {status} {short_key}: {mars_msg[bound_key]:.6f} vs {covjson_msg[bound_key]:.6f}") + if diff >= 0.001: + all_match = False + + return all_match, report + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + + +@_skip_integration +class TestConfigDrivenGridToGrib: + """Full pipeline: config → grid_metadata → mars:grid → to_grib() vs MARS. + + Run: RUN_INTEGRATION_TESTS=1 pytest tests/test_grib_config_driven.py -v -s + """ + + def test_octahedral_config_vs_mars_area(self): + """Octahedral N=1280 config (polytope-od production) → GRIB matches MARS area. + + Pipeline: + 1. _get_grid_metadata_from_config(octahedral, resolution=1280) + → {"gridType": "reduced_gg", "N": 1280} + 2. Retrieve CovJSON via feature request + 3. Inject mars:grid (simulating encoder) + 4. to_grib() → GRIB + 5. Compare vs MARS area retrieve + """ + # Step 1: Extract grid metadata from config (as polytope-mars would) + grid_metadata = _get_grid_metadata_from_config(AXIS_CONFIG_OCTAHEDRAL) + print(f"\nGrid metadata from config: {grid_metadata}") + assert grid_metadata == {"gridType": "reduced_gg", "N": 1280} + + # Step 2: Retrieve reference GRIB from MARS + mars_grib = tempfile.NamedTemporaryFile(suffix="_mars.grib", delete=False) + mars_grib.close() + + # Step 3: Retrieve CoverageJSON via feature request + covjson_grib = tempfile.NamedTemporaryFile(suffix="_config.grib", delete=False) + covjson_grib.close() + + try: + _retrieve_mars_area(mars_grib.name) + covjson = _retrieve_boundingbox_covjson() + + # Step 4: Inject mars:grid (simulating what encoder does with grid_metadata) + print(f"\n{'='*60}") + print("INJECT: Adding mars:grid from config to each coverage") + print(f"{'='*60}") + covjson = _inject_grid_metadata(covjson, grid_metadata) + + # Verify injection + for i, cov in enumerate(covjson.get("coverages", [])): + assert "mars:grid" in cov, f"Coverage {i} missing mars:grid" + assert cov["mars:grid"]["gridType"] == "reduced_gg" + assert cov["mars:grid"]["N"] == 1280 + print(f" Coverage {i}: mars:grid = {cov['mars:grid']}") + + # Step 5: Convert to GRIB + print(f"\n{'='*60}") + print("CONVERT: CoverageJSON (with mars:grid) → GRIB") + print(f"{'='*60}") + _convert_to_grib(covjson, covjson_grib.name) + + # Step 6: Compare + print(f"\n{'='*60}") + print("COMPARE: MARS area vs config-driven to_grib()") + print(f"{'='*60}") + mars_msgs = _read_all_grib_messages(mars_grib.name) + covjson_msgs = _read_all_grib_messages(covjson_grib.name) + + match, report = _compare_messages(mars_msgs, covjson_msgs) + print("\n".join(report)) + + assert match, "Config-driven GRIB does not match MARS area output" + + finally: + for path in [mars_grib.name, covjson_grib.name]: + if os.path.exists(path): + os.unlink(path) + + def test_without_config_uses_defaults_still_matches(self): + """Without mars:grid (no config), fallback defaults should still produce matching values. + + This confirms backward compatibility — even without the config pipeline, + to_grib() uses _OPER_GRID_DEFAULTS and still produces correct values. + """ + mars_grib = tempfile.NamedTemporaryFile(suffix="_mars.grib", delete=False) + mars_grib.close() + covjson_grib = tempfile.NamedTemporaryFile(suffix="_noconfig.grib", delete=False) + covjson_grib.close() + + try: + _retrieve_mars_area(mars_grib.name) + covjson = _retrieve_boundingbox_covjson() + + # No mars:grid injection — uses defaults + assert "mars:grid" not in covjson.get("coverages", [{}])[0] + + print(f"\n{'='*60}") + print("CONVERT: CoverageJSON (no mars:grid, using defaults) → GRIB") + print(f"{'='*60}") + _convert_to_grib(covjson, covjson_grib.name) + + mars_msgs = _read_all_grib_messages(mars_grib.name) + covjson_msgs = _read_all_grib_messages(covjson_grib.name) + + match, report = _compare_messages(mars_msgs, covjson_msgs) + print("\n".join(report)) + + assert match, "Default-driven GRIB does not match MARS area output" + + finally: + for path in [mars_grib.name, covjson_grib.name]: + if os.path.exists(path): + os.unlink(path) + + +class TestGridMetadataExtraction: + """Unit tests for _get_grid_metadata_from_config() — no network needed.""" + + def test_octahedral_config(self): + """Octahedral mapper → reduced_gg with N.""" + result = _get_grid_metadata_from_config(AXIS_CONFIG_OCTAHEDRAL) + assert result == {"gridType": "reduced_gg", "N": 1280} + + def test_regular_ll_config(self): + """local_regular mapper → regular_ll with Ni, Nj, area.""" + result = _get_grid_metadata_from_config(AXIS_CONFIG_REGULAR_LL) + assert result == { + "gridType": "regular_ll", + "Nj": 2969, + "Ni": 4529, + "area": [22.76, 72.24, -25.24, 50.24], + } + + def test_no_mapper_transform(self): + """Config without mapper transform → empty dict.""" + config = { + "axis_name": "values", + "transformations": [{"name": "type_change", "type": "int"}], + } + result = _get_grid_metadata_from_config(config) + assert result == {} + + def test_empty_transformations(self): + """Config with no transformations → empty dict.""" + config = {"axis_name": "values", "transformations": []} + result = _get_grid_metadata_from_config(config) + assert result == {} + + def test_unknown_mapper_type(self): + """Unknown mapper type → returns mapperType key.""" + config = { + "axis_name": "values", + "transformations": [{"name": "mapper", "type": "healpix", "resolution": 512}], + } + result = _get_grid_metadata_from_config(config) + assert result == {"mapperType": "healpix"} + + def test_injection_adds_mars_grid(self): + """_inject_grid_metadata adds mars:grid to every coverage.""" + covjson = { + "coverages": [ + {"mars:metadata": {"class": "od"}}, + {"mars:metadata": {"class": "od"}}, + ] + } + grid_meta = {"gridType": "reduced_gg", "N": 1280} + result = _inject_grid_metadata(covjson, grid_meta) + + for cov in result["coverages"]: + assert cov["mars:grid"] == grid_meta + + def test_injection_empty_metadata_no_op(self): + """Empty grid_metadata → no mars:grid injected.""" + covjson = {"coverages": [{"mars:metadata": {"class": "od"}}]} + result = _inject_grid_metadata(covjson, {}) + + assert "mars:grid" not in result["coverages"][0] diff --git a/tests/test_grib_mars_vs_bbox_backends.py b/tests/test_grib_mars_vs_bbox_backends.py new file mode 100644 index 0000000..c21a9e8 --- /dev/null +++ b/tests/test_grib_mars_vs_bbox_backends.py @@ -0,0 +1,302 @@ +"""Integration test: MARS `area` retrieve vs polytope BoundingBox→to_grib(), +compared for **both** GRIB backends (eccodes and mars2grib). + +For each backend, this test: + 1. Prints the exact MARS request and polytope BoundingBox request + 2. Retrieves the reference GRIB via `polytope retrieve` with `area` + 3. Retrieves CoverageJSON via `polytope retrieve` with `feature=boundingbox` + 4. Converts CovJSON → GRIB via `Covjsonkit.decode(...).to_grib(backend=...)` + 5. Compares the two GRIB files (metadata + values) and reports diffs + 6. Asserts they match within tolerance + +Prerequisites: + - ~/.polytopeapirc with valid credentials + - Network access to polytope.ecmwf.int + - eccodes (Python) always; pymars2grib for the mars2grib backend + - For mars2grib: `source mars2grib-bundle/env.sh` before invoking pytest + +Run: + export RUN_INTEGRATION_TESTS=1 + source mars2grib-bundle/env.sh # only needed for the mars2grib parametrisation + .venv311/bin/pytest tests/test_grib_mars_vs_bbox_backends.py -v -s +""" + +import json +import os +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required") + +pytestmark = pytest.mark.skipif( + os.environ.get("RUN_INTEGRATION_TESTS", "0") != "1", + reason="Integration test — set RUN_INTEGRATION_TESTS=1 to run", +) + +# --------------------------------------------------------------------------- +# Request configuration — identical MARS keys used across both paths. +# Small equatorial box keeps GRIB messages tiny and easy to diff. +# --------------------------------------------------------------------------- + +AREA = [0.18, 0.0, 0.03, 0.15] # N / W / S / E (degrees) + +MARS_KEYS = { + "class": "od", + "stream": "oper", + "type": "fc", + "levtype": "sfc", + "param": "2t/10v", + "step": "0", + "date": "-1", # yesterday + "time": "0000", +} + +COLLECTION = "ecmwf-mars" + +# Value-comparison tolerance. GRIB2 with bitsPerValue=16 quantises to +# ~2e-5 K; both backends use identical section-3+7 encoding paths so +# they must land in the same bucket. +VALUE_TOLERANCE = 1e-4 + + +# --------------------------------------------------------------------------- +# GRIB helpers +# --------------------------------------------------------------------------- + + +def _read_all_grib_messages(filepath): + messages = [] + with open(filepath, "rb") as f: + while True: + gid = eccodes.codes_grib_new_from_file(f) + if gid is None: + break + try: + msg = { + "paramId": eccodes.codes_get(gid, "paramId"), + "shortName": eccodes.codes_get(gid, "shortName"), + "dataDate": eccodes.codes_get(gid, "dataDate"), + "dataTime": eccodes.codes_get(gid, "dataTime"), + "stepRange": eccodes.codes_get(gid, "stepRange"), + "gridType": eccodes.codes_get(gid, "gridType"), + "numberOfDataPoints": eccodes.codes_get(gid, "numberOfDataPoints"), + "values": eccodes.codes_get_values(gid).tolist(), + } + for k in ( + "latitudeOfFirstGridPointInDegrees", + "longitudeOfFirstGridPointInDegrees", + "latitudeOfLastGridPointInDegrees", + "longitudeOfLastGridPointInDegrees", + ): + try: + msg[k] = eccodes.codes_get(gid, k) + except Exception: + pass + messages.append(msg) + finally: + eccodes.codes_release(gid) + return messages + + +# --------------------------------------------------------------------------- +# Polytope retrieval helpers +# --------------------------------------------------------------------------- + + +def _mars_area_request(): + """Exact dict sent to polytope for the MARS `area` path.""" + return {**MARS_KEYS, "area": f"{AREA[0]}/{AREA[1]}/{AREA[2]}/{AREA[3]}"} + + +def _boundingbox_request(): + """Exact dict sent to polytope for the BoundingBox feature path.""" + return { + **MARS_KEYS, + "feature": { + "type": "boundingbox", + "points": [[AREA[0], AREA[1]], [AREA[2], AREA[3]]], + }, + } + + +def _retrieve_mars_area(output_path): + from polytope.api import Client + + request = _mars_area_request() + print("\n" + "=" * 70) + print("PATH 1 — MARS `area` retrieve") + print("=" * 70) + print(f"Collection: {COLLECTION}") + print(f"Request: {json.dumps(request, indent=2)}") + print(f"Output: {output_path}") + + Client().retrieve(COLLECTION, request, output_file=output_path, pointer=False) + print(f"Retrieved {os.path.getsize(output_path)} bytes") + + +def _retrieve_boundingbox_covjson(): + from polytope.api import Client + + request = _boundingbox_request() + print("\n" + "=" * 70) + print("PATH 2 — polytope BoundingBox → CoverageJSON → to_grib()") + print("=" * 70) + print(f"Collection: {COLLECTION}") + print(f"Request: {json.dumps(request, indent=2)}") + + with tempfile.NamedTemporaryFile(suffix=".covjson", delete=False, mode="w") as tmp: + cj_path = tmp.name + try: + Client().retrieve(COLLECTION, request, output_file=cj_path) + print(f"Retrieved {os.path.getsize(cj_path)} bytes CoverageJSON") + with open(cj_path) as f: + return json.load(f) + finally: + if os.path.exists(cj_path): + os.unlink(cj_path) + + +def _covjson_to_grib(covjson, output_path, backend): + from covjsonkit.api import Covjsonkit + + print(f"\nConverting CoverageJSON → GRIB (backend={backend!r})") + decoder = Covjsonkit().decode(covjson) + decoder.to_grib(output_path, backend=backend) + print(f"Wrote {os.path.getsize(output_path)} bytes to {output_path}") + + +# --------------------------------------------------------------------------- +# Comparison +# --------------------------------------------------------------------------- + + +def _compare_grib_files(mars_path, ours_path, backend): + mars = _read_all_grib_messages(mars_path) + ours = _read_all_grib_messages(ours_path) + + lines = [] + lines.append("\n" + "=" * 70) + lines.append(f"COMPARISON — backend={backend!r}") + lines.append("=" * 70) + lines.append(f"MARS messages: {len(mars)}") + lines.append(f"CovJSON→GRIB msgs: {len(ours)}") + + all_match = True + if len(mars) != len(ours): + lines.append(f"❌ message count mismatch: {len(mars)} vs {len(ours)}") + all_match = False + + mars_params = sorted({m["shortName"] for m in mars}) + ours_params = sorted({m["shortName"] for m in ours}) + lines.append(f"MARS params: {mars_params}") + lines.append(f"Ours params: {ours_params}") + if mars_params != ours_params: + lines.append("❌ parameter mismatch") + all_match = False + + for mm in mars: + key = (mm["shortName"], mm["stepRange"]) + cand = [m for m in ours if (m["shortName"], m["stepRange"]) == key] + if not cand: + lines.append(f"\n❌ no match for {key}") + all_match = False + continue + om = cand[0] + lines.append(f"\n--- {key[0]} step={key[1]} ---") + + if mm["numberOfDataPoints"] != om["numberOfDataPoints"]: + lines.append(f" ❌ point count: MARS={mm['numberOfDataPoints']} vs ours={om['numberOfDataPoints']}") + all_match = False + continue + lines.append(f" ✓ point count: {mm['numberOfDataPoints']}") + + if mm["gridType"] != om["gridType"]: + lines.append(f" ⚠ gridType differs: MARS={mm['gridType']} vs ours={om['gridType']}") + + diffs = [abs(a - b) for a, b in zip(mm["values"], om["values"])] + max_d = max(diffs) if diffs else 0.0 + mean_d = (sum(diffs) / len(diffs)) if diffs else 0.0 + n_exact = sum(1 for d in diffs if d == 0.0) + lines.append(f" values: max_diff={max_d:.3e} mean={mean_d:.3e} " f"exact={n_exact}/{len(diffs)}") + if max_d > VALUE_TOLERANCE: + lines.append(f" ❌ max_diff {max_d:.3e} exceeds tolerance {VALUE_TOLERANCE:.0e}") + all_match = False + + for bk in ( + "latitudeOfFirstGridPointInDegrees", + "longitudeOfFirstGridPointInDegrees", + "latitudeOfLastGridPointInDegrees", + "longitudeOfLastGridPointInDegrees", + ): + if bk in mm and bk in om: + d = abs(mm[bk] - om[bk]) + mark = "✓" if d < 1e-3 else "⚠" + lines.append(f" {mark} {bk}: MARS={mm[bk]:.6f} ours={om[bk]:.6f} (Δ={d:.2e})") + + lines.append("\n" + ("✅ MATCH" if all_match else "❌ DIFFERENCES") + f" (backend={backend!r})") + return all_match, "\n".join(lines) + + +# --------------------------------------------------------------------------- +# Parametrised test — one row per backend +# --------------------------------------------------------------------------- + + +def _backend_available(name): + if name == "eccodes": + return True + if name in ("mars2grib", "mars2grib_native"): + try: + import pymars2grib # noqa: F401 + + return True + except ImportError: + return False + return False + + +@pytest.mark.parametrize("backend", ["eccodes", "mars2grib", "mars2grib_native"]) +def test_mars_area_vs_boundingbox_to_grib(backend): + """MARS `area` output must equal polytope BoundingBox → to_grib(backend).""" + if not _backend_available(backend): + pytest.skip(f"backend '{backend}' unavailable (source mars2grib-bundle/env.sh)") + + mars_f = tempfile.NamedTemporaryFile(suffix=f"_mars_{backend}.grib", delete=False) + ours_f = tempfile.NamedTemporaryFile(suffix=f"_ours_{backend}.grib", delete=False) + mars_f.close() + ours_f.close() + + try: + _retrieve_mars_area(mars_f.name) + covjson = _retrieve_boundingbox_covjson() + _covjson_to_grib(covjson, ours_f.name, backend=backend) + + match, report = _compare_grib_files(mars_f.name, ours_f.name, backend) + print(report) + assert match, f"backend={backend}: GRIB files differ — see report above" + finally: + for p in (mars_f.name, ours_f.name): + if os.path.exists(p): + os.unlink(p) + + +if __name__ == "__main__": + """Run directly (bypasses pytest for quick manual iteration).""" + print("MARS area vs BoundingBox→to_grib comparison") + print(f"AREA: N={AREA[0]} W={AREA[1]} S={AREA[2]} E={AREA[3]}") + + os.environ["RUN_INTEGRATION_TESTS"] = "1" + mars_path = "/tmp/mars_area_ref.grib" + _retrieve_mars_area(mars_path) + covjson = _retrieve_boundingbox_covjson() + + for backend in ("eccodes", "mars2grib", "mars2grib_native"): + if not _backend_available(backend): + print(f"\n[skip] backend={backend} unavailable") + continue + ours_path = f"/tmp/covjson_{backend}.grib" + _covjson_to_grib(covjson, ours_path, backend=backend) + match, report = _compare_grib_files(mars_path, ours_path, backend) + print(report) + print(f"\n→ files preserved: {mars_path}, {ours_path}") diff --git a/tests/test_grib_snap_out_comparison.py b/tests/test_grib_snap_out_comparison.py new file mode 100644 index 0000000..0addad5 --- /dev/null +++ b/tests/test_grib_snap_out_comparison.py @@ -0,0 +1,340 @@ +"""Snap-out comparison: polytope BoundingBox vs MARS `area`. + +Both paths take an N/W/S/E box and return a subset of the ECMWF reduced-gg +grid. They should agree on WHICH grid points fall inside the box, but they +use different snap-out conventions: + + * MARS `area` — eckit-geo `BoundingBox` crop applied server-side, then + eccodes writes the exact first/last lat/lon corners it kept. + * Polytope BoundingBox — floor-truncated coordinate walk over the + underlying datacube axes; returns the coverage's composite coordinates. + +This test characterises the divergence across representative boxes so we +can detect regressions when either side changes its snap-out logic. + +For each bbox: + 1. Retrieve MARS area GRIB → extract grid points via eccodes iterator + 2. Retrieve polytope BoundingBox → extract composite coordinates + 3. Compare: point count, per-point max lat/lon Δ, symmetric-diff count + +Run: + RUN_INTEGRATION_TESTS=1 pytest tests/test_grib_snap_out_comparison.py -v -s +""" + +import json +import os +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required") + +pytestmark = pytest.mark.skipif( + os.environ.get("RUN_INTEGRATION_TESTS", "0") != "1", + reason="Integration test — set RUN_INTEGRATION_TESTS=1 to run", +) + + +BASE_MARS_KEYS = { + "class": "od", + "stream": "oper", + "type": "fc", + "levtype": "sfc", + "param": "2t", + "step": "0", + "date": "-1", + "time": "0000", +} +COLLECTION = "ecmwf-mars" + +# Tolerance for pairing "the same" grid point across the two sources. +# At O1280, grid spacing is ~0.07° near the equator, so 0.02° is a safe +# nearest-neighbour bucket that won't cross-pair distinct grid points. +POINT_MATCH_TOL = 0.02 # degrees + + +# --------------------------------------------------------------------------- +# Bbox test cases — each stresses a different edge of the snap-out logic +# --------------------------------------------------------------------------- + +BBOX_CASES = [ + pytest.param( + [0.18, 0.00, 0.03, 0.15], # N W S E + id="tiny-equatorial-off-grid", # edges land between grid points + ), + pytest.param( + [0.50, 0.00, 0.00, 0.50], + id="small-equatorial-round-numbers", # nominally clean bounds + ), + pytest.param( + [45.10, 10.10, 44.90, 10.30], + id="mid-latitude-tiny", # denser grid; N pole reduction different + ), + pytest.param( + [1.00, -0.10, -1.00, 0.10], + id="crosses-prime-meridian-signed-lons", + # Meridian-straddling region expressed with signed longitudes + # (W < E numerically, W negative). This is the canonical polytope + # convention: bbox uses [[SW],[NE]] with signed lons — wrap syntax + # (W > E) is rejected/undefined at the polytope-fe layer. + ), +] + + +# --------------------------------------------------------------------------- +# Retrieval helpers +# --------------------------------------------------------------------------- + + +def _retrieve_mars_area(bbox, output_path): + from polytope.api import Client + + req = {**BASE_MARS_KEYS, "area": f"{bbox[0]}/{bbox[1]}/{bbox[2]}/{bbox[3]}"} + print(f"\n MARS area request: {json.dumps(req)}") + Client(quiet=True).retrieve(COLLECTION, req, output_file=output_path, pointer=False) + + +def _retrieve_boundingbox_covjson(bbox): + from polytope.api import Client + + req = { + **BASE_MARS_KEYS, + "feature": { + "type": "boundingbox", + "points": [[bbox[0], bbox[1]], [bbox[2], bbox[3]]], + }, + } + print(f" BBox request: {json.dumps(req)}") + with tempfile.NamedTemporaryFile(suffix=".covjson", delete=False, mode="w") as tmp: + p = tmp.name + try: + Client(quiet=True).retrieve(COLLECTION, req, output_file=p) + with open(p) as f: + return json.load(f) + finally: + if os.path.exists(p): + os.unlink(p) + + +# --------------------------------------------------------------------------- +# Point extraction +# --------------------------------------------------------------------------- + + +def _mars_grid_points(grib_path): + """Extract (lat, lon) of every point in a GRIB message, with fallbacks. + + Returns (pts, corners, npoints). `pts` may be None if eccodes can't + build a geometry for the sub-area (which is itself a diagnostic worth + reporting). + """ + with open(grib_path, "rb") as f: + gid = eccodes.codes_grib_new_from_file(f) + try: + corners = { + "N": eccodes.codes_get(gid, "latitudeOfFirstGridPointInDegrees"), + "W": eccodes.codes_get(gid, "longitudeOfFirstGridPointInDegrees"), + "S": eccodes.codes_get(gid, "latitudeOfLastGridPointInDegrees"), + "E": eccodes.codes_get(gid, "longitudeOfLastGridPointInDegrees"), + } + npoints = eccodes.codes_get(gid, "numberOfDataPoints") + + # Try iterator first + pts = None + try: + iterid = eccodes.codes_grib_iterator_new(gid, 0) + try: + pts = [] + while True: + r = eccodes.codes_grib_iterator_next(iterid) + if not r: + break + lat, lon, _ = r + pts.append((lat, lon)) + finally: + eccodes.codes_grib_iterator_delete(iterid) + except Exception as e: + print(f" (iterator failed: {e}; trying latitudes/longitudes arrays)") + pts = None + + # Fallback to raw lat/lon arrays + if pts is None: + try: + lats = eccodes.codes_get_array(gid, "latitudes") + lons = eccodes.codes_get_array(gid, "longitudes") + pts = list(zip(lats.tolist(), lons.tolist())) + except Exception as e: + print( + f" (latitudes/longitudes also failed: {e}; " + f"eccodes cannot decode geometry — reporting counts only)" + ) + pts = None + + return pts, corners, npoints + finally: + eccodes.codes_release(gid) + + +def _covjson_grid_points(covjson): + """Extract (lat, lon) tuples from the first coverage's composite axis.""" + cov = covjson["coverages"][0] + coords = cov["domain"]["axes"]["composite"]["values"] + return [(c[0], c[1]) for c in coords] + + +# --------------------------------------------------------------------------- +# Comparison +# --------------------------------------------------------------------------- + + +def _nearest_match(target, candidates, tol): + """Return (best_match, distance) or (None, inf) if no candidate within tol.""" + best = None + best_d = float("inf") + for c in candidates: + d = max(abs(target[0] - c[0]), abs(target[1] - c[1])) + if d < best_d: + best_d = d + best = c + return (best, best_d) if best_d <= tol else (None, best_d) + + +def _normalise_lon(lon): + """Map lon into (-180, 180] so [0,360) and [-180,180] representations compare equal.""" + x = ((lon + 180.0) % 360.0) - 180.0 + return x if x != -180.0 else 180.0 + + +def _normalise_pts(pts): + return [(lat, _normalise_lon(lon)) for lat, lon in pts] + + +def _compare_point_sets(mars_pts, mars_corners, mars_npoints, bbox_pts, bbox): + lines = [] + lines.append("") + lines.append(f" bbox (N/W/S/E): {bbox}") + lines.append( + f" MARS corners kept: N={mars_corners['N']:.6f} W={mars_corners['W']:.6f} " + f"S={mars_corners['S']:.6f} E={mars_corners['E']:.6f}" + ) + lines.append(f" MARS numberOfDataPoints: {mars_npoints}") + lines.append(f" BBox point count: {len(bbox_pts)}") + + if mars_pts is None: + lines.append(" MARS points: ") + return { + "mars_count": mars_npoints, + "bbox_count": len(bbox_pts), + "matched": 0, + "unmatched_mars": None, + "unmatched_bbox": None, + "max_delta": None, + "count_diff": abs(mars_npoints - len(bbox_pts)), + "report": "\n".join(lines), + } + + lines.append(f" MARS point count: {len(mars_pts)}") + + if mars_pts: + lats = [p[0] for p in mars_pts] + lons = [p[1] for p in mars_pts] + lines.append(f" MARS lat range: [{min(lats):.6f}, {max(lats):.6f}]") + lines.append(f" MARS lon range: [{min(lons):.6f}, {max(lons):.6f}]") + + if bbox_pts: + lats = [p[0] for p in bbox_pts] + lons = [p[1] for p in bbox_pts] + lines.append(f" BBox lat range: [{min(lats):.6f}, {max(lats):.6f}]") + lines.append(f" BBox lon range: [{min(lons):.6f}, {max(lons):.6f}]") + + # Pair up points by nearest-neighbour with a tight tolerance. + # Normalise longitudes to (-180, 180] so [0,360) vs [-180,180] + # representations don't spuriously mismatch across the prime meridian. + matched = 0 + max_delta = 0.0 + unmatched_mars = [] + remaining = _normalise_pts(bbox_pts) + for mp in _normalise_pts(mars_pts): + m, d = _nearest_match(mp, remaining, POINT_MATCH_TOL) + if m is None: + unmatched_mars.append(mp) + else: + matched += 1 + max_delta = max(max_delta, d) + remaining.remove(m) + + lines.append(f" Matched pairs: {matched}") + lines.append(f" Unmatched in MARS: {len(unmatched_mars)}") + lines.append(f" Unmatched in BBox: {len(remaining)}") + lines.append(f" Max per-pair Δ: {max_delta:.6f}° (tol {POINT_MATCH_TOL}°)") + + if unmatched_mars: + lines.append(f" MARS-only sample: {unmatched_mars[:3]}") + if remaining: + lines.append(f" BBox-only sample: {remaining[:3]}") + + return { + "mars_count": len(mars_pts), + "bbox_count": len(bbox_pts), + "matched": matched, + "unmatched_mars": len(unmatched_mars), + "unmatched_bbox": len(remaining), + "max_delta": max_delta, + "count_diff": abs(len(mars_pts) - len(bbox_pts)), + "report": "\n".join(lines), + } + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("bbox", BBOX_CASES) +def test_snap_out_matches(bbox): + """polytope BoundingBox and MARS `area` should agree on grid points kept.""" + print(f"\n{'=' * 70}\n Snap-out case: {bbox}\n{'=' * 70}") + mars_f = tempfile.NamedTemporaryFile(suffix="_mars.grib", delete=False) + mars_f.close() + try: + _retrieve_mars_area(bbox, mars_f.name) + mars_pts, mars_corners, mars_npoints = _mars_grid_points(mars_f.name) + + covjson = _retrieve_boundingbox_covjson(bbox) + bbox_pts = _covjson_grid_points(covjson) + + result = _compare_point_sets(mars_pts, mars_corners, mars_npoints, bbox_pts, bbox) + print(result["report"]) + + # If eccodes can't decode the MARS geometry we can only compare counts. + if result["unmatched_mars"] is None: + if result["count_diff"] != 0: + pytest.fail( + f"snap-out count divergence: MARS={result['mars_count']} " + f"vs BBox={result['bbox_count']} (Δ={result['count_diff']})" + ) + return + + if result["unmatched_mars"] or result["unmatched_bbox"]: + pytest.fail( + f"snap-out divergence: " + f"{result['unmatched_mars']} MARS-only, " + f"{result['unmatched_bbox']} BBox-only " + f"(tol {POINT_MATCH_TOL}°) — see printed report above" + ) + finally: + if os.path.exists(mars_f.name): + os.unlink(mars_f.name) + + +if __name__ == "__main__": + os.environ["RUN_INTEGRATION_TESTS"] = "1" + for bbox in [c.values[0] for c in BBOX_CASES]: + print(f"\n{'=' * 70}\n Snap-out case: {bbox}\n{'=' * 70}") + mars_path = f"/tmp/mars_snap_{bbox[0]}_{bbox[1]}_{bbox[2]}_{bbox[3]}.grib" + _retrieve_mars_area(bbox, mars_path) + mars_pts, mars_corners = _mars_grid_points(mars_path) + covjson = _retrieve_boundingbox_covjson(bbox) + bbox_pts = _covjson_grid_points(covjson) + res = _compare_point_sets(mars_pts, mars_corners, bbox_pts, bbox) + print(res["report"]) diff --git a/tests/test_grib_snap_out_polygon.py b/tests/test_grib_snap_out_polygon.py new file mode 100644 index 0000000..6466094 --- /dev/null +++ b/tests/test_grib_snap_out_polygon.py @@ -0,0 +1,190 @@ +"""Polygon-vs-MARS snap-out comparison at the prime meridian. + +The BoundingBox path fails to wrap across longitude 0° (see +`test_grib_snap_out_comparison.py::crosses-prime-meridian`). This test +checks whether the **polygon** feature handles the same geometry +correctly, using two equivalent shapes: + + 1. Signed longitudes — vertices in [-0.1, +0.1], unambiguous. + 2. Wrap longitudes — vertices at 359.9 and 0.1, matching the bbox + form; tests whether the polygon walker treats + the 359.9→0.1 edge as a short arc across 0° + or a long arc across 180°. + +Reference: MARS `area=1.0/359.9/-1.0/0.1` returns 84 points in a narrow +~0.14°-wide band straddling the meridian. + +Run: + RUN_INTEGRATION_TESTS=1 pytest tests/test_grib_snap_out_polygon.py -v -s +""" + +import json +import os +import tempfile + +import pytest + +eccodes = pytest.importorskip("eccodes", reason="eccodes required") + +pytestmark = pytest.mark.skipif( + os.environ.get("RUN_INTEGRATION_TESTS", "0") != "1", + reason="Integration test — set RUN_INTEGRATION_TESTS=1 to run", +) + + +BASE_MARS_KEYS = { + "class": "od", + "stream": "oper", + "type": "fc", + "levtype": "sfc", + "param": "2t", + "step": "0", + "date": "-1", + "time": "0000", +} +COLLECTION = "ecmwf-mars" + +# MARS `area=1.0/359.9/-1.0/0.1` returns this many points (empirical). +MARS_WRAP_POINT_COUNT = 84 +POINT_MATCH_TOL = 0.02 # degrees + + +# Polygon geometry equivalent to bbox N=1, W=-0.1, S=-1, E=0.1 — the +# meridian-straddling region expressed with signed longitudes (the +# canonical polytope convention; wrap longitudes are invalid input). +POLYGON_CASES = [ + pytest.param( + [[1.0, -0.1], [1.0, 0.1], [-1.0, 0.1], [-1.0, -0.1], [1.0, -0.1]], + id="signed-lons-straddling-zero", + ), +] + + +def _retrieve_polygon_covjson(shape): + from polytope.api import Client + + req = { + **BASE_MARS_KEYS, + "feature": {"type": "polygon", "shape": shape}, + } + print(f" Polygon request: {json.dumps(req)}") + with tempfile.NamedTemporaryFile(suffix=".covjson", delete=False, mode="w") as tmp: + p = tmp.name + try: + Client(quiet=True).retrieve(COLLECTION, req, output_file=p) + with open(p) as f: + return json.load(f) + finally: + if os.path.exists(p): + os.unlink(p) + + +def _retrieve_mars_area_points(bbox_area): + """Retrieve MARS area and return (lat, lon) point list + count.""" + from polytope.api import Client + + req = {**BASE_MARS_KEYS, "area": bbox_area} + print(f" MARS area request: {json.dumps(req)}") + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + p = tmp.name + try: + Client(quiet=True).retrieve(COLLECTION, req, output_file=p, pointer=False) + with open(p, "rb") as f: + gid = eccodes.codes_grib_new_from_file(f) + try: + npoints = eccodes.codes_get(gid, "numberOfDataPoints") + pts = None + try: + it = eccodes.codes_grib_iterator_new(gid, 0) + pts = [] + try: + while True: + r = eccodes.codes_grib_iterator_next(it) + if not r: + break + pts.append((r[0], r[1])) + finally: + eccodes.codes_grib_iterator_delete(it) + except Exception as e: + print(f" (MARS iterator failed: {e}; count-only)") + return pts, npoints + finally: + eccodes.codes_release(gid) + finally: + if os.path.exists(p): + os.unlink(p) + + +def _covjson_points(covjson): + cov = covjson["coverages"][0] + coords = cov["domain"]["axes"]["composite"]["values"] + return [(c[0], c[1]) for c in coords] + + +def _normalise_lon(lon): + """Map lon into (-180, 180] for comparison.""" + x = ((lon + 180.0) % 360.0) - 180.0 + return x if x != -180.0 else 180.0 + + +def _nearest_match(target, candidates, tol): + best_d = float("inf") + best = None + for c in candidates: + d = max(abs(target[0] - c[0]), abs(target[1] - c[1])) + if d < best_d: + best_d = d + best = c + return (best, best_d) if best_d <= tol else (None, best_d) + + +@pytest.mark.parametrize("shape", POLYGON_CASES) +def test_polygon_wrap_matches_mars_area(shape): + """Polygon equivalents of the failing bbox should match MARS 84-pt wrap.""" + print(f"\n{'=' * 70}\n Polygon shape: {shape}\n{'=' * 70}") + + mars_pts, mars_npoints = _retrieve_mars_area_points("1.0/359.9/-1.0/0.1") + covjson = _retrieve_polygon_covjson(shape) + poly_pts = _covjson_points(covjson) + + print(f" MARS numberOfDataPoints: {mars_npoints}") + print(f" Polygon coverage count: {len(poly_pts)}") + + if mars_pts is not None: + mars_norm = [(lat, _normalise_lon(lon)) for lat, lon in mars_pts] + poly_norm = [(lat, _normalise_lon(lon)) for lat, lon in poly_pts] + print( + f" MARS lon range (norm): " f"[{min(l for _, l in mars_norm):.6f}, {max(l for _, l in mars_norm):.6f}]" + ) + print( + f" Polygon lon range (norm): " f"[{min(l for _, l in poly_norm):.6f}, {max(l for _, l in poly_norm):.6f}]" + ) + + # Pair by nearest neighbour on normalised coords + matched = 0 + remaining = list(poly_norm) + unmatched_mars = [] + max_delta = 0.0 + for mp in mars_norm: + m, d = _nearest_match(mp, remaining, POINT_MATCH_TOL) + if m is None: + unmatched_mars.append(mp) + else: + matched += 1 + max_delta = max(max_delta, d) + remaining.remove(m) + + print(f" Matched pairs: {matched}") + print(f" Unmatched in MARS: {len(unmatched_mars)}") + print(f" Unmatched in Polygon: {len(remaining)}") + print(f" Max per-pair Δ: {max_delta:.6f}°") + if unmatched_mars: + print(f" MARS-only sample: {unmatched_mars[:3]}") + if remaining: + print(f" Polygon-only sample: {remaining[:3]}") + + assert len(unmatched_mars) == 0, f"{len(unmatched_mars)} MARS points not covered by polygon" + assert len(remaining) == 0, f"{len(remaining)} polygon points not covered by MARS wrap area" + else: + # Fallback: count-only assertion + assert len(poly_pts) == mars_npoints, f"Point count mismatch: polygon={len(poly_pts)} vs MARS={mars_npoints}" diff --git a/tests/test_pymars2grib_smoke.py b/tests/test_pymars2grib_smoke.py new file mode 100644 index 0000000..553bc34 --- /dev/null +++ b/tests/test_pymars2grib_smoke.py @@ -0,0 +1,299 @@ +"""Smoke tests for pymars2grib backend. + +Tests the hybrid approach: mars2grib generates the GRIB template (sections +1, 4, 5 from MARS keys), then eccodes overwrites section 3 (grid definition) +with sub-area geometry and fills in the data values. + +NOTE: `skipSection3=True` is unusable with the current eccodes build — it +triggers the `gridSpec` string setter which is not implemented +(`GridSpec::pack_string not available`). Instead we let mars2grib produce a +normal template and then overwrite section 3 via eccodes. + +Run from the Python 3.11 venv where pymars2grib is importable: + source mars2grib-bundle/env.sh + python tests/test_pymars2grib_smoke.py +""" + +import os +import tempfile + + +def test_pymars2grib_basic_encode(): + """Verify pymars2grib can encode and eccodes can overwrite section 3.""" + from pymars2grib import Mars2Grib + + # Plain Mars2Grib (no skipSection3 — that path hits an unimplemented + # gridSpec::pack_string in eccodes). mars2grib builds sections 1,4,5; + # eccodes then overwrites section 3 with sub-area geometry. + encoder = Mars2Grib() + + values = [280.5, 281.0, 279.3, 282.1, 280.0, 281.5, 279.8, 280.2, 281.3] + mars = { + "origin": "ecmf", + "class": "od", + "stream": "oper", + "type": "fc", + "expver": "0001", + "date": 20250702, + "time": 0, + "step": 0, + "levtype": "sfc", + "param": 167, + "packing": "ccsds", + "grid": "O1280", + } + + # mars2grib produces a small template; eccodes will overwrite section 3 + data = encoder.encode(values, mars) + print(f" mars2grib encoded → {len(data)} bytes (template)") + assert len(data) > 0, "Empty GRIB output" + + # Now use eccodes to fill in section 3 (grid geometry) for our sub-area + import eccodes + + gid = eccodes.codes_new_from_message(data) + + # Set sub-area grid geometry + eccodes.codes_set(gid, "gridType", "reduced_gg") + eccodes.codes_set_long(gid, "numberOfParallelsBetweenAPoleAndTheEquator", 1280) + eccodes.codes_set_long(gid, "Nj", 3) # 3 latitude rows in our sub-area + eccodes.codes_set_long_array(gid, "pl", [3, 3, 3]) # 3 points per row + eccodes.codes_set_double(gid, "latitudeOfFirstGridPointInDegrees", 0.176) + eccodes.codes_set_double(gid, "latitudeOfLastGridPointInDegrees", 0.035) + eccodes.codes_set_double(gid, "longitudeOfFirstGridPointInDegrees", 0.0) + eccodes.codes_set_double(gid, "longitudeOfLastGridPointInDegrees", 0.141) + + # Set the actual data values + eccodes.codes_set_long(gid, "bitsPerValue", 16) + eccodes.codes_set_values(gid, values) + + # Write final GRIB + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as f: + path = f.name + eccodes.codes_write(gid, open(path, "wb")) + eccodes.codes_release(gid) + + # Read back and verify + try: + with open(path, "rb") as f: + gid = eccodes.codes_grib_new_from_file(f) + param_id = eccodes.codes_get(gid, "paramId") + grid_type = eccodes.codes_get(gid, "gridType") + n_points = eccodes.codes_get(gid, "numberOfDataPoints") + decoded = eccodes.codes_get_values(gid).tolist() + eccodes.codes_release(gid) + + print(f" paramId: {param_id}") + print(f" gridType: {grid_type}") + print(f" numberOfDataPoints: {n_points}") + + assert param_id == 167, f"Expected paramId=167, got {param_id}" + assert grid_type == "reduced_gg", f"Expected reduced_gg, got {grid_type}" + assert n_points == len(values), f"Expected {len(values)} points, got {n_points}" + + max_diff = max(abs(a - b) for a, b in zip(values, decoded)) + print(f" max value diff: {max_diff:.2e}") + assert max_diff < 0.01, f"Value mismatch too large: {max_diff}" + + print("✓ Hybrid mars2grib+eccodes encode works") + finally: + os.unlink(path) + + +def test_pymars2grib_skip_section3_ensemble(): + """Verify ensemble fields resolve correctly (perturbationNumber).""" + from pymars2grib import Mars2Grib + + encoder = Mars2Grib() + + values = [280.5, 281.0, 279.3] + mars = { + "origin": "ecmf", + "class": "od", + "stream": "enfo", + "type": "pf", + "expver": "0001", + "date": 20250702, + "time": 0, + "step": 0, + "levtype": "sfc", + "param": 167, + "packing": "ccsds", + "grid": "O1280", + "number": 1, + } + + # Ensemble encoding requires `numberOfForecastsInEnsemble` in misc — + # mars2grib has no default deduction for it. Aborts the process otherwise. + misc = {"numberOfForecastsInEnsemble": 51} + + data = encoder.encode(values, mars, misc) + print(f" Ensemble template: {len(data)} bytes") + assert len(data) > 0 + + # Verify mars2grib set the ensemble keys correctly + import eccodes + + gid = eccodes.codes_new_from_message(data) + + try: + pdt = eccodes.codes_get(gid, "productDefinitionTemplateNumber") + print(f" productDefinitionTemplateNumber: {pdt}") + # Template 1 = individual ensemble forecast member + assert pdt == 1, f"Expected template 1 for ensemble, got {pdt}" + + perturb = eccodes.codes_get(gid, "perturbationNumber") + print(f" perturbationNumber: {perturb}") + assert perturb == 1, f"Expected perturbationNumber=1, got {perturb}" + + print("✓ Ensemble template correctly resolved by mars2grib") + finally: + eccodes.codes_release(gid) + + +def test_pymars2grib_mars_keys_only(): + """Verify mars2grib correctly resolves MARS keys to GRIB keys (sections 1,4,5).""" + from pymars2grib import Mars2Grib + + encoder = Mars2Grib() + + mars = { + "origin": "ecmf", + "class": "od", + "stream": "oper", + "type": "fc", + "expver": "0001", + "date": 20250702, + "time": 1200, + "step": 6, + "levtype": "pl", + "levelist": 500, + "param": 130, # temperature + "packing": "ccsds", + "grid": "O1280", + } + + data = encoder.encode([250.0] * 10, mars) + assert len(data) > 0 + + import eccodes + + gid = eccodes.codes_new_from_message(data) + try: + print(f" paramId: {eccodes.codes_get(gid, 'paramId')}") + print(f" dataDate: {eccodes.codes_get(gid, 'dataDate')}") + print(f" dataTime: {eccodes.codes_get(gid, 'dataTime')}") + print(f" stepRange: {eccodes.codes_get(gid, 'stepRange')}") + + assert eccodes.codes_get(gid, "paramId") == 130 + assert eccodes.codes_get(gid, "dataDate") == 20250702 + assert eccodes.codes_get(gid, "dataTime") in (12, 1200) + assert eccodes.codes_get(gid, "level") in (5, 500) + + print("✓ MARS key resolution correct (sections 1, 4, 5)") + finally: + eccodes.codes_release(gid) + + +def test_pymars2grib_via_backend_factory(): + """Verify the covjsonkit backend factory picks up pymars2grib.""" + from covjsonkit.decoder.grib_backends import get_backend + from covjsonkit.decoder.grib_backends.mars2grib_backend import Mars2GribBackend + + backend = get_backend("auto") + print(f"\n get_backend('auto') → {type(backend).__name__}") + assert isinstance(backend, Mars2GribBackend), f"Expected Mars2GribBackend, got {type(backend).__name__}" + print("✓ Auto-detection picks mars2grib when available") + + backend = get_backend("mars2grib") + assert isinstance(backend, Mars2GribBackend) + print("✓ Explicit 'mars2grib' works") + + +def test_pymars2grib_full_covjson_roundtrip(): + """Full test: load CovJSON fixture → to_grib(backend='mars2grib') → verify.""" + import json + + from covjsonkit.api import Covjsonkit + + fixture_path = os.path.join(os.path.dirname(__file__), "data", "test_oper_multipoint_coverage.json") + with open(fixture_path) as f: + covjson = json.load(f) + + decoder = Covjsonkit().decode(covjson) + + with tempfile.NamedTemporaryFile(suffix=".grib", delete=False) as tmp: + output = tmp.name + + try: + result = decoder.to_grib(output, backend="mars2grib") + assert result == output + file_size = os.path.getsize(output) + print(f"\n to_grib(backend='mars2grib') wrote {file_size} bytes") + + # Read back and verify + import eccodes + + messages = [] + with open(output, "rb") as f: + while True: + gid = eccodes.codes_grib_new_from_file(f) + if gid is None: + break + msg = { + "paramId": eccodes.codes_get(gid, "paramId"), + "shortName": eccodes.codes_get(gid, "shortName"), + "gridType": eccodes.codes_get(gid, "gridType"), + "numberOfDataPoints": eccodes.codes_get(gid, "numberOfDataPoints"), + "values": eccodes.codes_get_values(gid).tolist(), + } + eccodes.codes_release(gid) + messages.append(msg) + + print(f" Messages: {len(messages)}") + assert len(messages) == 2, f"Expected 2 messages, got {len(messages)}" + + params = sorted(m["shortName"] for m in messages) + print(f" Params: {params}") + assert "2t" in params + assert "10v" in params + + for msg in messages: + assert msg["numberOfDataPoints"] == 9 + assert msg["gridType"] == "reduced_gg" + + print("✓ Full CovJSON → mars2grib → GRIB roundtrip works") + finally: + os.unlink(output) + + +if __name__ == "__main__": + print("=" * 60) + print("pymars2grib Smoke Test") + print("=" * 60) + + tests = [ + test_pymars2grib_basic_encode, + test_pymars2grib_skip_section3_ensemble, + test_pymars2grib_mars_keys_only, + test_pymars2grib_via_backend_factory, + test_pymars2grib_full_covjson_roundtrip, + ] + + passed = 0 + failed = 0 + for test in tests: + print(f"\n--- {test.__name__} ---") + try: + test() + passed += 1 + except Exception as e: + print(f" ✗ FAILED: {e}") + import traceback + + traceback.print_exc() + failed += 1 + + print(f"\n{'='*60}") + print(f"Results: {passed} passed, {failed} failed") + print("=" * 60)