diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index 26fef03..819f03e 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -4,7 +4,7 @@ import pandas as pd -from .encoder import Encoder +from .encoder import Encoder, normalize_step_value class Position(Encoder): @@ -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 = {} @@ -256,6 +255,111 @@ def from_polytope(self, result, date_key: str = "date") -> dict: return self.covjson + def from_polytope_reforecast(self, result) -> dict: + """Encode reforecast/reanalysis data into a PointSeries (Position) collection. + + For legacy merged trees (no independent ``time`` axis) this delegates to + :meth:`from_polytope` with ``date_key="hdate"``. For separate-datetime + (``class=ce``) trees where ``date``, ``hdate`` and ``time`` are independent + axes, the reference datetime is ``hdate + time`` and one coverage is + produced per ``(point, levelist, number, reference)`` with a ``t`` axis + spanning ``reference + step`` across all steps. + """ + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + self.add_reference( + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ) + + exclude_meta = { + "latitude", + "longitude", + "hdate", + "time", + "step", + "param", + "number", + "levelist", + } + + coverages = {} + coverage_order = [] + param_order = [] + + for rec in self._reforecast_records(result): + value = float(rec["__value__"]) + lat = float(rec["latitude"]) + lon = float(rec["longitude"]) + level = rec.get("levelist", 0) + try: + level = int(level) + except (TypeError, ValueError): + pass + number = rec.get("number", 0) + try: + number = int(number) + except (TypeError, ValueError): + pass + para = rec.get("param") + step = normalize_step_value(rec.get("step", 0)) + ref = self._reforecast_reference(rec) + valid = ref + self._reforecast_step_timedelta(step) + valid_iso = valid.isoformat() + "Z" + + key = (lat, lon, level, number, ref.isoformat()) + if key not in coverages: + meta = {} + for name in rec: + if name == "__value__" or name in exclude_meta: + continue + meta[name] = self._reforecast_stringify(rec[name]) + meta["number"] = number + meta["Forecast date"] = ref.isoformat() + "Z" + coverages[key] = { + "lat": lat, + "lon": lon, + "level": level, + "meta": meta, + "values": {}, + "times": {}, + } + coverage_order.append(key) + + if para not in param_order: + param_order.append(para) + cov = coverages[key] + cov["times"][valid_iso] = None + cov["values"].setdefault(para, {})[valid_iso] = value + + if not coverages: + raise ValueError("No data was returned.") + + for para in param_order: + self.add_parameter(para) + + for key in coverage_order: + cov = coverages[key] + times = sorted(cov["times"].keys()) + coords = { + "latitude": [cov["lat"]], + "longitude": [cov["lon"]], + "levelist": [cov["level"]], + "t": times, + } + val_dict = {} + for para, time_vals in cov["values"].items(): + val_dict[para] = [time_vals[t] for t in times] + self.add_coverage(cov["meta"], coords, val_dict) + + return self.covjson + def from_polytope_month(self, result): coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/TimeSeries.py b/covjsonkit/encoder/TimeSeries.py index 9f3f5eb..d6a6c89 100644 --- a/covjsonkit/encoder/TimeSeries.py +++ b/covjsonkit/encoder/TimeSeries.py @@ -54,15 +54,24 @@ def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata @staticmethod - def _hdate_step_timestamp(date, step): + def _hdate_step_timestamp(date, step, time_offset=None): """Return the valid-time (as a datetime) for a given hdate and step. Mirrors the per-hdate stamp computation used in ``from_polytope`` so the collapsed reanalysis path produces identical timestamps. + + ``time_offset`` is an optional time-of-day offset (a ``timedelta``) coming + from an independent ``time`` axis (the separate-datetime reforecast + representation). When ``None`` the timestamp reduces to ``hdate + step``, + preserving the legacy merged-tree behaviour. """ date_format = "%Y%m%dT%H%M%S" + if isinstance(date, str) and date.endswith("Z"): + date = date[:-1] new_date = pd.Timestamp(date).strftime(date_format) start_time = datetime.strptime(new_date, date_format) + if time_offset is not None: + start_time = start_time + time_offset if isinstance(step, timedelta): return start_time + step try: @@ -87,11 +96,14 @@ def _collapse_reanalysis(self, fields, coords, mars_metadata, range_dict): points = len(coords[fields["dates"][0]]["composite"]) first_date = fields["dates"][0] + # Optional independent time-of-day offset (separate-datetime reforecast). + time_offset = fields.get("time_offset") + # Ordered list of (stamp, date, step) across every hdate/step combination. stamp_order = [] for date in fields["dates"]: for step in fields["step"]: - stamp = self._hdate_step_timestamp(date, step) + stamp = self._hdate_step_timestamp(date, step, time_offset) stamp_order.append((stamp, date, step)) # Stable sort by valid-time so ties preserve insertion order. stamp_order.sort(key=lambda x: x[0]) @@ -175,7 +187,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 = {} @@ -200,13 +211,16 @@ 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", reforecast: bool = False) -> dict: """Encode a polytope ``TensorIndexTree`` result into a PointSeries CoverageJSON collection. Args: result: The polytope ``TensorIndexTree`` containing the data to be converted. date_key: Tree axis name to treat as the time dimension (``"date"`` for forecasts, ``"hdate"`` for hindcast/reforecast). + reforecast: When ``True`` the separate-datetime reforecast walker is used, + which keeps an independent ``time`` axis as a scalar time-of-day offset + instead of folding it into the ``hdate`` time dimension. Returns: dict: The CoverageJSON representation of the coverage collection. """ @@ -223,7 +237,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 - self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) + if reforecast: + self.walk_tree_reforecast(result, fields, coords, mars_metadata, range_dict) + else: + self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) end = time.time() delta = end - start logging.debug("Tree walking ends at: %s", end) # noqa: E501 @@ -341,6 +358,208 @@ def from_polytope(self, result, date_key: str = "date") -> dict: return self.covjson + def from_polytope_reforecast(self, result) -> dict: + """Encode separate-datetime reforecast/reanalysis data (``class=ce``). + + Handles a result tree in which ``date``, ``hdate`` and ``time`` are + **independent axes** (rather than a single merged datetime). ``hdate`` is + the reforecast time dimension and ``time`` an independent time-of-day + axis; both may be multi-valued and compressed into a single tree node. + + Every ``(hdate, time, step)`` combination for a spatial point yields a + valid-time ``hdate + time + step``; all such valid-times for a point are + collapsed into a single PointSeries coverage whose ``t``-axis holds them + chronologically sorted, with the parameter values in matching order. + + Backward compatible with the legacy tree where ``time`` was pre-merged + into ``hdate`` (no separate ``time`` node) and ``hdate`` was branched: + the timestamps then reduce to ``hdate + step``. + """ + import itertools + + import numpy as np + + from .encoder import is_merged_node + + # Spatial reference system (temporal RS is implied by the t-axis). + self.add_reference( + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ) + + # Axes that must not leak into the per-coverage mars:metadata block. + exclude_meta = { + "latitude", + "longitude", + "hdate", + "time", + "step", + "param", + "number", + "levelist", + } + + # coverage_key -> accumulated coverage state. + coverages = {} + coverage_order = [] + param_order = [] + # Set True as soon as a forecast (efas) coverage is created, so the final + # emission can be reordered date -> time -> point (see below). + forecast_result = False + + def stringify(value): + if isinstance(value, (np.datetime64, np.timedelta64)): + # np.datetime64 already stringifies as ISO-8601 (e.g. 2024-03-01 + # or 2024-03-01T00:00:00); np.timedelta64 has no ISO date form. + return str(value) + if isinstance(value, (pd.Timestamp, datetime)): + # pandas/py datetimes render as "YYYY-MM-DD HH:MM:SS" via str(); + # emit ISO-8601 ("YYYY-MM-DDTHH:MM:SS") so metadata dates are + # spec-compliant. + return value.isoformat() + if isinstance(value, timedelta): + return str(value) + return value + + def emit(full_path, flat_result): + nonlocal forecast_result + axis_names = [name for name, _ in full_path] + axis_values = [values for _, values in full_path] + for idx, combo in enumerate(itertools.product(*axis_values)): + value = flat_result[idx] + if value is None: + continue + d = dict(zip(axis_names, combo)) + lat = float(d["latitude"]) + lon = float(d["longitude"]) + level = d.get("levelist", 0) + try: + level = int(level) + except (TypeError, ValueError): + pass + number = d.get("number", 0) + try: + number = int(number) + except (TypeError, ValueError): + pass + para = d.get("param") + step = d.get("step", 0) + hdate = d.get("hdate", d.get("date")) + time_off = d.get("time") + if isinstance(time_off, np.timedelta64): + time_off = pd.to_timedelta(time_off).to_pytimedelta() + stamp = self._hdate_step_timestamp(hdate, step, time_off) + + # Stream determines the coverage grouping for class=ce: + # * efcl (reforecast/reanalysis): collapse all hdates into a + # single coverage whose t-axis holds the hdates. + # * efas (forecast): date+time is the forecast *run* reference; + # each (date, time) run is its own coverage and the steps + # become that coverage's t-axis. + # * anything else (e.g. enfh): legacy one-coverage-per-hdate + # with a scalar "Forecast date". + stream = d.get("stream") + is_ce = d.get("class") == "ce" + collapse = is_ce and stream == "efcl" + forecast = is_ce and stream == "efas" + if collapse: + key = (lat, lon, level, number) + elif forecast: + # Reference datetime of the forecast run = date + time. + reference = self._hdate_step_timestamp(hdate, 0, time_off) + key = (lat, lon, level, number, reference) + else: + key = (lat, lon, level, number, stringify(hdate)) + if key not in coverages: + meta = {} + for name in axis_names: + if name in exclude_meta: + continue + meta[name] = stringify(d[name]) + meta["number"] = number + meta["levelist"] = level + if forecast: + # date+time is folded into the run reference; drop the + # raw date axis so it doesn't duplicate "Forecast date". + meta.pop("date", None) + meta["Forecast date"] = reference.isoformat() + "Z" + elif not collapse: + meta["Forecast date"] = pd.Timestamp(hdate).isoformat() + "Z" + coverages[key] = { + "lat": lat, + "lon": lon, + "level": level, + "number": number, + "meta": meta, + "params": {}, + } + if forecast: + forecast_result = True + # Ordering key so coverages come out date -> time -> point + # (reference = date + time). Insertion/tree order is + # otherwise time-major for multi-date forecast requests. + coverages[key]["sort_key"] = (reference, lat, lon, level, number) + coverage_order.append(key) + + if para not in param_order: + param_order.append(para) + params = coverages[key]["params"] + params.setdefault(para, []).append((stamp, float(value))) + + def recurse(node, path): + children = node.children + if len(children) == 0: + # Leaf longitude node: ``path`` already carries latitude/longitude. + emit(path, node.result) + return + for child in children: + if is_merged_node(child): + # Compacted lat/lon leaf: carries its own (lat, lon) + result. + lat, lon = child.values[0], child.values[1] + merged_path = path + [ + ("latitude", (lat,)), + ("longitude", (lon,)), + ] + emit(merged_path, child.result) + continue + recurse(child, path + [(child.axis.name, tuple(child.values))]) + + recurse(result, []) + + if not coverages: + raise ValueError("No data was returned.") + + for para in param_order: + self.add_parameter(para) + + # Forecast (efas) coverages: emit in date -> time -> point order. + if forecast_result: + coverage_order.sort(key=lambda k: coverages[k]["sort_key"]) + + for key in coverage_order: + cov = coverages[key] + params = cov["params"] + paras = list(params.keys()) + reference = params[paras[0]] + # Stable sort by valid-time; ties preserve insertion order. + order = sorted(range(len(reference)), key=lambda i: reference[i][0]) + t_values = [reference[i][0].isoformat() + "Z" for i in order] + val_dict = {para: [params[para][i][1] for i in order] for para in paras} + coord_entry = { + "latitude": [cov["lat"]], + "longitude": [cov["lon"]], + "levelist": [cov["level"]], + "t": t_values, + } + self.add_coverage(cov["meta"], coord_entry, val_dict) + + return self.covjson + def from_polytope_month(self, result): """Convert a Polytope result for monthly-mean streams (e.g. clmn) into CovJSON. diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 0781f89..70cc0e5 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: @@ -240,6 +239,114 @@ def from_polytope(self, result, date_key: str = "date") -> dict: return self.covjson + def from_polytope_reforecast(self, result) -> dict: + """Encode reforecast/reanalysis data into a VerticalProfile collection. + + For legacy merged trees (no independent ``time`` axis) this delegates to + :meth:`from_polytope` with ``date_key="hdate"``. For separate-datetime + (``class=ce``) trees where ``date``, ``hdate`` and ``time`` are independent + axes, the reference datetime is ``hdate + time`` and one coverage is + produced per ``(point, number, reference, step)`` holding all levels in + its ``levelist`` axis and a single valid time ``reference + step``. + """ + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + self.add_reference( + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ) + + exclude_meta = { + "latitude", + "longitude", + "hdate", + "time", + "step", + "param", + "number", + "levelist", + } + + coverages = {} + coverage_order = [] + param_order = [] + + for rec in self._reforecast_records(result): + value = float(rec["__value__"]) + lat = float(rec["latitude"]) + lon = float(rec["longitude"]) + level = rec.get("levelist", 0) + try: + level = int(level) + except (TypeError, ValueError): + pass + number = rec.get("number", 0) + try: + number = int(number) + except (TypeError, ValueError): + pass + para = rec.get("param") + step = normalize_step_value(rec.get("step", 0)) + ref = self._reforecast_reference(rec) + valid = ref + self._reforecast_step_timedelta(step) + valid_iso = valid.isoformat() + "Z" + + key = (lat, lon, number, ref.isoformat(), str(step)) + if key not in coverages: + meta = {} + for name in rec: + if name == "__value__" or name in exclude_meta: + continue + meta[name] = self._reforecast_stringify(rec[name]) + meta["number"] = number + meta["step"] = step + meta["Forecast date"] = ref.isoformat() + "Z" + coverages[key] = { + "lat": lat, + "lon": lon, + "t": valid_iso, + "meta": meta, + # per param: {level: value} + "values": {}, + # level -> None (ordered set) + "levels": {}, + } + coverage_order.append(key) + + if para not in param_order: + param_order.append(para) + cov = coverages[key] + cov["levels"][level] = None + cov["values"].setdefault(para, {})[level] = value + + if not coverages: + raise ValueError("No data was returned.") + + for para in param_order: + self.add_parameter(para) + + for key in coverage_order: + cov = coverages[key] + levels = sorted(cov["levels"].keys()) + coords = { + "latitude": [cov["lat"]], + "longitude": [cov["lon"]], + "levelist": list(levels), + "t": [cov["t"]], + } + val_dict = {} + for para, level_vals in cov["values"].items(): + val_dict[para] = [level_vals[lev] for lev in levels] + self.add_coverage(cov["meta"], coords, val_dict) + + return self.covjson + def from_polytope_month(self, result): coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/encoder.py b/covjsonkit/encoder/encoder.py index 9e97760..6f74b3f 100644 --- a/covjsonkit/encoder/encoder.py +++ b/covjsonkit/encoder/encoder.py @@ -1,7 +1,7 @@ from __future__ import annotations from abc import ABC, abstractmethod -from datetime import timedelta +from datetime import datetime, timedelta from typing import Any import numpy as np @@ -471,6 +471,138 @@ def emit_leaf(lat, lon_values, result): else: emit_leaf(fields["lat"], tree.values, tree.result) + def walk_tree_reforecast(self, tree, fields, coords, mars_metadata, range_dict): + """Walk the result tree for reforecast/reanalysis with an independent ``time`` axis. + + Unlike :meth:`walk_tree` (which, for the legacy merged representation, + folds any ``time`` axis into the ``hdate`` time dimension), this walker + treats ``hdate`` as the branching time axis and captures the separate + ``time`` axis as a scalar time-of-day offset stored in + ``fields["time_offset"]``. The valid-time for each hdate is then + ``hdate + time_offset + step`` (computed downstream in the collapse step). + + For efcl the ``date`` and ``step`` axes are single-valued and ``time`` is a + single timedelta, so the offset is a scalar. Backward compatibility with + merged trees (no separate ``time`` node) is preserved: ``time_offset`` + simply stays absent and the timestamps reduce to ``hdate + step``. + """ + date_key = "hdate" + + def create_composite_key(date, level, num, para, s): + return (date, level, num, para, s) + + def handle_non_leaf_node(child): + non_leaf_axes = ["latitude", "longitude", "param", date_key, "time"] + if child.axis.name not in non_leaf_axes: + val = child.values[0] + if isinstance(val, np.datetime64): + val = str(val) + elif isinstance(val, timedelta): + val = timedelta_to_step_string(val) + elif child.axis.name == "step": + # Step is not a timedelta! Need to normalize it + val = normalize_step_value(val) + mars_metadata[child.axis.name] = val + + def handle_specific_axes(child): + if child.axis.name == "latitude": + return child.values[0] + if child.axis.name == "levelist": + return child.values + if child.axis.name == "param": + return child.values + if child.axis.name == date_key: + dates = [f"{date}Z" for date in child.values] + mars_metadata["Forecast date"] = str(child.values[0]) + for date in dates: + coords[date] = {} + coords[date]["composite"] = [] + coords[date]["t"] = [date] + return dates + if child.axis.name == "time": + # Independent time-of-day axis: capture as a scalar offset rather + # than folding it into the hdate time dimension. efcl guarantees a + # single time value; take the first if a span is ever returned. + fields["time_offset"] = child.values[0] + return None + if child.axis.name == "number": + return child.values + if child.axis.name == "step": + return child.values + return None + + def calculate_index_bounds(level_len, num_len, para_len, step_len, l, i, j, k): # noqa: E741 + start_index = int(l * level_len) + int(i * num_len) + int(j * para_len) + int(k * step_len) + end_index = start_index + int(step_len) + return start_index, end_index + + def append_composite_coords(dates, tree_values, lat, coords): + for value in tree_values: + coords[dates]["composite"].append([lat, value]) + + def emit_leaf(lat, lon_values, result): + lon_values = [float(val) for val in lon_values] + if all(val is None for val in result): + fields["dates"] = fields["dates"][:-1] + for date in fields["dates"]: + for level in fields["levels"]: + for num in fields["number"]: + for para in fields["param"]: + for s in fields["step"]: + key = create_composite_key(date, level, num, para, s) + if key in range_dict: + del range_dict[key] + else: + result = [float(val) if val is not None else val for val in result] + level_len = len(result) / len(fields["levels"]) + num_len = level_len / len(fields["number"]) + para_len = num_len / len(fields["param"]) + step_len = para_len / len(fields["step"]) + + append_composite_coords(fields["dates"][-1], lon_values, lat, coords) + + for l, level in enumerate(fields["levels"]): # noqa: E741 + for i, num in enumerate(fields["number"]): + for j, para in enumerate(fields["param"]): + for k, s in enumerate(fields["step"]): + start_index, end_index = calculate_index_bounds( + level_len, num_len, para_len, step_len, l, i, j, k + ) + key = create_composite_key(fields["dates"][-1], level, num, para, s) + if key not in range_dict: + range_dict[key] = [] + range_dict[key].extend(result[start_index:end_index]) + + if len(tree.children) != 0: + for child in tree.children: + # Compacted unstructured leaf: values=(lat, lon), own result. Emit directly. + if is_merged_node(child): + emit_leaf(child.values[0], [child.values[1]], child.result) + continue + handle_non_leaf_node(child) + result = handle_specific_axes(child) + if result is not None: + if child.axis.name == "latitude": + fields["lat"] = result + elif child.axis.name == "levelist": + fields["levels"] = result + if "l" in fields: + fields["l"].extend(result) + elif child.axis.name == "param": + fields["param"] = result + elif child.axis.name == date_key: + fields["dates"].extend(result) + elif child.axis.name == "number": + fields["number"] = result + elif child.axis.name == "step": + fields["step"] = result + if "s" in fields: + fields["s"].extend(result) + + self.walk_tree_reforecast(child, fields, coords, mars_metadata, range_dict) + else: + emit_leaf(fields["lat"], tree.values, tree.result) + def walk_tree_step(self, tree, fields, coords, mars_metadata, range_dict): def create_composite_key_step(date, level, num, para): return (date, level, num, para) @@ -792,11 +924,224 @@ def from_xarray(self, dataset): def from_polytope(self, result, date_key: str = "date") -> dict: pass + @staticmethod + def _reforecast_stringify(value): + """Coerce datetime/timedelta-like values to strings for mars:metadata.""" + if isinstance(value, (pd.Timestamp, datetime)): + # str() renders "YYYY-MM-DD HH:MM:SS"; emit ISO-8601 instead. + return value.isoformat() + if isinstance(value, (np.datetime64, np.timedelta64, timedelta)): + return str(value) + return value + + @staticmethod + def _reforecast_timedelta(value) -> timedelta: + """Coerce a ``time``/offset value to a :class:`datetime.timedelta`.""" + if value is None: + return timedelta(0) + if isinstance(value, timedelta): + return value + return pd.to_timedelta(value).to_pytimedelta() + + @staticmethod + def _reforecast_step_timedelta(step) -> timedelta: + """Coerce a normalized ``step`` value to a :class:`datetime.timedelta`.""" + if isinstance(step, timedelta): + return step + if isinstance(step, np.timedelta64): + return pd.to_timedelta(step).to_pytimedelta() + if isinstance(step, str): + return timedelta(hours=parse_step_string(step)) + try: + return timedelta(hours=float(step)) + except (TypeError, ValueError): + return timedelta(0) + + @staticmethod + def _reforecast_reference(rec): + """Reference datetime for a record: ``hdate (or date) + time``, ISO ``...Z``.""" + hdate = rec.get("hdate", rec.get("date")) + ref = pd.Timestamp(hdate) + Encoder._reforecast_timedelta(rec.get("time")) + return ref + + @staticmethod + def _tree_has_axis(tree, axis_name) -> bool: + stack = [tree] + while stack: + node = stack.pop() + for child in getattr(node, "children", []): + if is_merged_node(child): + continue + if getattr(getattr(child, "axis", None), "name", None) == axis_name: + return True + stack.append(child) + return False + + def _reforecast_records(self, result): + """Flatten a reforecast result tree into per-value records. + + Every leaf value is expanded into a dict of ``{axis_name: value}`` for + the full root-to-leaf path (latitude/longitude included), using the same + ``itertools.product`` layout the compressed leaf ``result`` array follows. + Handles both the compacted ``MergedTensorIndexNode`` lat/lon leaves and + the classic nested latitude/longitude branches. + """ + import itertools + + records = [] + + def emit(full_path, flat_result): + axis_names = [name for name, _ in full_path] + axis_values = [values for _, values in full_path] + for idx, combo in enumerate(itertools.product(*axis_values)): + value = flat_result[idx] + if value is None: + continue + records.append(dict(zip(axis_names, combo))) + records[-1]["__value__"] = value + + def recurse(node, path): + children = node.children + if len(children) == 0: + emit(path, node.result) + return + for child in children: + if is_merged_node(child): + lat, lon = child.values[0], child.values[1] + emit( + path + [("latitude", (lat,)), ("longitude", (lon,))], + child.result, + ) + continue + recurse(child, path + [(child.axis.name, tuple(child.values))]) + + recurse(result, []) + return records + def from_polytope_reforecast(self, result) -> dict: """Encode reforecast/reanalysis data that uses ``"hdate"`` as the time axis. - Delegates to :meth:`from_polytope` with ``date_key="hdate"``. - Each hdate produces a separate coverage; steps within a single - hdate become that coverage's t-axis values. + Two representations are supported: + + * **Legacy merged** trees (no independent ``time`` node): ``hdate`` is the + branching time axis and each hdate/step produces its own coverage. This + is delegated to :meth:`from_polytope` with ``date_key="hdate"``. + * **Separate-datetime** trees (``class=ce``) where ``date``, ``hdate`` and + ``time`` are independent axes: the reforecast reference datetime is + ``hdate + time`` and one coverage is produced per + ``(reference-datetime, step, number)`` combination, holding all spatial + points in its ``composite`` axis. """ - return self.from_polytope(result, date_key="hdate") + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + self.add_reference( + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ) + + # Axes that should not leak into the per-coverage mars:metadata block. + exclude_meta = { + "latitude", + "longitude", + "hdate", + "time", + "step", + "param", + "number", + "levelist", + } + + def stringify(value): + return self._reforecast_stringify(value) + + def to_timedelta(value): + return self._reforecast_timedelta(value) + + coverages = {} + coverage_order = [] + param_order = [] + + for rec in self._reforecast_records(result): + value = rec["__value__"] + lat = float(rec["latitude"]) + lon = float(rec["longitude"]) + level = rec.get("levelist", 0) + try: + level = int(level) + except (TypeError, ValueError): + pass + number = rec.get("number", 0) + try: + number = int(number) + except (TypeError, ValueError): + pass + para = rec.get("param") + step = normalize_step_value(rec.get("step", 0)) + hdate = rec.get("hdate", rec.get("date")) + ref = pd.Timestamp(hdate) + to_timedelta(rec.get("time")) + valid = ref + self._reforecast_step_timedelta(step) + + # One coverage per (reference, step, number); reference is + # hdate + time for efcl and date + time (the run) for efas. + key = (ref, step, number) + if key not in coverages: + meta = {} + for name in rec: + if name == "__value__" or name in exclude_meta: + continue + meta[name] = stringify(rec[name]) + meta["number"] = number + meta["step"] = step + is_ce = rec.get("class") == "ce" + if is_ce and rec.get("stream") == "efas": + # date + time is the run reference; drop the raw date so it + # doesn't duplicate "Forecast date". + meta.pop("date", None) + meta["Forecast date"] = ref.isoformat() + "Z" + elif not (is_ce and rec.get("stream") == "efcl"): + # efcl coverages are delineated by their valid time, which + # the t-axis already carries, so they get no "Forecast date". + meta["Forecast date"] = ref.isoformat() + "Z" + coverages[key] = { + "valid": valid.isoformat() + "Z", + "points": [], + "point_index": {}, + "values": {}, + "meta": meta, + } + coverage_order.append(key) + + cov = coverages[key] + point = (lat, lon, level) + if point not in cov["point_index"]: + cov["point_index"][point] = len(cov["points"]) + cov["points"].append(point) + if para not in param_order: + param_order.append(para) + cov["values"].setdefault(para, {})[point] = float(value) + + if not coverages: + raise ValueError("No data was returned.") + + for para in param_order: + self.add_parameter(para) + + # Emit date -> time -> step -> number regardless of tree axis order. + coverage_order.sort(key=lambda k: (k[0], self._reforecast_step_timedelta(k[1]), k[2])) + + for key in coverage_order: + cov = coverages[key] + composite = [[lat, lon, level] for (lat, lon, level) in cov["points"]] + coords = {"composite": composite, "t": [cov["valid"]]} + val_dict = {} + for para, point_vals in cov["values"].items(): + val_dict[para] = [point_vals[pt] for pt in cov["points"]] + self.add_coverage(cov["meta"], coords, val_dict) + + return self.covjson diff --git a/covjsonkit/version.py b/covjsonkit/version.py index 8e523c9..3b9e925 100644 --- a/covjsonkit/version.py +++ b/covjsonkit/version.py @@ -1 +1 @@ -__version__ = "0.2.25" +__version__ = "0.2.26" diff --git a/pyproject.toml b/pyproject.toml index 533274e..99c0bd2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,7 +22,7 @@ markers = ["data: uses test data (deselect with '-m \"not data\"')",] [project] name = "covjsonkit" -version = "0.2.25" +version = "0.2.26" dependencies = [ "pandas<3", "orjson", diff --git a/tests/conftest.py b/tests/conftest.py index cb3041f..763fcf1 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -79,7 +79,7 @@ def make_merged_point(lat, lon, result): """ if MergedTensorIndexNode is None: raise RuntimeError( - "MergedTensorIndexNode is not available in this polytope build; " "merged-node tests should be skipped." + "MergedTensorIndexNode is not available in this polytope build; merged-node tests should be skipped." ) lat_ax = IntDatacubeAxis() lat_ax.name = "latitude" @@ -177,3 +177,141 @@ def reforecast_tree(branches, date=np.datetime64("2024-03-01")): for b in branches: root.add_child(b) return tree + + +def reforecast_separate_datetime_tree( + points, + hdates, + times, + date=np.datetime64("2024-03-01"), + step=(0,), + param="167", + point_factory=make_point, +): + """Build a separate-datetime reforecast tree (class=ce). + + ``date``, ``hdate`` and ``time`` are all present as independent axes, mirroring + the EFAS layout where the reforecast reference datetime is ``hdate + time``. + + Args: + points: list of (lat, lon, result_list) tuples. Each ``result_list`` must + be laid out as the flattened product over (hdate, step, time) in that + order (hdate-major), matching the compressed leaf ordering. + hdates: tuple of ``np.datetime64`` hdate values. + times: tuple of ``np.timedelta64``/``timedelta`` time-of-day offsets. + step: tuple of step values. + """ + branch = chain( + node("hdate", hdates), + node("domain", ("g",)), + node("expver", ("4321",)), + node("levtype", ("sfc",)), + node("param", (param,)), + node("step", step), + node("stream", ("efcl",)), + node("time", times), + node("type", ("sfo",)), + ) + parent = tip(branch) + for lat, lon, result in points: + parent.add_child(point_factory(lat, lon, result)) + + tree = chain( + TensorIndexTree(), + node("class", ("ce",)), + node("date", (date,)), + ) + tip(tree).add_child(branch) + return tree + + +def forecast_separate_datetime_tree( + points, + times, + date=np.datetime64("2026-05-01"), + step=(6, 24), + param="240023", + point_factory=make_point, +): + """Build a separate-datetime *forecast* tree (class=ce, stream=efas). + + Unlike :func:`reforecast_separate_datetime_tree` there is **no** ``hdate`` + axis: ``date`` + ``time`` is the forecast *run* reference and ``step`` is the + lead time. Each ``(date, time)`` run should become its own coverage with the + steps forming its ``t``-axis. + + Args: + points: list of (lat, lon, result_list) tuples. Each ``result_list`` must + be laid out as the flattened product over (date, step, time) in that + order, matching the compressed leaf ordering (date is single-valued + here, so effectively (step, time)). + times: tuple of ``np.timedelta64``/``timedelta`` time-of-day offsets. + step: tuple of step values. + """ + branch = chain( + node("date", (date,)), + node("domain", ("g",)), + node("expver", ("8888",)), + node("levtype", ("sfc",)), + node("model", ("lisflood",)), + node("origin", ("ecmf",)), + node("param", (param,)), + node("step", step), + node("stream", ("efas",)), + node("time", times), + node("type", ("cf",)), + ) + parent = tip(branch) + for lat, lon, result in points: + parent.add_child(point_factory(lat, lon, result)) + + tree = chain( + TensorIndexTree(), + node("class", ("ce",)), + ) + tip(tree).add_child(branch) + return tree + + +def reforecast_separate_datetime_vertical_tree( + points, + hdates, + times, + levels, + date=np.datetime64("2024-03-01"), + step=(0,), + param="130", + point_factory=make_point, +): + """Build a separate-datetime reforecast tree (class=ce) with a ``levelist`` axis. + + Mirrors :func:`reforecast_separate_datetime_tree` but inserts a multi-valued + ``levelist`` axis, for exercising the VerticalProfile reforecast path. + + Each ``result_list`` must be laid out as the flattened product over + (hdate, levelist, step, time) in that order, matching the compressed leaf + ordering (top-down axis order in the tree). + """ + branch = chain( + node("hdate", hdates), + node("domain", ("g",)), + node("expver", ("4321",)), + node("levtype", ("pl",)), + node("param", (param,)), + node("levelist", levels), + node("step", step), + node("stream", ("efcl",)), + node("time", times), + node("type", ("sfo",)), + ) + parent = tip(branch) + for lat, lon, result in points: + parent.add_child(point_factory(lat, lon, result)) + + tree = chain( + TensorIndexTree(), + node("class", ("ce",)), + node("date", (date,)), + ) + tip(tree).add_child(branch) + return tree diff --git a/tests/test_encoder_forecast_separate_datetime.py b/tests/test_encoder_forecast_separate_datetime.py new file mode 100644 index 0000000..9e14cce --- /dev/null +++ b/tests/test_encoder_forecast_separate_datetime.py @@ -0,0 +1,110 @@ +"""Tests for the separate-datetime *forecast* (class=ce, stream=efas) path in +TimeSeries.from_polytope_reforecast. + +Forecast data has no ``hdate`` axis: ``date + time`` is the forecast run +reference and ``step`` is the lead time. Each ``(date, time)`` run must become +its own coverage, with the steps forming that coverage's ``t``-axis -- unlike +the efcl reforecast path which collapses everything into a single coverage. +""" + +import numpy as np +from conftest import chain, forecast_separate_datetime_tree, make_point, node, tip +from polytope_feature.datacube.tensor_index_tree import TensorIndexTree + +from covjsonkit.api import Covjsonkit + +TIMES = (np.timedelta64(0, "h"), np.timedelta64(12, "h")) + + +def test_efas_forecast_splits_runs_by_datetime(): + # result layout: product(date, step, time) with single date -> (step, time): + # idx 0: step6 time0 -> 1 + # idx 1: step6 time1 -> 2 + # idx 2: step24 time0 -> 3 + # idx 3: step24 time1 -> 4 + points = [(50.73746373267438, 7.107723168102066, [1, 2, 3, 4])] + tree = forecast_separate_datetime_tree(points, TIMES, step=(6, 24)) + covjson = Covjsonkit().encode("CoverageCollection", "timeseries").from_polytope_reforecast(tree) + + # Two runs (2026-05-01 00:00 and 12:00), NOT a single collapsed coverage. + assert len(covjson["coverages"]) == 2 + + by_ref = {c["mars:metadata"]["Forecast date"]: c for c in covjson["coverages"]} + assert set(by_ref) == {"2026-05-01T00:00:00Z", "2026-05-01T12:00:00Z"} + + # Run 00:00 -> steps 6, 24 -> valid times 06:00 and next-day 00:00. + run0 = by_ref["2026-05-01T00:00:00Z"] + assert run0["domain"]["axes"]["t"]["values"] == [ + "2026-05-01T06:00:00Z", + "2026-05-02T00:00:00Z", + ] + (param_name,) = run0["ranges"].keys() + assert run0["ranges"][param_name]["values"] == [1, 3] + + # Run 12:00 -> steps 6, 24 -> valid times 18:00 and next-day 12:00. + run1 = by_ref["2026-05-01T12:00:00Z"] + assert run1["domain"]["axes"]["t"]["values"] == [ + "2026-05-01T18:00:00Z", + "2026-05-02T12:00:00Z", + ] + assert run1["ranges"][param_name]["values"] == [2, 4] + + for cov in covjson["coverages"]: + assert cov["mars:metadata"]["class"] == "ce" + assert cov["mars:metadata"]["stream"] == "efas" + # The raw date/time axes must not leak into metadata; the run reference + # is captured by "Forecast date" and the t-axis. + assert "date" not in cov["mars:metadata"] + assert "time" not in cov["mars:metadata"] + + +def test_efas_forecast_coverages_ordered_date_then_time(): + """A time-major tree (time above date) must still emit coverages ordered + date -> time -> point. + + We build the tree manually with the ``time`` axis *above* ``date`` so the + natural insertion/tree-traversal order is time-major + (May29 00Z, May30 00Z, May29 12Z, May30 12Z). The encoder must reorder the + output by the run reference (date + time) to + (May29 00Z, May29 12Z, May30 00Z, May30 12Z). + """ + dates = (np.datetime64("2026-05-29"), np.datetime64("2026-05-30")) + times = (np.timedelta64(0, "h"), np.timedelta64(12, "h")) + + # Single step so each run is a single value. Result is the flattened product + # over (time, date) in tree order: + # idx 0: time0 (00Z) date0 (May29) -> 10 + # idx 1: time0 (00Z) date1 (May30) -> 20 + # idx 2: time1 (12Z) date0 (May29) -> 30 + # idx 3: time1 (12Z) date1 (May30) -> 40 + tree = chain( + TensorIndexTree(), + node("class", ("ce",)), + node("time", times), + node("date", dates), + node("domain", ("g",)), + node("expver", ("8888",)), + node("levtype", ("sfc",)), + node("model", ("lisflood",)), + node("origin", ("ecmf",)), + node("param", ("240023",)), + node("step", (6,)), + node("stream", ("efas",)), + node("type", ("cf",)), + ) + tip(tree).add_child(make_point(50.0, 7.0, [10, 20, 30, 40])) + + covjson = Covjsonkit().encode("CoverageCollection", "timeseries").from_polytope_reforecast(tree) + + assert len(covjson["coverages"]) == 4 + + refs = [c["mars:metadata"]["Forecast date"] for c in covjson["coverages"]] + assert refs == [ + "2026-05-29T00:00:00Z", + "2026-05-29T12:00:00Z", + "2026-05-30T00:00:00Z", + "2026-05-30T12:00:00Z", + ] + + values = [next(iter(c["ranges"].values()))["values"][0] for c in covjson["coverages"]] + assert values == [10, 30, 20, 40] diff --git a/tests/test_encoder_separate_datetime_pointwise.py b/tests/test_encoder_separate_datetime_pointwise.py new file mode 100644 index 0000000..a03f2dc --- /dev/null +++ b/tests/test_encoder_separate_datetime_pointwise.py @@ -0,0 +1,120 @@ +"""Tests for the separate-datetime (independent date/hdate/time axes) reforecast +path in the point-wise encoders that override ``from_polytope_reforecast``: +Position (PointSeries) and VerticalProfile. + +These encoders don't use the composite/multi-point domain; instead they emit one +coverage per spatial point (Position: t across steps; VerticalProfile: levelist +across levels, single valid time per step). +""" + +import numpy as np +from conftest import ( + reforecast_separate_datetime_tree, + reforecast_separate_datetime_vertical_tree, +) + +from covjsonkit.api import Covjsonkit + +HDATES = (np.datetime64("2025-07-14T00:00:00"), np.datetime64("2025-07-15T00:00:00")) +TIMES = (np.timedelta64(0, "h"), np.timedelta64(12, "h")) + +# Two spatial points. Each result is the flattened product over (hdate, step, time). +# With a single step this is [h0t0, h0t1, h1t0, h1t1]. +POINTS = [ + (48.0, 11.0, [1.0, 2.0, 3.0, 4.0]), + (50.0, 12.0, [10.0, 20.0, 30.0, 40.0]), +] + +# (point lat/lon, reference datetime, value) +POSITION_EXPECTED = [ + ((48.0, 11.0), "2025-07-14T00:00:00Z", 1.0), + ((48.0, 11.0), "2025-07-14T12:00:00Z", 2.0), + ((48.0, 11.0), "2025-07-15T00:00:00Z", 3.0), + ((48.0, 11.0), "2025-07-15T12:00:00Z", 4.0), + ((50.0, 12.0), "2025-07-14T00:00:00Z", 10.0), + ((50.0, 12.0), "2025-07-14T12:00:00Z", 20.0), + ((50.0, 12.0), "2025-07-15T00:00:00Z", 30.0), + ((50.0, 12.0), "2025-07-15T12:00:00Z", 40.0), +] + + +def test_position_separate_datetime_single_step(): + tree = reforecast_separate_datetime_tree(POINTS, HDATES, TIMES) + covjson = Covjsonkit().encode("CoverageCollection", "position").from_polytope_reforecast(tree) + + assert covjson["domainType"] == "PointSeries" + assert len(covjson["coverages"]) == len(POSITION_EXPECTED) + + seen = set() + for cov in covjson["coverages"]: + axes = cov["domain"]["axes"] + lat = axes["latitude"]["values"][0] + lon = axes["longitude"]["values"][0] + assert axes["levelist"]["values"] == [0] + # Single step -> single valid time equal to the reference datetime. + assert axes["t"]["values"] == [cov["mars:metadata"]["Forecast date"]] + (param_name,) = cov["ranges"].keys() + (value,) = cov["ranges"][param_name]["values"] + assert "step" not in cov["mars:metadata"] + assert cov["mars:metadata"]["class"] == "ce" + seen.add(((lat, lon), cov["mars:metadata"]["Forecast date"], value)) + + assert seen == set(POSITION_EXPECTED) + + +def test_position_separate_datetime_multi_step_collapses_into_t(): + # 2 hdates x 2 times x 2 steps. result layout: product(hdate, step, time). + # idx: 0 h0s0t0, 1 h0s0t1, 2 h0s1t0, 3 h0s1t1, 4 h1s0t0, 5 h1s0t1, 6 h1s1t0, 7 h1s1t1 + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + tree = reforecast_separate_datetime_tree(points, HDATES, TIMES, step=(0, 6)) + covjson = Covjsonkit().encode("CoverageCollection", "position").from_polytope_reforecast(tree) + + # One coverage per (point, level, number, reference) = 4 references, single point. + assert len(covjson["coverages"]) == 4 + + by_ref = {c["mars:metadata"]["Forecast date"]: c for c in covjson["coverages"]} + # Reference 2025-07-14T00:00 (h0, t0): step0 -> idx0=1, step6 -> idx2=3. + cov = by_ref["2025-07-14T00:00:00Z"] + assert cov["domain"]["axes"]["t"]["values"] == [ + "2025-07-14T00:00:00Z", + "2025-07-14T06:00:00Z", + ] + (param_name,) = cov["ranges"].keys() + assert cov["ranges"][param_name]["values"] == [1, 3] + + +def test_verticalprofile_separate_datetime(): + # 2 levels, layout product(hdate, levelist, step, time) with single step: + # idx: 0 h0 L0 t0, 1 h0 L0 t1, 2 h0 L1 t0, 3 h0 L1 t1, + # 4 h1 L0 t0, 5 h1 L0 t1, 6 h1 L1 t0, 7 h1 L1 t1 + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + levels = (500, 850) + tree = reforecast_separate_datetime_vertical_tree(points, HDATES, TIMES, levels) + covjson = Covjsonkit().encode("CoverageCollection", "verticalprofile").from_polytope_reforecast(tree) + + assert covjson["domainType"] == "VerticalProfile" + # One coverage per (point, number, reference, step) = 4 references, single point/step. + assert len(covjson["coverages"]) == 4 + + by_ref = {c["mars:metadata"]["Forecast date"]: c for c in covjson["coverages"]} + cov = by_ref["2025-07-14T00:00:00Z"] + axes = cov["domain"]["axes"] + assert axes["levelist"]["values"] == [500, 850] + assert axes["t"]["values"] == ["2025-07-14T00:00:00Z"] + (param_name,) = cov["ranges"].keys() + # L0 t0 -> idx0=1, L1 t0 -> idx2=3 + assert cov["ranges"][param_name]["values"] == [1, 3] + assert cov["ranges"][param_name]["axisNames"] == ["levelist"] + assert cov["mars:metadata"]["step"] == 0 + + +def test_pointwise_reforecast_backward_compat_delegates(): + # A legacy merged tree without a separate ``time`` axis should delegate to + # from_polytope(date_key="hdate") and still produce PointSeries coverages. + from conftest import reforecast_branch, reforecast_tree + + branch = reforecast_branch(np.datetime64("2025-07-14T00:00:00"), POINTS) + tree = reforecast_tree([branch]) + covjson = Covjsonkit().encode("CoverageCollection", "position").from_polytope_reforecast(tree) + assert covjson["domainType"] == "PointSeries" + assert len(covjson["coverages"]) > 0 diff --git a/tests/test_encoder_separate_datetime_reforecast.py b/tests/test_encoder_separate_datetime_reforecast.py new file mode 100644 index 0000000..2679e84 --- /dev/null +++ b/tests/test_encoder_separate_datetime_reforecast.py @@ -0,0 +1,112 @@ +"""Tests for the separate-datetime (independent date/hdate/time axes) reforecast +path in the base encoder's ``from_polytope_reforecast``. + +Unlike the legacy merged representation (one coverage per hdate), separate-datetime +class=ce data has independent ``date``, ``hdate`` and ``time`` axes and produces one +coverage per ``(reference-datetime = hdate + time, step, number)``. +""" + +import numpy as np +import pytest +from conftest import forecast_separate_datetime_tree, reforecast_separate_datetime_tree + +from covjsonkit.api import Covjsonkit + +# Two spatial points. Each result is the flattened product over (hdate, step, time) +# with a single step -> [h0t0, h0t1, h1t0, h1t1]. +POINTS = [ + (48.0, 11.0, [1.0, 2.0, 3.0, 4.0]), + (50.0, 12.0, [10.0, 20.0, 30.0, 40.0]), +] +HDATES = (np.datetime64("2025-07-14T00:00:00"), np.datetime64("2025-07-15T00:00:00")) +TIMES = (np.timedelta64(0, "h"), np.timedelta64(12, "h")) + +EXPECTED = [ + ("2025-07-14T00:00:00Z", [1.0, 10.0]), + ("2025-07-14T12:00:00Z", [2.0, 20.0]), + ("2025-07-15T00:00:00Z", [3.0, 30.0]), + ("2025-07-15T12:00:00Z", [4.0, 40.0]), +] + + +@pytest.mark.parametrize("feature", ["BoundingBox", "Frame", "Circle", "Shapefile", "polygon"]) +def test_separate_datetime_reforecast_composite_encoders(feature): + tree = reforecast_separate_datetime_tree(POINTS, HDATES, TIMES) + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope_reforecast(tree) + + assert len(covjson["coverages"]) == len(EXPECTED) + for cov, (ref, vals) in zip(covjson["coverages"], EXPECTED): + assert cov["domain"]["axes"]["t"]["values"] == [ref] + assert cov["domain"]["axes"]["composite"]["values"] == [ + [48.0, 11.0, 0], + [50.0, 12.0, 0], + ] + # Single parameter (167 -> 2t). + (param_name,) = cov["ranges"].keys() + assert cov["ranges"][param_name]["values"] == vals + # efcl coverages are delineated by their valid time (the t-axis). + assert "Forecast date" not in cov["mars:metadata"] + assert cov["mars:metadata"]["step"] == 0 + assert cov["mars:metadata"]["class"] == "ce" + assert cov["mars:metadata"]["stream"] == "efcl" + + +def test_separate_datetime_reforecast_single_point_two_steps(): + # 2 hdates x 2 times x 2 steps = 8 coverages for a single point. + # result layout: product(hdate, step, time) -> h0s0t0, h0s0t1, h0s1t0, h0s1t1, ... + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + tree = reforecast_separate_datetime_tree(points, HDATES, TIMES, step=(0, 6)) + covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope_reforecast(tree) + + assert len(covjson["coverages"]) == 8 + # t = hdate + time + step; ordered by reference (hdate + time), then step. + seen = [ + (c["domain"]["axes"]["t"]["values"][0], c["mars:metadata"]["step"], c["ranges"]["2t"]["values"][0]) + for c in covjson["coverages"] + ] + assert seen == [ + ("2025-07-14T00:00:00Z", 0, 1.0), + ("2025-07-14T06:00:00Z", 6, 3.0), + ("2025-07-14T12:00:00Z", 0, 2.0), + ("2025-07-14T18:00:00Z", 6, 4.0), + ("2025-07-15T00:00:00Z", 0, 5.0), + ("2025-07-15T06:00:00Z", 6, 7.0), + ("2025-07-15T12:00:00Z", 0, 6.0), + ("2025-07-15T18:00:00Z", 6, 8.0), + ] + + +def test_separate_datetime_reforecast_polygon_metadata_date_is_iso(): + # Climatology data carries ``date`` as a pd.Timestamp; it must render ISO-8601. + import pandas as pd + + tree = reforecast_separate_datetime_tree(POINTS, HDATES, TIMES, date=pd.Timestamp("2023-01-01")) + covjson = Covjsonkit().encode("CoverageCollection", "polygon").from_polytope_reforecast(tree) + for cov in covjson["coverages"]: + assert cov["mars:metadata"]["date"] == "2023-01-01T00:00:00" + + +def test_efas_forecast_polygon_valid_time_and_run_reference(): + # No hdate axis: reference = date + time (the run), t = date + time + step. + # result layout: product(date, step, time) -> s6t0, s6t12, s24t0, s24t12. + tree = forecast_separate_datetime_tree([(50.0, 7.0, [1, 2, 3, 4]), (50.1, 7.1, [5, 6, 7, 8])], TIMES, step=(6, 24)) + covjson = Covjsonkit().encode("CoverageCollection", "polygon").from_polytope_reforecast(tree) + + seen = [ + ( + c["mars:metadata"]["Forecast date"], + c["mars:metadata"]["step"], + c["domain"]["axes"]["t"]["values"], + next(iter(c["ranges"].values()))["values"], + ) + for c in covjson["coverages"] + ] + assert seen == [ + ("2026-05-01T00:00:00Z", 6, ["2026-05-01T06:00:00Z"], [1.0, 5.0]), + ("2026-05-01T00:00:00Z", 24, ["2026-05-02T00:00:00Z"], [3.0, 7.0]), + ("2026-05-01T12:00:00Z", 6, ["2026-05-01T18:00:00Z"], [2.0, 6.0]), + ("2026-05-01T12:00:00Z", 24, ["2026-05-02T12:00:00Z"], [4.0, 8.0]), + ] + for cov in covjson["coverages"]: + # Raw date is folded into "Forecast date". + assert "date" not in cov["mars:metadata"] diff --git a/tests/test_encoder_time_series_from_polytope.py b/tests/test_encoder_time_series_from_polytope.py index 39e91db..dba4c45 100644 --- a/tests/test_encoder_time_series_from_polytope.py +++ b/tests/test_encoder_time_series_from_polytope.py @@ -523,4 +523,39 @@ def test_non_ce_efcl_hdate_not_collapsed(self): assert len(covjson["coverages"]) == 2 for cov in covjson["coverages"]: assert "Forecast date" in cov["mars:metadata"] + # Forecast date must be ISO-8601 with a trailing Z (e.g. 2025-07-14T06:00:00Z), + # not the pandas str() form "2025-07-14 06:00:00". + fd = cov["mars:metadata"]["Forecast date"] + assert "T" in fd and fd.endswith("Z") assert len(cov["domain"]["axes"]["t"]["values"]) == 1 + + def test_metadata_date_is_iso_compliant(self): + """A pd.Timestamp ``date`` axis (as climatology data yields) must render + in metadata as ISO-8601 (2023-01-01T00:00:00), not "2023-01-01 00:00:00".""" + import pandas as pd + + suffix = [ + ("domain", ("g",)), + ("expver", ("8888",)), + ("levtype", ("sfc",)), + ("model", ("lisflood",)), + ("origin", ("ecmf",)), + ("param", ("240023",)), + ("step", (6,)), + ("stream", ("efcl",)), + ("type", ("sfo",)), + ] + tree = chain( + TensorIndexTree(), + node("class", ("ce",)), + node("date", (pd.Timestamp("2023-01-01"),)), + node("hdate", (np.datetime64("2025-07-14T06:00:00"),)), + *[node(n, v) for n, v in suffix], + make_point(51.5, 6.5, [42.17]), + ) + + covjson = Covjsonkit().encode("CoverageCollection", "PointSeries").from_polytope_reforecast(tree) + + date_meta = covjson["coverages"][0]["mars:metadata"]["date"] + assert date_meta == "2023-01-01T00:00:00" + assert " " not in date_meta