From 7144ddc64bee79e50210cbf19bfc8c298016c5a0 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 3 Sep 2026 16:13:46 +0200 Subject: [PATCH 01/10] Add updated covjson output for vertical profile and position for encoder and decoder use --- covjsonkit/decoder/Position.py | 120 +++++++++------- covjsonkit/decoder/VerticalProfile.py | 35 ++--- covjsonkit/encoder/Position.py | 128 ++++++++---------- covjsonkit/encoder/VerticalProfile.py | 85 +++++------- tests/test_encoder_position_from_polytope.py | 68 ++++++---- ..._encoder_vertical_profile_from_polytope.py | 72 ++++++---- tests/test_spec_compliance.py | 99 ++++++++++++++ tests/test_xarray_vp.py | 4 +- 8 files changed, 367 insertions(+), 244 deletions(-) diff --git a/covjsonkit/decoder/Position.py b/covjsonkit/decoder/Position.py index 4331582..fa7cc34 100644 --- a/covjsonkit/decoder/Position.py +++ b/covjsonkit/decoder/Position.py @@ -9,18 +9,18 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() - if "x" in self.covjson["coverages"][0]["domain"]["axes"]: - self.x_name = "x" - else: - self.x_name = "latitude" - if "y" in self.covjson["coverages"][0]["domain"]["axes"]: - self.y_name = "y" - else: - self.y_name = "longitude" - if "z" in self.covjson["coverages"][0]["domain"]["axes"]: + first_axes = self.covjson["coverages"][0]["domain"]["axes"] + # Backwards-compatible axis-name detection: read both spec-compliant + # coverages (x/y/z) and legacy coverages (longitude/latitude/levelist). + # Semantics are preserved: x == longitude, y == latitude. + self.x_name = "x" if "x" in first_axes else "longitude" + self.y_name = "y" if "y" in first_axes else "latitude" + if "z" in first_axes: self.z_name = "z" - else: + elif "levelist" in first_axes: self.z_name = "levelist" + else: + self.z_name = None def get_domains(self): domains = [] @@ -48,9 +48,9 @@ def get_coordinates(self): coord_dict[param] = [] # Get x,y,z,t coords and unpack t coords and match to x,y,z coords for ind, domain in enumerate(self.domains): - x = domain["axes"][self.x_name]["values"][0] - y = domain["axes"][self.y_name]["values"][0] - z = domain["axes"][self.z_name]["values"][0] + longitude = domain["axes"][self.x_name]["values"][0] + latitude = domain["axes"][self.y_name]["values"][0] + level = domain["axes"][self.z_name]["values"][0] if self.z_name else None fct = domain["axes"]["t"]["values"][0] ts = domain["axes"]["t"]["values"] if "number" in self.mars_metadata[ind]: @@ -61,8 +61,7 @@ def get_coordinates(self): coords = [] for t in ts: # Have to replicate these coords for each parameter - # coordinates.append([x, y, z, t]) - coords.append([x, y, z, fct, t, num]) + coords.append([longitude, latitude, level, fct, t, num]) coord_dict[param].append(coords) return coord_dict @@ -75,13 +74,16 @@ def to_geotiff(self): def to_geojson(self): features = [] for coverage in self.covjson["coverages"]: - lat = coverage["domain"]["axes"][self.x_name]["values"][0] - lon = coverage["domain"]["axes"][self.y_name]["values"][0] - z = coverage["domain"]["axes"][self.z_name]["values"][0] + longitude = coverage["domain"]["axes"][self.x_name]["values"][0] + latitude = coverage["domain"]["axes"][self.y_name]["values"][0] datetimes = coverage["domain"]["axes"]["t"]["values"] if "mars:metadata" in coverage: mars_metadata = coverage["mars:metadata"] + geom_coords = [longitude, latitude] + if self.z_name: + geom_coords.append(coverage["domain"]["axes"][self.z_name]["values"][0]) + values = {} for key in coverage["ranges"]: values[key] = coverage["ranges"][key]["values"] @@ -96,7 +98,7 @@ def to_geojson(self): features.append( { "type": "Feature", - "geometry": {"type": "Point", "coordinates": [lon, lat, z]}, + "geometry": {"type": "Point", "coordinates": geom_coords}, "properties": param_vals, } ) @@ -106,7 +108,10 @@ def to_geojson(self): # function to convert covjson to xarray dataset def to_xarray(self): - dims = ["latitude", "longitude", "levelist", "number", "datetime", "t"] + if self.z_name: + dims = ["latitude", "longitude", "levelist", "number", "datetime", "t"] + else: + dims = ["latitude", "longitude", "number", "datetime", "t"] ds = [] unique_coords = set() # To track unique coordinate tuples @@ -114,13 +119,15 @@ def to_xarray(self): for domain in self.domains: # Extract coordinate values - x = domain["axes"][self.x_name]["values"][0] - y = domain["axes"][self.y_name]["values"][0] - z = domain["axes"][self.z_name]["values"][0] + longitude = domain["axes"][self.x_name]["values"][0] + latitude = domain["axes"][self.y_name]["values"][0] t = tuple(domain["axes"]["t"]["values"]) # Use tuple for hashable type - # Create a unique identifier for the domain - coord_tuple = (x, y, z, t) + if self.z_name: + z = domain["axes"][self.z_name]["values"][0] + coord_tuple = (longitude, latitude, z, t) + else: + coord_tuple = (longitude, latitude, t) # Check if this coordinate combination is already seen if coord_tuple not in unique_coords: @@ -140,23 +147,32 @@ def to_xarray(self): # Process each coordinate domain for coords in all_coords: dataarraydict = {} - x = coords["axes"][self.x_name]["values"] - y = coords["axes"][self.y_name]["values"] - z = coords["axes"][self.z_name]["values"] + longitude = coords["axes"][self.x_name]["values"] + latitude = coords["axes"][self.y_name]["values"] steps = coords["axes"]["t"]["values"] steps = [step.replace("Z", "") for step in steps] steps = pd.to_datetime(steps) - cov_idx_list = self._find_coverages(nums, datetime, x, y, z) - - coords = { - "latitude": x, - "longitude": y, - "levelist": z, - "number": nums, - "datetime": datetime, - "t": steps, - } + if self.z_name: + z = coords["axes"][self.z_name]["values"] + cov_idx_list = self._find_coverages(nums, datetime, longitude, latitude, z) + coord_dict = { + "latitude": latitude, + "longitude": longitude, + "levelist": z, + "number": nums, + "datetime": datetime, + "t": steps, + } + else: + cov_idx_list = self._find_coverages(nums, datetime, longitude, latitude, None) + coord_dict = { + "latitude": latitude, + "longitude": longitude, + "number": nums, + "datetime": datetime, + "t": steps, + } for parameter in self.parameters: param_values = [[[] for _ in range(len(datetime))] for _ in range(len(nums))] @@ -175,13 +191,12 @@ def to_xarray(self): "units": self.get_parameter_metadata(parameter)["unit"]["symbol"], "long_name": long_name, } - dataarraydict[long_name] = ( - dims, - [[[param_values]]], - attrs, - ) + if self.z_name: + dataarraydict[long_name] = (dims, [[[param_values]]], attrs) + else: + dataarraydict[long_name] = (dims, [[param_values]], attrs) - ds.append(xr.Dataset(data_vars=dataarraydict, coords=coords)) + ds.append(xr.Dataset(data_vars=dataarraydict, coords=coord_dict)) # Combine all DataArrays into a Dataset for mars_metadata in self.mars_metadata[0]: @@ -194,22 +209,25 @@ def to_xarray(self): return ds - def _find_coverages(self, nums, datetime, x, y, z): + def _find_coverages(self, nums, datetime, longitude, latitude, z): """Find coverages matching domain parameters and return with indices.""" result = [] for i, num in enumerate(nums): for j, date in enumerate(datetime): for coverage in self.covjson["coverages"]: - if self._covers_domain(coverage, num, date, x, y, z): + if self._covers_domain(coverage, num, date, longitude, latitude, z): result.append((i, j, coverage)) return result - def _covers_domain(self, coverage, num, date, x, y, z): + def _covers_domain(self, coverage, num, date, longitude, latitude, z): """Check if coverage matches the given domain parameters.""" - return ( + axes = coverage["domain"]["axes"] + match = ( coverage["mars:metadata"]["number"] == num and coverage["mars:metadata"]["Forecast date"] == date - and coverage["domain"]["axes"][self.x_name]["values"] == x - and coverage["domain"]["axes"][self.y_name]["values"] == y - and coverage["domain"]["axes"][self.z_name]["values"] == z + and axes[self.x_name]["values"] == longitude + and axes[self.y_name]["values"] == latitude ) + if match and self.z_name and z is not None: + match = axes[self.z_name]["values"] == z + return match diff --git a/covjsonkit/decoder/VerticalProfile.py b/covjsonkit/decoder/VerticalProfile.py index 23dc1f4..591a588 100644 --- a/covjsonkit/decoder/VerticalProfile.py +++ b/covjsonkit/decoder/VerticalProfile.py @@ -12,18 +12,13 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() - if "x" in self.covjson["coverages"][0]["domain"]["axes"]: - self.x_name = "x" - else: - self.x_name = "latitude" - if "y" in self.covjson["coverages"][0]["domain"]["axes"]: - self.y_name = "y" - else: - self.y_name = "longitude" - if "z" in self.covjson["coverages"][0]["domain"]["axes"]: - self.z_name = "z" - else: - self.z_name = "levelist" + # Backwards-compatible axis-name detection: read both spec-compliant + # coverages (x/y/z) and legacy coverages (longitude/latitude/levelist). + # Semantics are preserved: x == longitude, y == latitude. + first_axes = self.covjson["coverages"][0]["domain"]["axes"] + self.x_name = "x" if "x" in first_axes else "longitude" + self.y_name = "y" if "y" in first_axes else "latitude" + self.z_name = "z" if "z" in first_axes else "levelist" def get_domains(self): domains = [] @@ -74,8 +69,8 @@ def to_geotiff(self): def to_geojson(self): features = [] for coverage in self.covjson["coverages"]: - lat = coverage["domain"]["axes"][self.x_name]["values"][0] - lon = coverage["domain"]["axes"][self.y_name]["values"][0] + lon = coverage["domain"]["axes"][self.x_name]["values"][0] + lat = coverage["domain"]["axes"][self.y_name]["values"][0] levels = coverage["domain"]["axes"][self.z_name]["values"] datetimes = coverage["domain"]["axes"]["t"]["values"] if "mars:metadata" in coverage: @@ -143,8 +138,8 @@ def to_xarray(self): dataarraydict = {} # Get coordinates - x = coords["axes"][self.x_name]["values"] - y = coords["axes"][self.y_name]["values"] + longitude = coords["axes"][self.x_name]["values"] + latitude = coords["axes"][self.y_name]["values"] levelist = coords["axes"][self.z_name]["values"] steps = coords["axes"]["t"]["values"] @@ -182,8 +177,8 @@ def to_xarray(self): if ( coverage["mars:metadata"]["number"] == num and coverage["mars:metadata"]["Forecast date"] == date - and coverage["domain"]["axes"][self.x_name]["values"][0] == x[0] - and coverage["domain"]["axes"][self.y_name]["values"][0] == y[0] + and coverage["domain"]["axes"][self.x_name]["values"][0] == longitude[0] + and coverage["domain"]["axes"][self.y_name]["values"][0] == latitude[0] and coverage["domain"]["axes"]["t"]["values"][0] == new_step ): param_values[parameter][domain_idx][i][j][k] = coverage["ranges"][parameter][ @@ -192,8 +187,8 @@ def to_xarray(self): for parameter in self.parameters: param_coords = { - "latitude": x, - "longitude": y, + "latitude": latitude, + "longitude": longitude, "number": nums, "datetime": datetime, "time": steps, diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index 26fef03..2068ac7 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -13,30 +13,28 @@ def __init__(self, type, domaintype): self.covjson["domainType"] = "PointSeries" self.covjson["coverages"] = [] - def add_coverage(self, mars_metadata, coords, values): + def add_coverage(self, mars_metadata, coords, values, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) # cov = Coverage.model_validate_json(json.dumps(new_coverage)) # self.pydantic_coverage.coverages.append(cov) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): coverage["domain"]["type"] = "Domain" - coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["latitude"] = {} - coverage["domain"]["axes"]["longitude"] = {} - coverage["domain"]["axes"]["levelist"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["latitude"]["values"] = coords["latitude"] - coverage["domain"]["axes"]["longitude"]["values"] = coords["longitude"] - coverage["domain"]["axes"]["levelist"]["values"] = coords["levelist"] - coverage["domain"]["axes"]["t"]["values"] = coords["t"] + axes = {} + axes["x"] = {"values": coords["longitude"]} + axes["y"] = {"values": coords["latitude"]} + if include_z: + axes["z"] = {"values": coords["levelist"]} + axes["t"] = {"values": coords["t"]} + coverage["domain"]["axes"] = axes def add_range(self, coverage, values): for parameter in values.keys(): @@ -45,7 +43,7 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] + coverage["ranges"][param]["axisNames"] = ["t"] coverage["ranges"][param]["values"] = values[ parameter ] # [values[parameter][val][0] for val in values[parameter].keys()] @@ -53,6 +51,29 @@ def add_range(self, coverage, values): def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, datasets): """ Converts an xarray dataset or a list of xarray datasets into an OGC CoverageJSON @@ -71,27 +92,10 @@ def from_xarray(self, datasets): self.covjson["domainType"] = "PointSeries" self.covjson["coverages"] = [] - if "latitude" in datasets[0].coords: - x_coord = "latitude" - elif "x" in datasets[0].coords: - x_coord = "x" - if "longitude" in datasets[0].coords: - y_coord = "longitude" - elif "y" in datasets[0].coords: - y_coord = "y" - if "levelist" in datasets[0].coords: - z_coord = "levelist" + include_z = "levelist" in datasets[0].coords # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + self._set_references(include_z) for data_var in datasets[0].data_vars: data_var = self.convert_param_to_param_id(data_var) @@ -115,10 +119,11 @@ def from_xarray(self, datasets): { "latitude": [float(x) for x in dataset["latitude"].values], "longitude": [float(x) for x in dataset["longitude"].values], - "levelist": [float(x) for x in dataset["levelist"].values], + "levelist": [float(x) for x in dataset["levelist"].values] if include_z else None, "t": [str(x) for x in dataset["t"].values], }, dv_dict, + include_z=include_z, ) return self.covjson @@ -143,6 +148,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = 0 fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 @@ -155,15 +161,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) coordinates = {} @@ -246,8 +245,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["Forecast date"] = date + mm["levelist"] = level + coordinates[date][i]["levelist"] = [level] del mm["step"] - self.add_coverage(mm, coordinates[date][i], val_dict) + self.add_coverage(mm, coordinates[date][i], val_dict, include_z) end = time.time() delta = end - start @@ -268,6 +269,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 @@ -279,15 +281,8 @@ def from_polytope_month(self, result): start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) if fields["param"] == 0: raise ValueError("No data was returned.") @@ -344,7 +339,7 @@ def from_polytope_month(self, result): mm = mars_metadata.copy() mm["number"] = num mm["levelist"] = level - self.add_coverage(mm, coord_entry, val_dict) + self.add_coverage(mm, coord_entry, val_dict, include_z) end = time.time() logging.debug("Coverage creation: %s", end) # noqa: E501 @@ -364,6 +359,7 @@ def from_polytope_step(self, result): fields["dates"] = [] fields["levels"] = [0] fields["times"] = [] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 @@ -376,15 +372,8 @@ def from_polytope_step(self, result): start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) coordinates = {} @@ -413,11 +402,6 @@ def from_polytope_step(self, result): for para in fields["param"]: for date in fields["dates"]: for times in fields["times"]: - # date_format = "%Y%m%dT%H%M%S" - # new_date = pd.Timestamp(date).strftime(date_format) - # start_time = datetime.strptime(new_date, date_format) - # add current date to list by converting it to iso format - # stamp = start_time + timedelta(hours=int(step)) datetime = pd.Timestamp(date) + times coordinates[fields["dates"][0]][(i * len(fields["levels"]) + j)]["t"].append( str(datetime).split("+")[0] + "Z" @@ -441,14 +425,16 @@ def from_polytope_step(self, result): val_dict[para] = [] for date in fields["dates"]: key = (date, level, num, para) - # for k, v in range_dict.items(): - # if k == key: - # val_dict[para].append(v[0]) val_dict[para].extend(range_dict[key][i]) mm = mars_metadata.copy() mm["number"] = num mm["Forecast date"] = date - self.add_coverage(mm, coordinates[fields["dates"][0]][(i * len(fields["levels"]) + j)], val_dict) + self.add_coverage( + mm, + coordinates[fields["dates"][0]][(i * len(fields["levels"]) + j)], + val_dict, + include_z, + ) end = time.time() delta = end - start diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 0781f89..2029570 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -25,16 +25,15 @@ def add_coverage(self, mars_metadata, coords, values): self.covjson["coverages"].append(new_coverage) def add_domain(self, coverage, coords): + # OGC-compliant VerticalProfile domain: x == longitude, y == latitude, + # z == the vertical (levelist) profile axis, t == a single time value. coverage["domain"]["type"] = "Domain" - coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["latitude"] = {} - coverage["domain"]["axes"]["longitude"] = {} - coverage["domain"]["axes"]["levelist"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["latitude"]["values"] = coords["latitude"] - coverage["domain"]["axes"]["longitude"]["values"] = coords["longitude"] - coverage["domain"]["axes"]["levelist"]["values"] = coords["levelist"] - coverage["domain"]["axes"]["t"]["values"] = coords["t"] + coverage["domain"]["axes"] = { + "x": {"values": coords["longitude"]}, + "y": {"values": coords["latitude"]}, + "z": {"values": coords["levelist"]}, + "t": {"values": coords["t"]}, + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -43,12 +42,33 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["levelist"] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["z"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self): + # VerticalProfile always carries a vertical (z) axis, so referencing is + # split into GeographicCRS (x, y), VerticalCRS (z) and TemporalRS (t). + self.covjson["referencing"] = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + def from_xarray(self, datasets): """ Converts an xarray dataset or a list of xarray datasets into an OGC CoverageJSON @@ -64,30 +84,11 @@ def from_xarray(self, datasets): datasets = [datasets] self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "VeticalProfile" + self.covjson["domainType"] = "VerticalProfile" self.covjson["coverages"] = [] - if "latitude" in datasets[0].coords: - x_coord = "latitude" - elif "x" in datasets[0].coords: - x_coord = "x" - if "longitude" in datasets[0].coords: - y_coord = "longitude" - elif "y" in datasets[0].coords: - y_coord = "y" - if "levelist" in datasets[0].coords: - z_coord = "levelist" - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + self._set_references() for data_var in datasets[0].data_vars: data_var = self.convert_param_to_param_id(data_var) @@ -145,15 +146,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + self._set_references() coordinates = {} @@ -263,15 +256,7 @@ def from_polytope_month(self, result): start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + self._set_references() levels = fields["levels"] if fields["param"] == 0: diff --git a/tests/test_encoder_position_from_polytope.py b/tests/test_encoder_position_from_polytope.py index 0c14a3f..c21d829 100644 --- a/tests/test_encoder_position_from_polytope.py +++ b/tests/test_encoder_position_from_polytope.py @@ -1,5 +1,10 @@ import numpy as np -from conftest import forecast_tree, reforecast_branch, reforecast_tree +from conftest import ( + assert_valid_covjson, + forecast_tree, + reforecast_branch, + reforecast_tree, +) from covjsonkit.api import Covjsonkit @@ -14,11 +19,22 @@ def test_single_point_two_steps(self): assert covjson["type"] == "CoverageCollection" assert covjson["domainType"] == "PointSeries" - - # Referencing - ref = covjson["referencing"][0] - assert ref["coordinates"] == ["latitude", "longitude", "levelist"] - assert ref["system"]["type"] == "GeographicCRS" + assert_valid_covjson(covjson) + + # Spec-compliant split referencing (surface -> no z). + assert covjson["referencing"] == [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] # Parameters assert "2t" in covjson["parameters"] @@ -28,18 +44,18 @@ def test_single_point_two_steps(self): cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [0]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, "t": {"values": ["2025-01-01T00:00:00Z", "2025-01-01T06:00:00Z"]}, } + assert "z" not in cov["domain"]["axes"] assert cov["ranges"] == { "2t": { "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["t"], "values": [264.9, 263.8], } } @@ -53,6 +69,7 @@ def test_single_point_two_steps(self): "stream": "oper", "type": "fc", "number": 0, + "levelist": 0, } def test_two_points_two_steps(self): @@ -62,6 +79,7 @@ def test_two_points_two_steps(self): step=(0, 6), ) covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope(tree) + assert_valid_covjson(covjson) shared_metadata = { "class": "od", @@ -72,6 +90,7 @@ def test_two_points_two_steps(self): "stream": "oper", "type": "fc", "number": 0, + "levelist": 0, } expected = [ @@ -81,9 +100,8 @@ def test_two_points_two_steps(self): assert len(covjson["coverages"]) == len(expected) for cov, (lat, lon, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "latitude": {"values": [lat]}, - "longitude": {"values": [lon]}, - "levelist": {"values": [0]}, + "x": {"values": [lon]}, + "y": {"values": [lat]}, "t": {"values": ["2025-01-01T00:00:00Z", "2025-01-01T06:00:00Z"]}, } assert cov["ranges"]["2t"]["values"] == vals @@ -93,14 +111,14 @@ def test_single_step(self): """Edge case: 1 step -> shape [1], single t-value.""" tree = forecast_tree([(48.0, 11.0, [264.9])], step=(0,)) covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope(tree) + assert_valid_covjson(covjson) assert len(covjson["coverages"]) == 1 cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [0]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, "t": {"values": ["2025-01-01T00:00:00Z"]}, } assert cov["ranges"] == { @@ -108,7 +126,7 @@ def test_single_step(self): "type": "NdArray", "dataType": "float", "shape": [1], - "axisNames": ["2t"], + "axisNames": ["t"], "values": [264.9], } } @@ -127,6 +145,7 @@ def test_reforecast_single_hdate_two_points(self): ) covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope_reforecast(tree) + assert_valid_covjson(covjson) shared_metadata = { "class": "ce", @@ -138,6 +157,7 @@ def test_reforecast_single_hdate_two_points(self): "stream": "efcl", "type": "sfo", "number": 0, + "levelist": 0, } expected = [ @@ -147,16 +167,15 @@ def test_reforecast_single_hdate_two_points(self): assert len(covjson["coverages"]) == len(expected) for cov, (lat, lon, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "latitude": {"values": [lat]}, - "longitude": {"values": [lon]}, - "levelist": {"values": [0]}, + "x": {"values": [lon]}, + "y": {"values": [lat]}, "t": {"values": ["2025-07-14T06:00:00Z"]}, } assert cov["ranges"]["2t"]["values"] == vals assert cov["mars:metadata"] == shared_metadata def test_reforecast_two_hdates_two_points(self): - """2 hdates x 2 points -> 4 coverages.""" + """2 hdates x 2 points -> 4 coverages (Position does not collapse reanalysis).""" tree = reforecast_tree( [ reforecast_branch(np.datetime64("2025-07-14T06:00:00"), [(48.0, 11.0, [264.9]), (50.0, 12.0, [265.1])]), @@ -165,6 +184,7 @@ def test_reforecast_two_hdates_two_points(self): ) covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope_reforecast(tree) + assert_valid_covjson(covjson) shared_metadata = { "class": "ce", @@ -175,6 +195,7 @@ def test_reforecast_two_hdates_two_points(self): "stream": "efcl", "type": "sfo", "number": 0, + "levelist": 0, } expected = [ @@ -186,9 +207,8 @@ def test_reforecast_two_hdates_two_points(self): assert len(covjson["coverages"]) == len(expected) for cov, (lat, lon, t, vals, fc_date) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "latitude": {"values": [lat]}, - "longitude": {"values": [lon]}, - "levelist": {"values": [0]}, + "x": {"values": [lon]}, + "y": {"values": [lat]}, "t": {"values": t}, } assert cov["ranges"]["2t"]["values"] == vals diff --git a/tests/test_encoder_vertical_profile_from_polytope.py b/tests/test_encoder_vertical_profile_from_polytope.py index 7380fd4..a443af0 100644 --- a/tests/test_encoder_vertical_profile_from_polytope.py +++ b/tests/test_encoder_vertical_profile_from_polytope.py @@ -1,5 +1,5 @@ import numpy as np -from conftest import chain, make_leaf, make_point, node, tip +from conftest import assert_valid_covjson, chain, make_leaf, make_point, node, tip from polytope_feature.datacube.tensor_index_tree import TensorIndexTree from covjsonkit.api import Covjsonkit @@ -44,11 +44,27 @@ def test_single_point_three_levels(self): assert covjson["type"] == "CoverageCollection" assert covjson["domainType"] == "VerticalProfile" - - # Referencing (folded from removed test_referencing) - ref = covjson["referencing"][0] - assert ref["coordinates"] == ["latitude", "longitude", "levelist"] - assert ref["system"]["type"] == "GeographicCRS" + assert_valid_covjson(covjson) + + # Spec-compliant split referencing: GeographicCRS (x, y), VerticalCRS (z), + # TemporalRS (t). + assert covjson["referencing"] == [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] # Parameters (folded from removed test_parameters_block) assert "t" in covjson["parameters"] @@ -59,9 +75,9 @@ def test_single_point_three_levels(self): cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [1000, 850, 500]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, + "z": {"values": [1000, 850, 500]}, "t": {"values": ["2025-01-01T00:00:00Z"]}, } @@ -70,7 +86,7 @@ def test_single_point_three_levels(self): "type": "NdArray", "dataType": "float", "shape": [3], - "axisNames": ["levelist"], + "axisNames": ["z"], "values": [290.1, 280.2, 250.3], } } @@ -109,6 +125,7 @@ def test_two_points_two_levels(self): lev_node.add_child(make_point(50.0, 13.0, [288.5, 248.7])) covjson = Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope(tree) + assert_valid_covjson(covjson) shared_metadata = { "class": "od", @@ -130,9 +147,9 @@ def test_two_points_two_levels(self): assert len(covjson["coverages"]) == len(expected) for cov, (lat, lon, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "latitude": {"values": [lat]}, - "longitude": {"values": [lon]}, - "levelist": {"values": [1000, 500]}, + "x": {"values": [lon]}, + "y": {"values": [lat]}, + "z": {"values": [1000, 500]}, "t": {"values": ["2025-06-15T12:00:00Z"]}, } assert cov["ranges"] == { @@ -140,7 +157,7 @@ def test_two_points_two_levels(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["levelist"], + "axisNames": ["z"], "values": vals, } } @@ -150,14 +167,15 @@ def test_step_offset(self): """Step=6 should shift the t coordinate by 6 hours.""" tree = self._build_vp_tree(step=6, levels_values={1000: 290.0}) covjson = Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope(tree) + assert_valid_covjson(covjson) assert len(covjson["coverages"]) == 1 cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [1000]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, + "z": {"values": [1000]}, "t": {"values": ["2025-01-01T06:00:00Z"]}, } @@ -166,7 +184,7 @@ def test_step_offset(self): "type": "NdArray", "dataType": "float", "shape": [1], - "axisNames": ["levelist"], + "axisNames": ["z"], "values": [290.0], } } @@ -227,14 +245,15 @@ def test_reforecast_single_hdate_three_levels(self): assert covjson["type"] == "CoverageCollection" assert covjson["domainType"] == "VerticalProfile" + assert_valid_covjson(covjson) assert len(covjson["coverages"]) == 1 cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [1000, 850, 500]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, + "z": {"values": [1000, 850, 500]}, "t": {"values": ["2025-07-14T12:00:00Z"]}, } @@ -243,7 +262,7 @@ def test_reforecast_single_hdate_three_levels(self): "type": "NdArray", "dataType": "float", "shape": [3], - "axisNames": ["levelist"], + "axisNames": ["z"], "values": [290.1, 280.2, 250.3], } } @@ -276,6 +295,7 @@ def test_reforecast_two_hdates_three_levels(self): date_node.add_child(branch) covjson = Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope_reforecast(tree) + assert_valid_covjson(covjson) expected = [ ("2025-07-14T12:00:00Z", [290.1, 280.2, 250.3], "2025-07-14T06:00:00Z"), @@ -284,9 +304,9 @@ def test_reforecast_two_hdates_three_levels(self): assert len(covjson["coverages"]) == len(expected) for cov, (t, vals, fc_date) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [1000, 850, 500]}, + "x": {"values": [11.0]}, + "y": {"values": [48.0]}, + "z": {"values": [1000, 850, 500]}, "t": {"values": [t]}, } assert cov["ranges"] == { @@ -294,7 +314,7 @@ def test_reforecast_two_hdates_three_levels(self): "type": "NdArray", "dataType": "float", "shape": [3], - "axisNames": ["levelist"], + "axisNames": ["z"], "values": vals, } } diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index 3782c5d..bab0373 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -91,6 +91,23 @@ def _step_tree(): ) +def _vertical_profile_tree(): + return chain( + TensorIndexTree(), + node("class", ("od",)), + node("date", (np.datetime64("2025-01-01T00:00:00"),)), + node("domain", ("g",)), + node("expver", ("0001",)), + node("levtype", ("pl",)), + node("param", ("130",)), + node("step", (0,)), + node("stream", ("oper",)), + node("type", ("an",)), + node("levelist", (1000, 850, 500)), + make_point(48.0, 11.0, [290.1, 280.2, 250.3]), + ) + + class TestSpecCompliantOutputValidates: """Encoder output must validate against covjson-pydantic (test-only guard).""" @@ -180,3 +197,85 @@ def test_legacy_month_xarray_equivalent(self): legacy_ds = Covjsonkit().decode(legacy_covjson).to_xarray() assert new_ds.identical(legacy_ds) + + +class TestPositionSpecCompliance: + """Position (PointSeries) encoder output validates and decodes with back-compat.""" + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_position_output_validates(self, tree_factory): + covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope(tree_factory()) + assert_valid_covjson(covjson) + + def test_position_month_validates(self): + covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope_month(_month_tree()) + assert_valid_covjson(covjson) + + def test_position_step_validates(self): + covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope_step(_step_tree()) + assert_valid_covjson(covjson) + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_position_legacy_and_new_geojson_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope(tree_factory()) + legacy_covjson = _to_legacy(new_covjson) + + new_gj = Covjsonkit().decode(new_covjson).to_geojson() + legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() + + assert legacy_gj == new_gj + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_position_legacy_and_new_xarray_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "Position").from_polytope(tree_factory()) + legacy_covjson = _to_legacy(new_covjson) + + new_ds = Covjsonkit().decode(new_covjson).to_xarray() + legacy_ds = Covjsonkit().decode(legacy_covjson).to_xarray() + + assert new_ds.identical(legacy_ds) + + assert "latitude" in legacy_ds.coords + assert "longitude" in legacy_ds.coords + has_z = any("z" in cov["domain"]["axes"] for cov in new_covjson["coverages"]) + assert ("levelist" in legacy_ds.coords) == has_z + + +class TestVerticalProfileSpecCompliance: + """VerticalProfile encoder output validates and decodes with back-compat.""" + + def test_vertical_profile_output_validates(self): + covjson = Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope(_vertical_profile_tree()) + assert_valid_covjson(covjson) + + def test_vertical_profile_legacy_and_new_geojson_equivalent(self): + new_covjson = ( + Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope(_vertical_profile_tree()) + ) + legacy_covjson = _to_legacy(new_covjson) + + new_gj = Covjsonkit().decode(new_covjson).to_geojson() + legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() + + assert legacy_gj == new_gj + + # The vertical (z) axis must survive into the GeoJSON geometry as the + # third coordinate, and lon/lat must be in the correct [lon, lat, z] order. + first = new_gj["features"][0]["geometry"]["coordinates"] + assert first == [11.0, 48.0, 1000] + + def test_vertical_profile_legacy_and_new_xarray_equivalent(self): + new_covjson = ( + Covjsonkit().encode("CoverageCollection", "VerticalProfile").from_polytope(_vertical_profile_tree()) + ) + legacy_covjson = _to_legacy(new_covjson) + + new_ds = Covjsonkit().decode(new_covjson).to_xarray() + legacy_ds = Covjsonkit().decode(legacy_covjson).to_xarray() + + assert new_ds.identical(legacy_ds) + + # Decoder read legacy coverage by axis NAME: lat/lon/levelist coords present. + assert "latitude" in legacy_ds.coords + assert "longitude" in legacy_ds.coords + assert "levelist" in legacy_ds.coords diff --git a/tests/test_xarray_vp.py b/tests/test_xarray_vp.py index d9c31b0..c6ce031 100644 --- a/tests/test_xarray_vp.py +++ b/tests/test_xarray_vp.py @@ -34,11 +34,11 @@ def test_from_xarray(self): assert covjson_result["type"] == self.test_covjson["type"] assert len(covjson_result["coverages"]) == len(self.test_covjson["coverages"]) assert ( - covjson_result["coverages"][0]["domain"]["axes"]["latitude"]["values"][0] + covjson_result["coverages"][0]["domain"]["axes"]["y"]["values"][0] == self.test_covjson["coverages"][0]["domain"]["axes"]["latitude"]["values"][0] ) assert ( - covjson_result["coverages"][0]["domain"]["axes"]["longitude"]["values"][0] + covjson_result["coverages"][0]["domain"]["axes"]["x"]["values"][0] == self.test_covjson["coverages"][0]["domain"]["axes"]["longitude"]["values"][0] ) assert ( From a2a3340060eb4a2dfa1fd3328548f6c09911b8ba Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 4 Sep 2026 11:07:42 +0200 Subject: [PATCH 02/10] Update composite axis features (bbox, frame, shapefile) to align with covjson spec --- covjsonkit/decoder/BoundingBox.py | 64 +++-- covjsonkit/decoder/Circle.py | 57 ++-- covjsonkit/decoder/Frame.py | 25 +- covjsonkit/decoder/Shapefile.py | 25 +- covjsonkit/decoder/Wkt.py | 25 +- covjsonkit/encoder/BoundingBox.py | 246 ++++++------------ covjsonkit/encoder/Circle.py | 142 +++++----- covjsonkit/encoder/Frame.py | 143 +++++----- covjsonkit/encoder/Shapefile.py | 144 +++++----- covjsonkit/encoder/Wkt.py | 163 ++++++------ tests/conftest.py | 6 +- ...test_encoder_bounding_box_from_polytope.py | 13 +- tests/test_encoder_circle_from_polytope.py | 12 +- tests/test_encoder_frame_from_polytope.py | 12 +- tests/test_encoder_shapefile_from_polytope.py | 12 +- tests/test_encoder_wkt_from_polytope.py | 12 +- tests/test_spec_compliance.py | 138 ++++++++++ tests/test_xarray_multipoint.py | 13 +- 18 files changed, 689 insertions(+), 563 deletions(-) diff --git a/covjsonkit/decoder/BoundingBox.py b/covjsonkit/decoder/BoundingBox.py index 34e0deb..0e91430 100644 --- a/covjsonkit/decoder/BoundingBox.py +++ b/covjsonkit/decoder/BoundingBox.py @@ -17,6 +17,25 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering: derive the position of the + # longitude (x), latitude (y) and optional vertical (z) components from + # the ``composite`` axis's ``coordinates`` labels. This reads both + # spec-compliant coverages (labels x/y[/z], values [lon, lat[, level]]) + # and legacy coverages (labels latitude/longitude/levelist, values + # [lat, lon[, level]]). + labels = self.domains[0]["axes"]["composite"]["coordinates"] + self.x_idx = self._label_index(labels, ("x", "longitude")) + self.y_idx = self._label_index(labels, ("y", "latitude")) + self.z_idx = self._label_index(labels, ("z", "levelist"), required=False) + + @staticmethod + def _label_index(labels, names, required=True): + for name in names: + if name in labels: + return labels.index(name) + if required: + raise ValueError(f"None of {names} found in composite coordinates {labels}") + return None def get_domains(self): domains = [] @@ -51,9 +70,9 @@ def to_geotiff(self, output_file="multipoint", resolution=0.01): if rasterio is None: raise ImportError("Please install 'rasterio' to use this feature: pip install covjsonkit[geo]") coords = self.covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] - x = [c[1] for c in coords] # longitude - y = [c[0] for c in coords] # latitude - # z = [c[2] for c in coords] # height/time/etc (not used yet) + x = [c[self.x_idx] for c in coords] # longitude + y = [c[self.y_idx] for c in coords] # latitude + # z = [c[self.z_idx] for c in coords] # height/time/etc (not used yet) # Define grid x_min, x_max = min(x), max(x) @@ -119,10 +138,13 @@ def to_geojson(self): param_vals["datetime"] = datetime if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata + geom_coords = [lonlat[self.x_idx], lonlat[self.y_idx]] + if self.z_idx is not None: + geom_coords.append(lonlat[self.z_idx]) features.append( { "type": "Feature", - "geometry": {"type": "Point", "coordinates": [lonlat[1], lonlat[0], lonlat[2]]}, + "geometry": {"type": "Point", "coordinates": geom_coords}, "properties": param_vals, } ) @@ -135,14 +157,15 @@ def to_xarray(self): dataarraydict = {} # Get coordinates - x = [] - y = [] - z = [] + longitude = [] + latitude = [] + levelist = [] datetimes = [] for coord in self.get_coordinates()["composite"]["values"]: - x.append(float(coord[0])) - y.append(float(coord[1])) - z.append(float(coord[2])) + longitude.append(float(coord[self.x_idx])) + latitude.append(float(coord[self.y_idx])) + if self.z_idx is not None: + levelist.append(float(coord[self.z_idx])) for datetime in self.get_coordinates()["t"]["values"]: datetimes.append(datetime) @@ -197,17 +220,20 @@ def to_xarray(self): dataarray.attrs["long_name"] = self.get_parameter_metadata(parameter)["observedProperty"]["id"] dataarraydict[dataarray.attrs["long_name"]] = dataarray + coords_dict = dict( + datetimes=(["datetimes"], datetimes), + number=(["number"], numbers), + steps=(["steps"], steps), + points=(["points"], list(range(0, len(longitude)))), + latitude=(["points"], latitude), + longitude=(["points"], longitude), + ) + if self.z_idx is not None: + coords_dict["levelist"] = (["points"], levelist) + ds = xr.Dataset( dataarraydict, - coords=dict( - datetimes=(["datetimes"], datetimes), - number=(["number"], numbers), - steps=(["steps"], steps), - points=(["points"], list(range(0, len(x)))), - latitude=(["points"], x), - longitude=(["points"], y), - levelist=(["points"], z), - ), + coords=coords_dict, ) for mars_metadata in self.mars_metadata[0]: ds.attrs[mars_metadata] = self.mars_metadata[0][mars_metadata] diff --git a/covjsonkit/decoder/Circle.py b/covjsonkit/decoder/Circle.py index a3a1ae9..adddfca 100644 --- a/covjsonkit/decoder/Circle.py +++ b/covjsonkit/decoder/Circle.py @@ -9,6 +9,24 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering: derive the position of the + # longitude (x), latitude (y) and optional vertical (z) components from + # the ``composite`` axis's ``coordinates`` labels. Reads both spec + # coverages (x/y[/z], values [lon, lat[, level]]) and legacy coverages + # (latitude/longitude/levelist, values [lat, lon[, level]]). + labels = self.domains[0]["axes"]["composite"]["coordinates"] + self.x_idx = self._label_index(labels, ("x", "longitude")) + self.y_idx = self._label_index(labels, ("y", "latitude")) + self.z_idx = self._label_index(labels, ("z", "levelist"), required=False) + + @staticmethod + def _label_index(labels, names, required=True): + for name in names: + if name in labels: + return labels.index(name) + if required: + raise ValueError(f"None of {names} found in composite coordinates {labels}") + return None def get_domains(self): domains = [] @@ -61,10 +79,13 @@ def to_geojson(self): param_vals["datetime"] = datetime if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata + geom_coords = [lonlat[self.x_idx], lonlat[self.y_idx]] + if self.z_idx is not None: + geom_coords.append(lonlat[self.z_idx]) features.append( { "type": "Feature", - "geometry": {"type": "Point", "coordinates": [lonlat[1], lonlat[0], lonlat[2]]}, + "geometry": {"type": "Point", "coordinates": geom_coords}, "properties": param_vals, } ) @@ -77,14 +98,15 @@ def to_xarray(self): dataarraydict = {} # Get coordinates - x = [] - y = [] - z = [] + longitude = [] + latitude = [] + levelist = [] datetimes = [] for coord in self.get_coordinates()["composite"]["values"]: - x.append(float(coord[0])) - y.append(float(coord[1])) - z.append(float(coord[2])) + longitude.append(float(coord[self.x_idx])) + latitude.append(float(coord[self.y_idx])) + if self.z_idx is not None: + levelist.append(float(coord[self.z_idx])) for datetime in self.get_coordinates()["t"]["values"]: datetimes.append(datetime) @@ -139,17 +161,20 @@ def to_xarray(self): dataarray.attrs["long_name"] = self.get_parameter_metadata(parameter)["observedProperty"]["id"] dataarraydict[dataarray.attrs["long_name"]] = dataarray + coords_dict = dict( + datetimes=(["datetimes"], datetimes), + number=(["number"], numbers), + steps=(["steps"], steps), + points=(["points"], list(range(0, len(longitude)))), + latitude=(["points"], latitude), + longitude=(["points"], longitude), + ) + if self.z_idx is not None: + coords_dict["levelist"] = (["points"], levelist) + ds = xr.Dataset( dataarraydict, - coords=dict( - datetimes=(["datetimes"], datetimes), - number=(["number"], numbers), - steps=(["steps"], steps), - points=(["points"], list(range(0, len(x)))), - latitude=(["points"], x), - longitude=(["points"], y), - levelist=(["points"], z), - ), + coords=coords_dict, ) for mars_metadata in self.mars_metadata[0]: ds.attrs[mars_metadata] = self.mars_metadata[0][mars_metadata] diff --git a/covjsonkit/decoder/Frame.py b/covjsonkit/decoder/Frame.py index 67a15b6..41348f6 100644 --- a/covjsonkit/decoder/Frame.py +++ b/covjsonkit/decoder/Frame.py @@ -8,6 +8,22 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering: derive x/y/z positions from + # the ``composite`` axis's ``coordinates`` labels (reads both x/y[/z] + # and legacy latitude/longitude/levelist labels). + labels = self.domains[0]["axes"]["composite"]["coordinates"] + self.x_idx = self._label_index(labels, ("x", "longitude")) + self.y_idx = self._label_index(labels, ("y", "latitude")) + self.z_idx = self._label_index(labels, ("z", "levelist"), required=False) + + @staticmethod + def _label_index(labels, names, required=True): + for name in names: + if name in labels: + return labels.index(name) + if required: + raise ValueError(f"None of {names} found in composite coordinates {labels}") + return None def get_domains(self): domains = [] @@ -57,12 +73,15 @@ def to_geojson(self): param_vals["datetime"] = datetime if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata + geom_coords = [lonlat[self.x_idx], lonlat[self.y_idx]] + if self.z_idx is not None: + geom_coords.append(lonlat[self.z_idx]) features.append( { "type": "Feature", "geometry": { "type": "Point", - "coordinates": [lonlat[1], lonlat[0], lonlat[2]], + "coordinates": geom_coords, }, "properties": param_vals, } @@ -91,8 +110,8 @@ def to_xarray(self): times.append(domain["axes"]["t"]["values"][0]) if i == 0: for coord in domain["axes"]["composite"]["values"]: - x.append(float(coord[0])) - y.append(float(coord[1])) + x.append(float(coord[self.x_idx])) + y.append(float(coord[self.y_idx])) n_points = len(x) diff --git a/covjsonkit/decoder/Shapefile.py b/covjsonkit/decoder/Shapefile.py index 3342e04..9b82b0b 100644 --- a/covjsonkit/decoder/Shapefile.py +++ b/covjsonkit/decoder/Shapefile.py @@ -8,6 +8,22 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering: derive x/y/z positions from + # the ``composite`` axis's ``coordinates`` labels (reads both x/y[/z] + # and legacy latitude/longitude/levelist labels). + labels = self.domains[0]["axes"]["composite"]["coordinates"] + self.x_idx = self._label_index(labels, ("x", "longitude")) + self.y_idx = self._label_index(labels, ("y", "latitude")) + self.z_idx = self._label_index(labels, ("z", "levelist"), required=False) + + @staticmethod + def _label_index(labels, names, required=True): + for name in names: + if name in labels: + return labels.index(name) + if required: + raise ValueError(f"None of {names} found in composite coordinates {labels}") + return None def get_domains(self): domains = [] @@ -60,12 +76,15 @@ def to_geojson(self): param_vals["datetime"] = datetime if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata + geom_coords = [lonlat[self.x_idx], lonlat[self.y_idx]] + if self.z_idx is not None: + geom_coords.append(lonlat[self.z_idx]) features.append( { "type": "Feature", "geometry": { "type": "Point", - "coordinates": [lonlat[1], lonlat[0], lonlat[2]], + "coordinates": geom_coords, }, "properties": param_vals, } @@ -94,8 +113,8 @@ def to_xarray(self): times.append(domain["axes"]["t"]["values"][0]) if i == 0: for coord in domain["axes"]["composite"]["values"]: - x.append(float(coord[0])) - y.append(float(coord[1])) + x.append(float(coord[self.x_idx])) + y.append(float(coord[self.y_idx])) n_points = len(x) diff --git a/covjsonkit/decoder/Wkt.py b/covjsonkit/decoder/Wkt.py index 8f16707..2cf512a 100644 --- a/covjsonkit/decoder/Wkt.py +++ b/covjsonkit/decoder/Wkt.py @@ -8,6 +8,22 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering: derive x/y/z positions from + # the ``composite`` axis's ``coordinates`` labels (reads both x/y[/z] + # and legacy latitude/longitude/levelist labels). + labels = self.domains[0]["axes"]["composite"]["coordinates"] + self.x_idx = self._label_index(labels, ("x", "longitude")) + self.y_idx = self._label_index(labels, ("y", "latitude")) + self.z_idx = self._label_index(labels, ("z", "levelist"), required=False) + + @staticmethod + def _label_index(labels, names, required=True): + for name in names: + if name in labels: + return labels.index(name) + if required: + raise ValueError(f"None of {names} found in composite coordinates {labels}") + return None def get_domains(self): domains = [] @@ -60,10 +76,13 @@ def to_geojson(self): param_vals["datetime"] = datetime if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata + geom_coords = [lonlat[self.x_idx], lonlat[self.y_idx]] + if self.z_idx is not None: + geom_coords.append(lonlat[self.z_idx]) features.append( { "type": "Feature", - "geometry": {"type": "Point", "coordinates": [lonlat[1], lonlat[0], lonlat[2]]}, + "geometry": {"type": "Point", "coordinates": geom_coords}, "properties": param_vals, } ) @@ -91,8 +110,8 @@ def to_xarray(self): times.append(domain["axes"]["t"]["values"][0]) if i == 0: for coord in domain["axes"]["composite"]["values"]: - x.append(float(coord[0])) - y.append(float(coord[1])) + x.append(float(coord[self.x_idx])) + y.append(float(coord[self.y_idx])) n_points = len(x) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index d80d96f..7fc169b 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -1,8 +1,6 @@ import logging import time -import pandas as pd - from .encoder import Encoder, normalize_step_value @@ -12,30 +10,30 @@ 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, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(json.dumps(new_coverage)) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): + # OGC-compliant MultiPoint domain: a single ``composite`` axis whose + # tuples are ordered [x, y] (== [longitude, latitude]) at the surface, + # or [x, y, z] (== [longitude, latitude, level]) when a vertical axis + # is present. ``t`` carries a single time value per coverage. coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0][ - "coordinates" - ] # self.pydantic_coverage.referencing[0].coordinates - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["t"] = {"values": coords["t"]} + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["x", "y", "z"] if include_z else ["x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -44,62 +42,70 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + # MultiPoint referencing is split into GeographicCRS (x, y), an optional + # VerticalCRS (z) and a TemporalRS (t) so each coordinate carries the + # correct coordinate reference system. + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a MultiPoint CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + # Add reference system (split GeographicCRS / TemporalRS / VerticalCRS) + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates: composite tuples ordered [x, y(, z)]. coords = { "composite": [], - "dataType": "tuple", "t": [str(x) for x in dataset["datetimes"].values], } for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([lon, lat, level]) + else: + coords["composite"].append([lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -112,7 +118,7 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) # Return the generated CoverageJSON return self.covjson @@ -129,20 +135,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = [0] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -156,15 +156,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) - # for k, v in range_dict.items(): - # if k == key: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = range_dict[key] else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += range_dict[key] levels = fields["levels"] @@ -178,12 +175,18 @@ def from_polytope(self, result, date_key: str = "date") -> dict: logging.debug("The fields retrieved were: %s", fields) # noqa: E501 logging.debug("The range_dict created was: %s", range_dict) # noqa: E501 + # Build composite tuples ordered [x, y] (== [lon, lat]) at the surface, + # or [x, y, z] (== [lon, lat, level]) when a vertical axis is present. + # walk_tree stores each spatial point as [lat, lon]. for date in coords.keys(): coord = coords[date]["composite"] coords[date]["composite"] = [] for level in levels: for cor in coord: - coords[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coords[date]["composite"].append([cor[1], cor[0], level]) + else: + coords[date]["composite"].append([cor[1], cor[0]]) for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -198,7 +201,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], include_z) return self.covjson @@ -214,6 +217,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) @@ -225,15 +229,8 @@ def from_polytope_month(self, result): start = time.time() logging.debug("Coords creation: %s", start) - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] if fields["param"] == 0: @@ -244,7 +241,7 @@ def from_polytope_month(self, result): logging.debug("The parameters added were: %s", self.parameters) # Build per-date composite coordinates, expanding each spatial point - # across all levels to produce [lat, lon, level] tuples. + # across all levels when a vertical axis is present. coordinates = {} for date in fields["dates"]: t_val = f"{date}-01T00:00:00Z" @@ -254,7 +251,10 @@ def from_polytope_month(self, result): coord = coords.get(date, {}).get("composite", []) for level in levels: for cor in coord: - coordinates[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[date]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[date]["composite"].append([cor[1], cor[0]]) end = time.time() logging.debug("Coords creation: %s", end) @@ -275,102 +275,10 @@ def from_polytope_month(self, result): mm = mars_metadata.copy() mm["number"] = num mm["Forecast date"] = f"{date}-01T00:00:00Z" - self.add_coverage(mm, coordinates[date], val_dict) + self.add_coverage(mm, coordinates[date], val_dict, include_z) end = time.time() logging.debug("Coverage creation: %s", end) logging.debug("Coverage creation takes: %s", end - start) return self.covjson - coords = {} - mars_metadata = {} - range_dict = {} - fields = {} - fields["lat"] = 0 - fields["param"] = 0 - fields["number"] = [0] - fields["step"] = [0] - fields["dates"] = [] - fields["levels"] = [0] - fields["times"] = [] - - start = time.time() - logging.debug("Tree walking starts at: %s", start) # noqa: E501 - self.walk_tree_step(result, fields, coords, mars_metadata, range_dict) - end = time.time() - delta = end - start - logging.debug("Tree walking ends at: %s", end) # noqa: E501 - logging.debug("Tree walking takes: %s", delta) # noqa: E501 - - start = time.time() - logging.debug("Coords creation: %s", start) # noqa: E501 - - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) - - coordinates = {} - - levels = fields["levels"] - if fields["param"] == 0: - raise ValueError("No data was returned.") - for para in fields["param"]: - self.add_parameter(para) - - logging.debug("The parameters added were: %s", self.parameters) # noqa: E501 - - # points = len(coords[fields["dates"][0]]["composite"]) - - coordinates = {} - - for date in coords.keys(): - for dt in coords[date]["t"]: - coordinates[dt] = {} - coordinates[dt]["composite"] = [] - coordinates[dt]["t"] = [dt] - coord = coords[date]["composite"] - for level in levels: - for cor in coord: - coordinates[dt]["composite"].append([cor[0], cor[1], level]) - - end = time.time() - delta = end - start - logging.debug("Coords creation: %s", end) # noqa: E501 - logging.debug("Coords creation: %s", delta) # noqa: E501 - - start = time.time() - logging.debug("Coverage creation: %s", start) # noqa: E501 - - val_dict = {} - for i, t in enumerate(fields["times"]): - val_dict[t] = {} - for level in fields["levels"]: - for num in fields["number"]: - for date in fields["dates"]: - for para in fields["param"]: - val_dict[t][para] = [] - key = (date, level, num, para) - vals = int(len(range_dict[key]) / len(fields["times"])) - # for val in range_dict[key]: - # vals.append(val[i]) - val_dict[t][para].extend(range_dict[key][i * vals : (i + 1) * vals]) - for para in fields["param"]: - val_dict[t][para] = [item for sublist in val_dict[t][para] for item in sublist] - mm = mars_metadata.copy() - mm["number"] = num - # mm["Forecast date"] = date - datetime = pd.Timestamp(date) + t - self.add_coverage(mm, coordinates[str(datetime).split("+")[0] + "Z"], val_dict[t]) - - end = time.time() - delta = end - start - logging.debug("Coverage creation: %s", end) # noqa: E501 - logging.debug("Coverage creation: %s", delta) # noqa: E501 - - return self.covjson diff --git a/covjsonkit/encoder/Circle.py b/covjsonkit/encoder/Circle.py index 7c2beb9..67c7700 100644 --- a/covjsonkit/encoder/Circle.py +++ b/covjsonkit/encoder/Circle.py @@ -9,30 +9,26 @@ 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, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(json.dumps(new_coverage)) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0][ - "coordinates" - ] # self.pydantic_coverage.referencing[0].coordinates - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["t"] = {"values": coords["t"]} + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["x", "y", "z"] if include_z else ["x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -41,62 +37,66 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a MultiPoint CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates: composite tuples ordered [x, y(, z)]. coords = { "composite": [], - "dataType": "tuple", "t": [str(x) for x in dataset["datetimes"].values], } for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([lon, lat, level]) + else: + coords["composite"].append([lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -109,9 +109,8 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) - # Return the generated CoverageJSON return self.covjson def from_polytope(self, result, date_key: str = "date") -> dict: @@ -126,21 +125,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = [0] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -154,15 +147,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) - # for k, v in range_dict.items(): - # if k == key: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = range_dict[key] else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += range_dict[key] levels = fields["levels"] @@ -181,7 +171,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: coords[date]["composite"] = [] for level in levels: for cor in coord: - coords[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coords[date]["composite"].append([cor[1], cor[0], level]) + else: + coords[date]["composite"].append([cor[1], cor[0]]) for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -196,7 +189,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], include_z) return self.covjson @@ -212,6 +205,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False logging.debug("Tree walking starts") self.walk_tree_month(result, fields, coords, mars_metadata, range_dict) @@ -220,15 +214,8 @@ def from_polytope_month(self, result): logging.debug("The values returned from walking tree: %s", range_dict) logging.debug("The coordinates returned from walking tree: %s", coords) - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] if fields["param"] == 0: @@ -241,7 +228,7 @@ def from_polytope_month(self, result): logging.debug("The range_dict created was: %s", range_dict) # Build per-date composite coordinates, expanding each spatial point - # across all levels to produce [lat, lon, level] tuples. + # across all levels when a vertical axis is present. coordinates = {} for date in fields["dates"]: t_val = f"{date}-01T00:00:00Z" @@ -251,7 +238,10 @@ def from_polytope_month(self, result): coord = coords.get(date, {}).get("composite", []) for level in levels: for cor in coord: - coordinates[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[date]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[date]["composite"].append([cor[1], cor[0]]) for num in fields["number"]: for date in fields["dates"]: @@ -265,6 +255,6 @@ def from_polytope_month(self, result): mm = mars_metadata.copy() mm["number"] = num mm["Forecast date"] = f"{date}-01T00:00:00Z" - self.add_coverage(mm, coordinates[date], val_dict) + self.add_coverage(mm, coordinates[date], val_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/Frame.py b/covjsonkit/encoder/Frame.py index 22bf4c6..309d41f 100644 --- a/covjsonkit/encoder/Frame.py +++ b/covjsonkit/encoder/Frame.py @@ -9,30 +9,26 @@ 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, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(cov) self.covjson["coverages"].append(new_coverage) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0][ - "coordinates" - ] # self.pydantic_coverage.referencing[0].coordinates - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["t"] = {"values": coords["t"]} + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["x", "y", "z"] if include_z else ["x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -41,62 +37,66 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a MultiPoint CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates: composite tuples ordered [x, y(, z)]. coords = { "composite": [], - "dataType": "tuple", "t": [str(x) for x in dataset["datetimes"].values], } for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([lon, lat, level]) + else: + coords["composite"].append([lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -109,9 +109,8 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) - # Return the generated CoverageJSON return self.covjson def from_polytope(self, result, date_key: str = "date") -> dict: @@ -126,21 +125,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = 0 fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -154,7 +147,6 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) for k, v in range_dict.items(): @@ -162,7 +154,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = v else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += v levels = fields["levels"] @@ -176,7 +168,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: coords[date]["composite"] = [] for level in levels: for cor in coord: - coords[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coords[date]["composite"].append([cor[1], cor[0], level]) + else: + coords[date]["composite"].append([cor[1], cor[0]]) for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -191,7 +186,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], include_z) return self.covjson @@ -207,6 +202,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False logging.debug("Tree walking starts") self.walk_tree_month(result, fields, coords, mars_metadata, range_dict) @@ -215,15 +211,8 @@ def from_polytope_month(self, result): logging.debug("The values returned from walking tree: %s", range_dict) logging.debug("The coordinates returned from walking tree: %s", coords) - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] for para in fields["param"]: @@ -241,7 +230,10 @@ def from_polytope_month(self, result): coord = coords.get(date, {}).get("composite", []) for level in levels: for cor in coord: - coordinates[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[date]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[date]["composite"].append([cor[1], cor[0]]) for num in fields["number"]: val_dict = {} @@ -255,6 +247,11 @@ def from_polytope_month(self, result): val_dict[para].extend([item for sublist in vals for item in sublist]) mm = mars_metadata.copy() mm["number"] = num - self.add_coverage(mm, coordinates[fields["dates"][0]] if fields["dates"] else {}, val_dict) + self.add_coverage( + mm, + coordinates[fields["dates"][0]] if fields["dates"] else {"composite": [], "t": []}, + val_dict, + include_z, + ) return self.covjson diff --git a/covjsonkit/encoder/Shapefile.py b/covjsonkit/encoder/Shapefile.py index 897b7fb..26b0a50 100644 --- a/covjsonkit/encoder/Shapefile.py +++ b/covjsonkit/encoder/Shapefile.py @@ -9,30 +9,26 @@ 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, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(cov) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0][ - "coordinates" - ] # self.pydantic_coverage.referencing[0].coordinates - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["t"] = {"values": coords["t"]} + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["x", "y", "z"] if include_z else ["x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -41,62 +37,66 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a MultiPoint CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates: composite tuples ordered [x, y(, z)]. coords = { "composite": [], - "dataType": "tuple", "t": [str(x) for x in dataset["datetimes"].values], } for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([lon, lat, level]) + else: + coords["composite"].append([lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -109,9 +109,8 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) - # Return the generated CoverageJSON return self.covjson def from_polytope(self, result, date_key: str = "date") -> dict: @@ -126,21 +125,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = 0 fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -154,7 +147,6 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) for k, v in range_dict.items(): @@ -162,7 +154,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = v else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += v levels = fields["levels"] @@ -176,7 +168,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: coords[date]["composite"] = [] for level in levels: for cor in coord: - coords[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coords[date]["composite"].append([cor[1], cor[0], level]) + else: + coords[date]["composite"].append([cor[1], cor[0]]) for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -191,7 +186,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], include_z) return self.covjson @@ -207,21 +202,15 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree_month(result, fields, coords, mars_metadata, range_dict) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] if fields["param"] == 0: @@ -235,7 +224,6 @@ def from_polytope_month(self, result): raise ValueError("No year/month data was found in the result tree.") # Build composite coordinates: expand spatial points across all levels per date. - # All dates are merged into a single coverage via the t axis. coordinates = {} for date in fields["dates"]: t_val = f"{date}-01T00:00:00Z" @@ -245,7 +233,10 @@ def from_polytope_month(self, result): coord = coords.get(date, {}).get("composite", []) for level in levels: for cor in coord: - coordinates[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[date]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[date]["composite"].append([cor[1], cor[0]]) for num in fields["number"]: val_dict = {} @@ -259,6 +250,11 @@ def from_polytope_month(self, result): val_dict[para].extend([item for sublist in vals for item in sublist]) mm = mars_metadata.copy() mm["number"] = num - self.add_coverage(mm, coordinates[fields["dates"][0]] if fields["dates"] else {}, val_dict) + self.add_coverage( + mm, + coordinates[fields["dates"][0]] if fields["dates"] else {"composite": [], "t": []}, + val_dict, + include_z, + ) return self.covjson diff --git a/covjsonkit/encoder/Wkt.py b/covjsonkit/encoder/Wkt.py index 931e6e3..86f84c3 100644 --- a/covjsonkit/encoder/Wkt.py +++ b/covjsonkit/encoder/Wkt.py @@ -13,30 +13,26 @@ 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, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(cov) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0][ - "coordinates" - ] # self.pydantic_coverage.referencing[0].coordinates - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["t"] = {"values": coords["t"]} + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["x", "y", "z"] if include_z else ["x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -45,62 +41,66 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a MultiPoint CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates: composite tuples ordered [x, y(, z)]. coords = { "composite": [], - "dataType": "tuple", "t": [str(x) for x in dataset["datetimes"].values], } for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([lon, lat, level]) + else: + coords["composite"].append([lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -113,9 +113,8 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) - # Return the generated CoverageJSON return self.covjson def from_polytope(self, result, date_key: str = "date") -> dict: @@ -130,21 +129,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = [0] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -158,15 +151,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) - # for k, v in range_dict.items(): - # if k == key: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = range_dict[key] else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += range_dict[key] levels = fields["levels"] @@ -185,7 +175,10 @@ def from_polytope(self, result, date_key: str = "date") -> dict: coords[date]["composite"] = [] for level in levels: for cor in coord: - coords[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coords[date]["composite"].append([cor[1], cor[0], level]) + else: + coords[date]["composite"].append([cor[1], cor[0]]) for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -200,7 +193,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], include_z) return self.covjson @@ -216,6 +209,7 @@ def from_polytope_step(self, result): fields["dates"] = [] fields["levels"] = [0] fields["times"] = [] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 @@ -228,17 +222,8 @@ def from_polytope_step(self, result): start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) - - coordinates = {} + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] if fields["param"] == 0: @@ -248,8 +233,6 @@ def from_polytope_step(self, result): logging.debug("The parameters added were: %s", self.parameters) # noqa: E501 - # points = len(coords[fields["dates"][0]]["composite"]) - coordinates = {} for date in coords.keys(): @@ -260,7 +243,10 @@ def from_polytope_step(self, result): coord = coords[date]["composite"] for level in levels: for cor in coord: - coordinates[dt]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[dt]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[dt]["composite"].append([cor[1], cor[0]]) end = time.time() delta = end - start @@ -286,7 +272,7 @@ def from_polytope_step(self, result): mm["number"] = num # mm["Forecast date"] = date datetime = pd.Timestamp(date) + t - self.add_coverage(mm, coordinates[str(datetime).split("+")[0] + "Z"], val_dict) + self.add_coverage(mm, coordinates[str(datetime).split("+")[0] + "Z"], val_dict, include_z) end = time.time() delta = end - start @@ -310,6 +296,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) # noqa: E501 @@ -321,15 +308,8 @@ def from_polytope_month(self, result): start = time.time() logging.debug("Coords creation: %s", start) # noqa: E501 - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) levels = fields["levels"] if fields["param"] == 0: @@ -352,7 +332,10 @@ def from_polytope_month(self, result): coord = coords.get(date, {}).get("composite", []) for level in levels: for cor in coord: - coordinates[date]["composite"].append([cor[0], cor[1], level]) + if include_z: + coordinates[date]["composite"].append([cor[1], cor[0], level]) + else: + coordinates[date]["composite"].append([cor[1], cor[0]]) end = time.time() logging.debug("Coords creation: %s", end) # noqa: E501 @@ -373,7 +356,7 @@ def from_polytope_month(self, result): mm = mars_metadata.copy() mm["number"] = num mm["Forecast date"] = f"{date}-01T00:00:00Z" - self.add_coverage(mm, coordinates[date], val_dict) + self.add_coverage(mm, coordinates[date], val_dict, include_z) end = time.time() logging.debug("Coverage creation: %s", end) # noqa: E501 diff --git a/tests/conftest.py b/tests/conftest.py index fbbd289..1d8b4cb 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -54,10 +54,10 @@ def assert_valid_covjson(covjson): "number": 0, } -COMPOSITE_TWO_POINTS_XYZ = { +COMPOSITE_TWO_POINTS = { "dataType": "tuple", - "coordinates": ["x", "y", "z"], - "values": [[48.0, 11.0, 0], [50.0, 12.0, 0]], + "coordinates": ["x", "y"], + "values": [[11.0, 48.0], [12.0, 50.0]], } diff --git a/tests/test_encoder_bounding_box_from_polytope.py b/tests/test_encoder_bounding_box_from_polytope.py index ec21c75..1df0548 100644 --- a/tests/test_encoder_bounding_box_from_polytope.py +++ b/tests/test_encoder_bounding_box_from_polytope.py @@ -12,12 +12,13 @@ from covjsonkit.api import Covjsonkit -# BoundingBox uses lat/lon/levelist coordinates (not x/y/z) and its -# reforecast metadata intentionally omits "step" (step varies per coverage). +# BoundingBox emits spec-compliant composite tuples ordered [x, y] (== [lon, lat]). +# Surface (sfc) data carries no vertical axis, so no z component. Reforecast +# metadata intentionally omits "step" (step varies per coverage). COMPOSITE_TWO_POINTS = { "dataType": "tuple", - "coordinates": ["latitude", "longitude", "levelist"], - "values": [[48.0, 11.0, 0], [50.0, 12.0, 0]], + "coordinates": ["x", "y"], + "values": [[11.0, 48.0], [12.0, 50.0]], } EXPECTED_REFORECAST_METADATA = { @@ -52,7 +53,7 @@ def test_single_date_single_step_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -144,7 +145,7 @@ def test_reforecast_single_hdate_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } diff --git a/tests/test_encoder_circle_from_polytope.py b/tests/test_encoder_circle_from_polytope.py index 1baceb0..4e94d5c 100644 --- a/tests/test_encoder_circle_from_polytope.py +++ b/tests/test_encoder_circle_from_polytope.py @@ -10,11 +10,11 @@ THREE_POINTS_COMPOSITE = { "dataType": "tuple", - "coordinates": ["latitude", "longitude", "levelist"], + "coordinates": ["x", "y"], "values": [ - [48.0, 11.0, 0], - [49.0, 11.5, 0], - [50.0, 12.0, 0], + [11.0, 48.0], + [11.5, 49.0], + [12.0, 50.0], ], } @@ -39,7 +39,7 @@ def test_single_date_single_step_three_points(self): "type": "NdArray", "dataType": "float", "shape": [3], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.5, 266.1], } } @@ -80,7 +80,7 @@ def test_reforecast_single_hdate_three_points(self): "type": "NdArray", "dataType": "float", "shape": [3], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.5, 266.1], } } diff --git a/tests/test_encoder_frame_from_polytope.py b/tests/test_encoder_frame_from_polytope.py index 0330dc4..743ab3e 100644 --- a/tests/test_encoder_frame_from_polytope.py +++ b/tests/test_encoder_frame_from_polytope.py @@ -1,6 +1,6 @@ import numpy as np from conftest import ( - COMPOSITE_TWO_POINTS_XYZ, + COMPOSITE_TWO_POINTS, REFORECAST_METADATA_BASE, forecast_tree, reforecast_branch, @@ -20,7 +20,7 @@ def test_single_date_single_step_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-01-01T00:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -28,7 +28,7 @@ def test_single_date_single_step_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -62,7 +62,7 @@ def test_reforecast_single_hdate_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-07-14T06:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -70,7 +70,7 @@ def test_reforecast_single_hdate_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -99,7 +99,7 @@ def test_reforecast_two_hdates_two_points(self): for cov, (fc_date, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { "t": {"values": [fc_date]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == vals assert cov["mars:metadata"] == { diff --git a/tests/test_encoder_shapefile_from_polytope.py b/tests/test_encoder_shapefile_from_polytope.py index 611d7f1..24289e0 100644 --- a/tests/test_encoder_shapefile_from_polytope.py +++ b/tests/test_encoder_shapefile_from_polytope.py @@ -1,6 +1,6 @@ import numpy as np from conftest import ( - COMPOSITE_TWO_POINTS_XYZ, + COMPOSITE_TWO_POINTS, REFORECAST_METADATA_BASE, forecast_tree, reforecast_branch, @@ -20,7 +20,7 @@ def test_single_date_single_step_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-01-01T00:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -28,7 +28,7 @@ def test_single_date_single_step_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -62,7 +62,7 @@ def test_reforecast_single_hdate_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-07-14T06:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -70,7 +70,7 @@ def test_reforecast_single_hdate_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -96,7 +96,7 @@ def test_reforecast_two_hdates_two_points(self): for cov, (t_vals, range_vals, fc_date) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { "t": {"values": t_vals}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == range_vals assert cov["mars:metadata"] == {"Forecast date": fc_date, **REFORECAST_METADATA_BASE} diff --git a/tests/test_encoder_wkt_from_polytope.py b/tests/test_encoder_wkt_from_polytope.py index 0862f3f..1c4721e 100644 --- a/tests/test_encoder_wkt_from_polytope.py +++ b/tests/test_encoder_wkt_from_polytope.py @@ -1,6 +1,6 @@ import numpy as np from conftest import ( - COMPOSITE_TWO_POINTS_XYZ, + COMPOSITE_TWO_POINTS, REFORECAST_METADATA_BASE, forecast_tree, reforecast_branch, @@ -20,7 +20,7 @@ def test_single_date_single_step_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-01-01T00:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -28,7 +28,7 @@ def test_single_date_single_step_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -62,7 +62,7 @@ def test_reforecast_single_hdate_two_points(self): assert cov["domain"]["axes"] == { "t": {"values": ["2025-07-14T06:00:00Z"]}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"] == { @@ -70,7 +70,7 @@ def test_reforecast_single_hdate_two_points(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -99,7 +99,7 @@ def test_reforecast_two_hdates_two_points(self): for cov, (t_vals, range_vals, fc_date) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { "t": {"values": t_vals}, - "composite": COMPOSITE_TWO_POINTS_XYZ, + "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == range_vals assert cov["mars:metadata"] == {"Forecast date": fc_date, **REFORECAST_METADATA_BASE} diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index bab0373..9d45784 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -279,3 +279,141 @@ def test_vertical_profile_legacy_and_new_xarray_equivalent(self): assert "latitude" in legacy_ds.coords assert "longitude" in legacy_ds.coords assert "levelist" in legacy_ds.coords + + +def _to_legacy_composite(covjson): + """Rewrite a spec-compliant MultiPoint (composite-axis) CoverageCollection + into the pre-#129 legacy shape: composite ``coordinates`` labelled + latitude/longitude[/levelist] with tuple values reordered to + [lat, lon[, level]], plus a single combined GeographicCRS referencing block. + Used to exercise composite decoder backwards-compat. + """ + legacy = copy.deepcopy(covjson) + + first_labels = legacy["coverages"][0]["domain"]["axes"]["composite"]["coordinates"] + + def _idx(labels, names): + for name in names: + if name in labels: + return labels.index(name) + return None + + include_z = _idx(first_labels, ("z", "levelist")) is not None + + for cov in legacy["coverages"]: + composite = cov["domain"]["axes"]["composite"] + labels = composite["coordinates"] + x_idx = _idx(labels, ("x", "longitude")) + y_idx = _idx(labels, ("y", "latitude")) + z_idx = _idx(labels, ("z", "levelist")) + new_values = [] + for tup in composite["values"]: + entry = [tup[y_idx], tup[x_idx]] + if z_idx is not None: + entry.append(tup[z_idx]) + new_values.append(entry) + composite["values"] = new_values + composite["coordinates"] = ["latitude", "longitude", "levelist"] if include_z else ["latitude", "longitude"] + + coords = ["latitude", "longitude"] + if include_z: + coords.append("levelist") + legacy["referencing"] = [ + { + "coordinates": coords, + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ] + return legacy + + +class TestBoundingBoxSpecCompliance: + """BoundingBox (MultiPoint) encoder output validates and decodes with back-compat.""" + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_bounding_box_output_validates(self, tree_factory): + covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(tree_factory()) + assert_valid_covjson(covjson) + + def test_bounding_box_month_validates(self): + covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope_month(_month_tree()) + assert_valid_covjson(covjson) + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_bounding_box_legacy_and_new_geojson_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(tree_factory()) + legacy_covjson = _to_legacy_composite(new_covjson) + + new_gj = Covjsonkit().decode(new_covjson).to_geojson() + legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() + + assert legacy_gj == new_gj + + def test_bounding_box_geojson_lonlat_order(self): + # Surface data: geometry must be [lon, lat] in the correct order. + covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(_surface_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0] + + def test_bounding_box_geojson_lonlat_z_order(self): + # Level data: geometry must be [lon, lat, level] in the correct order. + covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(_level_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0, 500] + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_bounding_box_legacy_and_new_xarray_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(tree_factory()) + legacy_covjson = _to_legacy_composite(new_covjson) + + new_ds = Covjsonkit().decode(new_covjson).to_xarray() + legacy_ds = Covjsonkit().decode(legacy_covjson).to_xarray() + + assert new_ds.identical(legacy_ds) + + # Decoder read the composite by coordinate LABEL, not position: the + # resulting dataset carries latitude/longitude, plus levelist when the + # source had a z component. + assert "latitude" in legacy_ds.coords + assert "longitude" in legacy_ds.coords + has_z = "z" in new_covjson["coverages"][0]["domain"]["axes"]["composite"]["coordinates"] + assert ("levelist" in legacy_ds.coords) == has_z + + +class TestCompositeFeatureSpecCompliance: + """Circle/Frame/Shapefile/Wkt (MultiPoint composite-axis) encoders must emit + spec-compliant output and decode with composite backwards-compat.""" + + FEATURES = ["Circle", "Frame", "Shapefile", "Polygon"] + + @pytest.mark.parametrize("feature", FEATURES) + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_output_validates(self, feature, tree_factory): + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope(tree_factory()) + assert_valid_covjson(covjson) + + @pytest.mark.parametrize("feature", FEATURES) + def test_geojson_lonlat_order(self, feature): + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope(_surface_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0] + + @pytest.mark.parametrize("feature", FEATURES) + def test_geojson_lonlat_z_order(self, feature): + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope(_level_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0, 500] + + @pytest.mark.parametrize("feature", FEATURES) + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_legacy_and_new_geojson_equivalent(self, feature, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope(tree_factory()) + legacy_covjson = _to_legacy_composite(new_covjson) + + new_gj = Covjsonkit().decode(new_covjson).to_geojson() + legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() + + assert legacy_gj == new_gj diff --git a/tests/test_xarray_multipoint.py b/tests/test_xarray_multipoint.py index 4cf8cee..47055bd 100644 --- a/tests/test_xarray_multipoint.py +++ b/tests/test_xarray_multipoint.py @@ -32,10 +32,15 @@ def test_from_xarray(self): assert covjson_result["type"] == self.test_covjson["type"] assert len(covjson_result["coverages"]) == len(self.test_covjson["coverages"]) - assert ( - covjson_result["coverages"][0]["domain"]["axes"]["composite"]["values"] - == self.test_covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] - ) + # Encoder now emits spec-compliant composite: coordinates ["x","y","z"] with + # values ordered [lon, lat, level]. The legacy fixture uses + # ["latitude","longitude","levelist"] with [lat, lon, level]; the decoder + # reads it back-compatibly, so the round-trip swaps lat/lon into spec order. + composite = covjson_result["coverages"][0]["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["x", "y", "z"] + legacy_values = self.test_covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] + expected_values = [[lat_lon_lvl[1], lat_lon_lvl[0], lat_lon_lvl[2]] for lat_lon_lvl in legacy_values] + assert composite["values"] == expected_values assert ( covjson_result["coverages"][0]["mars:metadata"]["number"] == self.test_covjson["coverages"][0]["mars:metadata"]["number"] From ac810c8e4e78e2ef9a29d6c572f35e4effd6ebed Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 4 Sep 2026 14:14:05 +0200 Subject: [PATCH 03/10] Update path to comply with ogc spec --- covjsonkit/decoder/Path.py | 69 +++++++-- covjsonkit/encoder/Path.py | 181 ++++++++++++----------- tests/test_encoder_path_from_polytope.py | 22 +-- tests/test_spec_compliance.py | 78 ++++++++++ tests/test_xarray_path.py | 13 +- 5 files changed, 249 insertions(+), 114 deletions(-) diff --git a/covjsonkit/decoder/Path.py b/covjsonkit/decoder/Path.py index 5001a24..6365aa7 100644 --- a/covjsonkit/decoder/Path.py +++ b/covjsonkit/decoder/Path.py @@ -9,6 +9,45 @@ def __init__(self, covjson): super().__init__(covjson) self.domains = self.get_domains() self.ranges = self.get_ranges() + # Backwards-compatible composite ordering for the Trajectory ``composite`` + # tuple. Both the spec-compliant output and the legacy output label the + # spatial components ``x``/``y``, but with SWAPPED meaning: + # * spec: values ordered [t, x=lon, y=lat(, z)] with split referencing + # (a GeographicCRS whose coordinates are exactly ["x", "y"]). + # * legacy: values ordered [t, lat, lon(, z)] with a single combined + # GeographicCRS whose coordinates are ["t", "x", "y", "z"]. + # Labels alone cannot disambiguate these, so the referencing structure is + # used: a split GeographicCRS(["x", "y"]) marks spec-compliant output. + self.t_idx, self.x_idx, self.y_idx, self.z_idx = self._composite_indices() + + def _composite_indices(self): + labels = self.domains[0]["axes"]["composite"]["coordinates"] + refs = self.covjson.get("referencing", []) + is_spec = any( + ref.get("system", {}).get("type") == "GeographicCRS" and set(ref.get("coordinates", [])) == {"x", "y"} + for ref in refs + ) + + t_idx = labels.index("t") if "t" in labels else None + z_idx = None + for name in ("z", "levelist"): + if name in labels: + z_idx = labels.index(name) + break + + if is_spec: + # Spec order: label ``x`` holds longitude, label ``y`` holds latitude. + x_idx = labels.index("x") if "x" in labels else labels.index("longitude") + y_idx = labels.index("y") if "y" in labels else labels.index("latitude") + else: + # Legacy order: values are [t, lat, lon(, z)]; the mislabeled ``x`` + # component actually holds latitude and ``y`` holds longitude. + lat_label = "x" if "x" in labels else "latitude" + lon_label = "y" if "y" in labels else "longitude" + y_idx = labels.index(lat_label) # latitude + x_idx = labels.index(lon_label) # longitude + + return t_idx, x_idx, y_idx, z_idx def get_domains(self): domains = [] @@ -55,13 +94,17 @@ def to_geojson(self): param_vals[key] = values[key][idx] if "mars:metadata" in coverage: param_vals["mars:metadata"] = mars_metadata - param_vals["datetime"] = coord[0] + if self.t_idx is not None: + param_vals["datetime"] = coord[self.t_idx] + geom_coords = [coord[self.x_idx], coord[self.y_idx]] # lon, lat + if self.z_idx is not None: + geom_coords.append(coord[self.z_idx]) features.append( { "type": "Feature", "geometry": { "type": "Point", - "coordinates": [[coord[1], coord[2], coord[3]]], # lon, lat, z + "coordinates": geom_coords, }, "properties": param_vals, } @@ -76,15 +119,17 @@ def to_xarray(self): dataarraydict = {} # Get coordinates - x = [] - y = [] + longitude = [] + latitude = [] levelist = [] time = [] for coord in self.get_coordinates()["composite"]["values"]: - x.append(float(coord[1])) - y.append(float(coord[2])) - levelist.append(float(coord[3])) - time.append(coord[0]) + longitude.append(float(coord[self.x_idx])) + latitude.append(float(coord[self.y_idx])) + if self.z_idx is not None: + levelist.append(float(coord[self.z_idx])) + if self.t_idx is not None: + time.append(coord[self.t_idx]) values = {} for parameter in self.parameters: @@ -143,10 +188,10 @@ def to_xarray(self): datetimes=(["datetimes"], datetimes), number=(["number"], numbers), steps=(["steps"], steps), - points=(["points"], list(range(0, len(x)))), - latitude=(["points"], x), - longitude=(["points"], y), - levelist=(["points"], levelist), + points=(["points"], list(range(0, len(longitude)))), + latitude=(["points"], latitude), + longitude=(["points"], longitude), + **({"levelist": (["points"], levelist)} if self.z_idx is not None else {}), time=(["points"], time), ), ) diff --git a/covjsonkit/encoder/Path.py b/covjsonkit/encoder/Path.py index bf5bed6..fe91f3f 100644 --- a/covjsonkit/encoder/Path.py +++ b/covjsonkit/encoder/Path.py @@ -9,26 +9,30 @@ def __init__(self, type, domaintype): self.covjson["domainType"] = "Trajectory" self.covjson["coverages"] = [] - def add_coverage(self, mars_metadata, coords, values): + def add_coverage(self, mars_metadata, coords, values, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) + self.add_domain(new_coverage, coords, include_z) self.add_range(new_coverage, values) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(cov) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): + # OGC-compliant Trajectory domain: a single ``composite`` axis whose + # tuples are ordered [t, x, y] (== [time, longitude, latitude]) at the + # surface, or [t, x, y, z] (== [time, longitude, latitude, level]) when a + # vertical axis is present. Time is embedded in the composite tuple (a + # Trajectory has no separate ``t`` axis). coverage["domain"]["type"] = "Domain" coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["composite"] = {} - coverage["domain"]["axes"]["composite"]["dataType"] = "tuple" - coverage["domain"]["axes"]["composite"]["coordinates"] = self.covjson["referencing"][0]["coordinates"] - coverage["domain"]["axes"]["composite"]["values"] = coords["composite"] + coverage["domain"]["axes"]["composite"] = { + "dataType": "tuple", + "coordinates": ["t", "x", "y", "z"] if include_z else ["t", "x", "y"], + "values": coords["composite"], + } def add_range(self, coverage, values): for parameter in values.keys(): @@ -37,66 +41,68 @@ def add_range(self, coverage, values): coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = [str(param)] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = ["composite"] + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata + def _set_references(self, include_z): + # Trajectory referencing is split into GeographicCRS (x, y), a TemporalRS + # (t) and an optional VerticalCRS (z) so each coordinate carries the + # correct coordinate reference system. + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a MultiPoint CoverageJSON format. + Converts an xarray dataset into a Trajectory CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" self.covjson["domainType"] = "Trajectory" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - elif "z" in dataset.coords: - z_coord = "z" - else: - z_coord = "level" - - # Add reference system - self.add_reference( - { - "coordinates": ["t", x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + # Add reference system (split GeographicCRS / TemporalRS / VerticalCRS) + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates - coords = { - "composite": [], - "dataType": "tuple", - } + # Prepare coordinates: composite tuples ordered [t, x, y(, z)]. + coords = {"composite": []} for point in dataset["points"].values: - coords["composite"].append( - [ - float(dataset.isel(points=point).time.values), - float(dataset.isel(points=point).latitude.values), - float(dataset.isel(points=point).longitude.values), - float(dataset.isel(points=point).levelist.values), - ] - ) + t = float(dataset.isel(points=point).time.values) + lon = float(dataset.isel(points=point).longitude.values) + lat = float(dataset.isel(points=point).latitude.values) + if include_z: + level = float(dataset.isel(points=point).levelist.values) + coords["composite"].append([t, lon, lat, level]) + else: + coords["composite"].append([t, lon, lat]) for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -109,7 +115,7 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) # Return the generated CoverageJSON return self.covjson @@ -128,6 +134,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["levels"] = [0] fields["s"] = [] fields["l"] = [] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) @@ -141,16 +148,13 @@ def from_polytope(self, result, date_key: str = "date") -> dict: logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 logging.debug("The fields: %s", fields) - self.add_reference( - { - "coordinates": ["t", "x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) + # Build trajectory composite tuples ordered [t, x, y] (== [step, lon, + # lat]) at the surface, or [t, x, y, z] (== [step, lon, lat, level]) when + # a vertical axis is present. walk_tree stores each spatial point as + # [lat, lon]. for date in coords.keys(): coord = coords[date]["composite"] coords[date]["composite"] = [] @@ -160,17 +164,22 @@ def from_polytope(self, result, date_key: str = "date") -> dict: if (date, level, fields["number"][0], fields["param"][0], s) in range_dict: cor_len = len(range_dict[(date, level, fields["number"][0], fields["param"][0], s)]) end = start + cor_len - # Normalize the step value before it is placed into the composite (t) - # coordinate so that timedelta/Timedelta values are JSON serialisable. - # The original ``s`` is kept for the range_dict lookups above. + # Normalize the step value before it is placed into the + # composite (t) coordinate so that timedelta/Timedelta + # values are JSON serialisable. The original ``s`` is kept + # for the range_dict lookups above. s_norm = normalize_step_value(s) - if len(fields["levels"]) != 1: - for lev in fields["levels"]: + if include_z: + if len(fields["levels"]) != 1: + for lev in fields["levels"]: + for cor in coord[int(start) : int(end)]: + coords[date]["composite"].append([s_norm, cor[1], cor[0], lev]) + else: for cor in coord[int(start) : int(end)]: - coords[date]["composite"].append([s_norm, cor[0], cor[1], lev]) + coords[date]["composite"].append([s_norm, cor[1], cor[0], level]) else: for cor in coord[int(start) : int(end)]: - coords[date]["composite"].append([s_norm, cor[0], cor[1], level]) + coords[date]["composite"].append([s_norm, cor[1], cor[0]]) start = end logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 @@ -186,17 +195,13 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in set(fields["s"]): key = (date, level, num, para, s) - # for k, v in range_dict.items(): - # if k == key: if s not in combined_dict[date][num][para]: if key in range_dict: combined_dict[date][num][para][s] = range_dict[key] - # combined_dict[date][num][para][s] = range_dict[key] else: - # Cocatenate arrays + # Concatenate arrays if key in range_dict: combined_dict[date][num][para][s] += range_dict[key] @@ -225,7 +230,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm["Forecast date"] = date if "levelist" in mm: del mm["levelist"] - self.add_coverage(mm, coords[date], val_dict) + self.add_coverage(mm, coords[date], val_dict, include_z) return self.covjson @@ -242,6 +247,7 @@ def from_polytope_month(self, result): fields["dates"] = [] fields["levels"] = [0] fields["l"] = [] + fields["has_level_axis"] = False self.walk_tree_month(result, fields, coords, mars_metadata, range_dict) @@ -252,18 +258,11 @@ def from_polytope_month(self, result): logging.debug("The coordinates returned from walking tree: %s", coords) logging.debug("The fields: %s", fields) - self.add_reference( - { - "coordinates": ["t", "x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) - # Build trajectory composite coords: [month_key, lat, lon, level] - # Each (year, month) pair acts as the "t" element of the 4-tuple. + # Build trajectory composite tuples: [month_key, x, y(, z)]. + # Each (year, month) key acts as the "t" element of the tuple. for date in coords.keys(): coord = coords[date]["composite"] coords[date]["composite"] = [] @@ -272,13 +271,17 @@ def from_polytope_month(self, result): if (date, level, fields["number"][0], fields["param"][0]) in range_dict: cor_len = len(range_dict[(date, level, fields["number"][0], fields["param"][0])]) end_idx = start_idx + cor_len - if len(fields["levels"]) != 1: - for lev in fields["levels"]: + if include_z: + if len(fields["levels"]) != 1: + for lev in fields["levels"]: + for cor in coord[int(start_idx) : int(end_idx)]: + coords[date]["composite"].append([date, cor[1], cor[0], lev]) + else: for cor in coord[int(start_idx) : int(end_idx)]: - coords[date]["composite"].append([date, cor[0], cor[1], lev]) + coords[date]["composite"].append([date, cor[1], cor[0], level]) else: for cor in coord[int(start_idx) : int(end_idx)]: - coords[date]["composite"].append([date, cor[0], cor[1], level]) + coords[date]["composite"].append([date, cor[1], cor[0]]) start_idx = end_idx logging.debug("The coordinates returned from walking tree: %s", coords) @@ -307,6 +310,6 @@ def from_polytope_month(self, result): mm["number"] = num if "levelist" in mm: del mm["levelist"] - self.add_coverage(mm, coords[date], val_dict) + self.add_coverage(mm, coords[date], val_dict, include_z) return self.covjson diff --git a/tests/test_encoder_path_from_polytope.py b/tests/test_encoder_path_from_polytope.py index 7f5e389..b77de66 100644 --- a/tests/test_encoder_path_from_polytope.py +++ b/tests/test_encoder_path_from_polytope.py @@ -29,8 +29,8 @@ def test_two_points_single_step(self): assert cov["domain"]["axes"] == { "composite": { "dataType": "tuple", - "coordinates": ["t", "x", "y", "z"], - "values": [[0, 48.0, 11.0, 0], [0, 49.0, 12.0, 0]], + "coordinates": ["t", "x", "y"], + "values": [[0, 11.0, 48.0], [0, 12.0, 49.0]], } } @@ -39,7 +39,7 @@ def test_two_points_single_step(self): "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } @@ -56,15 +56,19 @@ def test_two_points_single_step(self): "number": 0, } - # Collection-level referencing + # Collection-level referencing (split GeographicCRS / TemporalRS) assert covjson["referencing"] == [ { - "coordinates": ["t", "x", "y", "z"], + "coordinates": ["x", "y"], "system": { "type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", }, - } + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, ] # Collection-level parameters @@ -79,8 +83,8 @@ class TestPathFromPolytopeReforecast: EXPECTED_AXES = { "composite": { "dataType": "tuple", - "coordinates": ["t", "x", "y", "z"], - "values": [[0, 48.0, 11.0, 0], [0, 50.0, 12.0, 0]], + "coordinates": ["t", "x", "y"], + "values": [[0, 11.0, 48.0], [0, 12.0, 50.0]], } } @@ -89,7 +93,7 @@ class TestPathFromPolytopeReforecast: "type": "NdArray", "dataType": "float", "shape": [2], - "axisNames": ["2t"], + "axisNames": ["composite"], "values": [264.9, 265.1], } } diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index 9d45784..5e3cc93 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -417,3 +417,81 @@ def test_legacy_and_new_geojson_equivalent(self, feature, tree_factory): legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() assert legacy_gj == new_gj + + +def _to_legacy_path(covjson): + """Convert a spec-compliant Trajectory (Path) covjson into the legacy + mislabeled form: a single combined GeographicCRS whose coordinates are + ["t", "x", "y", "z"], with composite tuples ordered [t, lat, lon, level]. + + Spec output orders composite tuples as [t, lon, lat] (surface) or + [t, lon, lat, level] (with a vertical axis); the legacy form always carries + a z component (0 at the surface). Used to assert the decoder reads both. + """ + legacy = copy.deepcopy(covjson) + legacy["referencing"] = [ + { + "coordinates": ["t", "x", "y", "z"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ] + for coverage in legacy["coverages"]: + composite = coverage["domain"]["axes"]["composite"] + labels = composite["coordinates"] + has_z = "z" in labels + new_values = [] + for tup in composite["values"]: + t, lon, lat = tup[0], tup[1], tup[2] + z = tup[3] if has_z else 0 + new_values.append([t, lat, lon, z]) + composite["coordinates"] = ["t", "x", "y", "z"] + composite["values"] = new_values + return legacy + + +class TestPathSpecCompliance: + """Path (Trajectory) encoder output validates and decodes with back-compat. + + The legacy Path form mislabels its spatial components (label ``x`` holds + latitude) and uses a single combined GeographicCRS; the spec form uses split + referencing with [t, x=lon, y=lat(, z)]. The decoder discriminates via the + referencing structure, so both must decode equivalently. + """ + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_path_output_validates(self, tree_factory): + covjson = Covjsonkit().encode("CoverageCollection", "Path").from_polytope(tree_factory()) + assert_valid_covjson(covjson) + + def test_path_month_validates(self): + covjson = Covjsonkit().encode("CoverageCollection", "Path").from_polytope_month(_month_tree()) + assert_valid_covjson(covjson) + + def test_path_geojson_lonlat_order(self): + covjson = Covjsonkit().encode("CoverageCollection", "Path").from_polytope(_surface_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0] + + def test_path_geojson_lonlat_z_order(self): + covjson = Covjsonkit().encode("CoverageCollection", "Path").from_polytope(_level_forecast_tree()) + gj = Covjsonkit().decode(covjson).to_geojson() + assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0, 500] + + @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) + def test_path_legacy_and_new_geojson_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "Path").from_polytope(tree_factory()) + legacy_covjson = _to_legacy_path(new_covjson) + + new_gj = Covjsonkit().decode(new_covjson).to_geojson() + legacy_gj = Covjsonkit().decode(legacy_covjson).to_geojson() + + # Legacy always carries a z component; drop it from the geometry so the + # surface case (spec output has no z) compares equal. + for gj in (new_gj, legacy_gj): + for feat in gj["features"]: + feat["geometry"]["coordinates"] = feat["geometry"]["coordinates"][:2] + + assert legacy_gj == new_gj diff --git a/tests/test_xarray_path.py b/tests/test_xarray_path.py index 8717648..54c860f 100644 --- a/tests/test_xarray_path.py +++ b/tests/test_xarray_path.py @@ -32,10 +32,15 @@ def test_from_xarray(self): assert covjson_result["type"] == self.test_covjson["type"] assert len(covjson_result["coverages"]) == len(self.test_covjson["coverages"]) - assert ( - covjson_result["coverages"][0]["domain"]["axes"]["composite"]["values"] - == self.test_covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] - ) + # Encoder now emits spec-compliant Trajectory composite: coordinates + # ["t","x","y","z"] with values ordered [t, lon, lat, level]. The legacy + # fixture is mislabeled ([t, lat, lon, level]); the decoder reads it + # back-compatibly, so the round-trip swaps lat/lon into spec order. + composite = covjson_result["coverages"][0]["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["t", "x", "y", "z"] + legacy_values = self.test_covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] + expected_values = [[t, lon, lat, z] for (t, lat, lon, z) in legacy_values] + assert composite["values"] == expected_values assert ( covjson_result["coverages"][0]["mars:metadata"]["number"] == self.test_covjson["coverages"][0]["mars:metadata"]["number"] From a48a5b3cc3b8798e90d7cd542732dfd10e42944d Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 4 Sep 2026 15:18:18 +0200 Subject: [PATCH 04/10] ALign grid with ogc spec --- covjsonkit/decoder/Grid.py | 101 +++++--- covjsonkit/encoder/Grid.py | 308 +++++++++-------------- tests/test_encoder_grid_from_polytope.py | 28 ++- tests/test_spec_compliance.py | 132 +++++++++- 4 files changed, 322 insertions(+), 247 deletions(-) diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index 1abf68f..a0cf722 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -46,13 +46,30 @@ def get_coordinates(self): def to_geopandas(self): pass - def to_geotiff(self, output_file="multipoint", resolution=0.01): + def _grid_axes(self, domain): + """Resolve (lon_values, lat_values, z_values_or_None) from a Grid domain, + reading both spec (x/y/z) and legacy (longitude/latitude/levelist) axes.""" + lon_axis = "x" if "x" in domain else "longitude" + lat_axis = "y" if "y" in domain else "latitude" + if "z" in domain: + z_axis = "z" + elif "levelist" in domain: + z_axis = "levelist" + else: + z_axis = None + lon = domain[lon_axis]["values"] + lat = domain[lat_axis]["values"] + z = domain[z_axis]["values"] if z_axis is not None else None + return lon, lat, z + + def to_geotiff(self, output_file="grid", resolution=0.01): if rasterio is None: raise ImportError("Please install 'rasterio' to use this feature: pip install covjsonkit[geo]") - coords = self.covjson["coverages"][0]["domain"]["axes"]["composite"]["values"] - x = [c[1] for c in coords] # longitude - y = [c[0] for c in coords] # latitude - # z = [c[2] for c in coords] # height/time/etc (not used yet) + domain = self.covjson["coverages"][0]["domain"]["axes"] + lon_vals, lat_vals, _ = self._grid_axes(domain) + # Full grid of (lon, lat) cell centres. + x = [lon for lat in lat_vals for lon in lon_vals] + y = [lat for lat in lat_vals for lon in lon_vals] # Define grid x_min, x_max = min(x), max(x) @@ -102,29 +119,36 @@ def to_geotiff(self, output_file="multipoint", resolution=0.01): def to_geojson(self): features = [] for coverage in self.covjson["coverages"]: - coords = coverage["domain"]["axes"]["composite"]["values"] - datetime = coverage["domain"]["axes"]["t"]["values"][0] - if "mars:metadata" in coverage: - mars_metadata = coverage["mars:metadata"] + domain = coverage["domain"]["axes"] + lon_vals, lat_vals, z_vals = self._grid_axes(domain) + datetime = domain["t"]["values"][0] + mars_metadata = coverage.get("mars:metadata") values = {} for key in coverage["ranges"]: values[key] = coverage["ranges"][key]["values"] - for idx, lonlat in enumerate(coords): - param_vals = {} - for key in values.keys(): - param_vals[key] = values[key][idx] - param_vals["datetime"] = datetime - if "mars:metadata" in coverage: - param_vals["mars:metadata"] = mars_metadata - features.append( - { - "type": "Feature", - "geometry": {"type": "Point", "coordinates": [lonlat[1], lonlat[0], lonlat[2]]}, - "properties": param_vals, - } - ) + # Grid points iterate latitude-major then longitude (matching the + # [t, (z), y, x] range value ordering). + idx = 0 + z = z_vals[0] if z_vals else None + for lat in lat_vals: + for lon in lon_vals: + param_vals = {} + for key in values.keys(): + param_vals[key] = values[key][idx] + param_vals["datetime"] = datetime + if mars_metadata is not None: + param_vals["mars:metadata"] = mars_metadata + geom_coords = [lon, lat] if z is None else [lon, lat, z] + features.append( + { + "type": "Feature", + "geometry": {"type": "Point", "coordinates": geom_coords}, + "properties": param_vals, + } + ) + idx += 1 geojson = {"type": "FeatureCollection", "features": features} return geojson @@ -160,23 +184,21 @@ def to_xarray(self): first_cov = self.covjson["coverages"][0] domain = first_cov["domain"]["axes"] - if "latitude" in domain: - x_coords = "latitude" - else: - x_coords = "x" - if "longitude" in domain: - y_coords = "longitude" - else: - y_coords = "y" - if "levelist" in domain: - z_coords = "levelist" + # Resolve named axes with back-compat: spec output uses x/y/z, legacy + # output uses longitude/latitude/levelist. x == longitude, y == latitude. + lon_axis = "x" if "x" in domain else "longitude" + lat_axis = "y" if "y" in domain else "latitude" + if "z" in domain: + z_axis = "z" + elif "levelist" in domain: + z_axis = "levelist" else: - z_coords = "z" + z_axis = None steps = np.array(domain.get("t", {}).get("values", [0])) - levels = np.array(domain.get(z_coords, {}).get("values", [0])) - lat = np.array(domain[x_coords]["values"]) - lon = np.array(domain[y_coords]["values"]) + levels = np.array(domain[z_axis]["values"]) if z_axis is not None else np.array([0]) + lat = np.array(domain[lat_axis]["values"]) + lon = np.array(domain[lon_axis]["values"]) # Prepare arrays for each parameter data_arrays = { @@ -199,6 +221,11 @@ def to_xarray(self): for pname, prange in coverage["ranges"].items(): arr = np.array(prange["values"]).reshape(prange["shape"]) + # Spec-compliant surface grids drop the z axis (shape + # [t, y, x]); insert a length-1 level dim so the target 4D + # (steps, levels, lat, lon) block assignment lines up. + if arr.ndim == 3: + arr = arr[:, np.newaxis, :, :] data_arrays[pname][t_idx, n_idx, :, :, :, :] = arr # Build xarray Dataset diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index 5f3d56d..0d9a821 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -2,8 +2,6 @@ import time from datetime import timedelta -import pandas as pd - from .encoder import Encoder, normalize_step_value @@ -13,40 +11,41 @@ def __init__(self, type, domaintype): self.covjson["domainType"] = "Grid" self.covjson["coverages"] = [] - def add_coverage(self, mars_metadata, coords, values): + def add_coverage(self, mars_metadata, coords, values, include_z=False): new_coverage = {} new_coverage["mars:metadata"] = {} new_coverage["type"] = "Coverage" new_coverage["domain"] = {} new_coverage["ranges"] = {} self.add_mars_metadata(new_coverage, mars_metadata) - self.add_domain(new_coverage, coords) - self.add_range(new_coverage, values) + self.add_domain(new_coverage, coords, include_z) + self.add_range(new_coverage, values, include_z) self.covjson["coverages"].append(new_coverage) - # cov = Coverage.model_validate_json(json.dumps(new_coverage)) - # self.pydantic_coverage.coverages.append(json.dumps(new_coverage)) - def add_domain(self, coverage, coords): + def add_domain(self, coverage, coords, include_z=False): + # OGC-compliant Grid domain: regularly-sampled named axes ``x`` + # (longitude), ``y`` (latitude), ``t`` (time) and an optional ``z`` + # (level). ``z`` is omitted at the surface. Coordinate values may be + # provided under the spec keys (x/y/z) or the legacy keys + # (longitude/latitude/levelist). + axes = {"t": {"values": coords["t"]}} + if include_z: + axes["z"] = {"values": coords.get("z", coords.get("levelist"))} + axes["y"] = {"values": coords.get("y", coords.get("latitude"))} + axes["x"] = {"values": coords.get("x", coords.get("longitude"))} coverage["domain"]["type"] = "Domain" - coverage["domain"]["axes"] = {} - coverage["domain"]["axes"]["t"] = {} - coverage["domain"]["axes"]["t"]["values"] = coords["t"] - coverage["domain"]["axes"]["latitude"] = {} - coverage["domain"]["axes"]["latitude"]["values"] = coords["latitude"] - coverage["domain"]["axes"]["longitude"] = {} - coverage["domain"]["axes"]["longitude"]["values"] = coords["longitude"] - coverage["domain"]["axes"]["levelist"] = {} - coverage["domain"]["axes"]["levelist"]["values"] = coords["levelist"] - - def add_range(self, coverage, values): + coverage["domain"]["axes"] = axes + + def add_range(self, coverage, values, include_z=False): + axis_names = ["t", "z", "y", "x"] if include_z else ["t", "y", "x"] for parameter in values.keys(): param = self.convert_param_id_to_param(parameter) coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" coverage["ranges"][param]["shape"] = self.shp - coverage["ranges"][param]["axisNames"] = ["t", "levelist", "latitude", "longitude"] - coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] + coverage["ranges"][param]["axisNames"] = axis_names + coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): coverage["mars:metadata"] = metadata @@ -55,50 +54,61 @@ def add_if_not_close(self, my_list, number, threshold=0.01): if all(abs(number - x) > threshold for x in my_list): my_list.append(number) + def _set_references(self, include_z): + # Grid referencing is split into GeographicCRS (x, y), a TemporalRS (t) + # and an optional VerticalCRS (z) so each coordinate carries the correct + # coordinate reference system. + refs = [ + { + "coordinates": ["x", "y"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, + ] + if include_z: + refs.append( + { + "coordinates": ["z"], + "system": {"type": "VerticalCRS"}, + } + ) + self.covjson["referencing"] = refs + def from_xarray(self, dataset): """ - Converts an xarray dataset into a grid CoverageJSON format. + Converts an xarray dataset into a Grid CoverageJSON collection. """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "Grid" self.covjson["coverages"] = [] - if "latitude" in dataset.coords: - x_coord = "latitude" - elif "x" in dataset.coords: - x_coord = "x" - if "longitude" in dataset.coords: - y_coord = "longitude" - elif "y" in dataset.coords: - y_coord = "y" - if "levelist" in dataset.coords: - z_coord = "levelist" - - # Add reference system - self.add_reference( - { - "coordinates": [x_coord, y_coord, z_coord], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = "levelist" in dataset.coords + + # Add reference system (split GeographicCRS / TemporalRS / VerticalCRS) + self._set_references(include_z) for data_var in dataset.data_vars: data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates + # Prepare coordinates (spec-named axes: x=lon, y=lat, z=level). coords = { "t": [str(x) for x in dataset["steps"].values], - "latitude": dataset["latitude"].values.tolist(), - "longitude": dataset["longitude"].values.tolist(), - "levelist": dataset["levelist"].values.tolist(), + "y": dataset["latitude"].values.tolist(), + "x": dataset["longitude"].values.tolist(), } - - self.shp = [len(coords["t"]), len(coords["levelist"]), len(coords["latitude"]), len(coords["longitude"])] + if include_z: + coords["z"] = dataset["levelist"].values.tolist() + self.shp = [len(coords["t"]), len(coords["z"]), len(coords["y"]), len(coords["x"])] + else: + self.shp = [len(coords["t"]), len(coords["y"]), len(coords["x"])] for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -110,13 +120,11 @@ def from_xarray(self, dataset): mars_metadata["Forecast date"] = str(datetime) for dv in dataset.data_vars: nested_list = dataset[dv].sel(datetimes=datetime, number=num, steps=step).values.tolist() - print(nested_list) flattened_list = [item for sublist in nested_list for item in sublist] flattened_list = [item for sublist in flattened_list for item in sublist] - print(flattened_list) dv_dict[dv] = flattened_list - self.add_coverage(mars_metadata, coords, dv_dict) + self.add_coverage(mars_metadata, coords, dv_dict, include_z) # Return the generated CoverageJSON return self.covjson @@ -133,20 +141,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: fields["step"] = [0] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False self.walk_tree(result, fields, coords, mars_metadata, range_dict, date_key=date_key) logging.debug("The values returned from walking tree: %s", range_dict) # noqa: E501 logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) combined_dict = {} @@ -160,15 +162,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for para in fields["param"]: if para not in combined_dict[date][num]: combined_dict[date][num][para] = {} - # for s, value in range_dict[date][level][num][para].items(): for s in fields["step"]: key = (date, level, num, para, s) - # for k, v in range_dict.items(): - # if k == key: if s not in combined_dict[date][num][para]: combined_dict[date][num][para][s] = range_dict[key] else: - # Cocatenate arrays + # Concatenate arrays combined_dict[date][num][para][s] += range_dict[key] if fields["param"] == 0: @@ -181,29 +180,36 @@ def from_polytope(self, result, date_key: str = "date") -> dict: logging.debug("The fields retrieved were: %s", fields) # noqa: E501 logging.debug("The range_dict created was: %s", range_dict) # noqa: E501 + # Build named grid axes: y (latitude) and x (longitude) are the unique + # sorted coordinate values; walk_tree stores each point as [lat, lon]. coordinates = {} - coordinates["t"] = [int(s.total_seconds() // 3600) if isinstance(s, timedelta) else s for s in fields["step"]] - for date in coords.keys(): coordinates[date] = {} coordinates[date]["t"] = [ int(s.total_seconds() // 3600) if isinstance(s, timedelta) else s for s in fields["step"] ] - coordinates[date]["levelist"] = list(fields["levels"]) - coordinates[date]["latitude"] = [] - coordinates[date]["longitude"] = [] + if include_z: + coordinates[date]["z"] = list(fields["levels"]) + coordinates[date]["y"] = [] + coordinates[date]["x"] = [] for cor in coords[date]["composite"]: - self.add_if_not_close(coordinates[date]["latitude"], cor[0]) - self.add_if_not_close(coordinates[date]["longitude"], cor[1]) - coordinates[date]["latitude"] = list(coordinates[date]["latitude"]) - coordinates[date]["longitude"] = list(coordinates[date]["longitude"]) - - self.shp = [ - len(coordinates[fields["dates"][0]]["t"]), - len(coordinates[fields["dates"][0]]["levelist"]), - len(coordinates[fields["dates"][0]]["latitude"]), - len(coordinates[fields["dates"][0]]["longitude"]), - ] + self.add_if_not_close(coordinates[date]["y"], cor[0]) # latitude + self.add_if_not_close(coordinates[date]["x"], cor[1]) # longitude + + first_date = fields["dates"][0] + if include_z: + self.shp = [ + len(coordinates[first_date]["t"]), + len(coordinates[first_date]["z"]), + len(coordinates[first_date]["y"]), + len(coordinates[first_date]["x"]), + ] + else: + self.shp = [ + len(coordinates[first_date]["t"]), + len(coordinates[first_date]["y"]), + len(coordinates[first_date]["x"]), + ] for date in combined_dict.keys(): for num in combined_dict[date].keys(): @@ -216,7 +222,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, coordinates[date], val_dict) + self.add_coverage(mm, coordinates[date], val_dict, include_z) return self.covjson @@ -232,6 +238,7 @@ def from_polytope_month(self, result): fields["months"] = [] fields["dates"] = [] fields["levels"] = [0] + fields["has_level_axis"] = False start = time.time() logging.debug("Tree walking starts at: %s", start) @@ -243,15 +250,8 @@ def from_polytope_month(self, result): start = time.time() logging.debug("Coords creation: %s", start) - self.add_reference( - { - "coordinates": ["latitude", "longitude", "levelist"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) + include_z = fields["has_level_axis"] + self._set_references(include_z) if fields["param"] == 0: raise ValueError("No data was returned.") @@ -268,21 +268,29 @@ def from_polytope_month(self, result): for date in fields["dates"]: coordinates[date] = {} coordinates[date]["t"] = [f"{date}-01T00:00:00Z"] - coordinates[date]["levelist"] = list(fields["levels"]) - coordinates[date]["latitude"] = [] - coordinates[date]["longitude"] = [] + if include_z: + coordinates[date]["z"] = list(fields["levels"]) + coordinates[date]["y"] = [] + coordinates[date]["x"] = [] for cor in coords.get(date, {}).get("composite", []): - self.add_if_not_close(coordinates[date]["latitude"], cor[0]) - self.add_if_not_close(coordinates[date]["longitude"], cor[1]) + self.add_if_not_close(coordinates[date]["y"], cor[0]) # latitude + self.add_if_not_close(coordinates[date]["x"], cor[1]) # longitude if fields["dates"]: first_date = fields["dates"][0] - self.shp = [ - len(coordinates[first_date]["t"]), - len(coordinates[first_date]["levelist"]), - len(coordinates[first_date]["latitude"]), - len(coordinates[first_date]["longitude"]), - ] + if include_z: + self.shp = [ + len(coordinates[first_date]["t"]), + len(coordinates[first_date]["z"]), + len(coordinates[first_date]["y"]), + len(coordinates[first_date]["x"]), + ] + else: + self.shp = [ + len(coordinates[first_date]["t"]), + len(coordinates[first_date]["y"]), + len(coordinates[first_date]["x"]), + ] end = time.time() logging.debug("Coords creation: %s", end) @@ -303,102 +311,10 @@ def from_polytope_month(self, result): mm = mars_metadata.copy() mm["number"] = num if fields["dates"]: - self.add_coverage(mm, coordinates[fields["dates"][0]], val_dict) + self.add_coverage(mm, coordinates[fields["dates"][0]], val_dict, include_z) end = time.time() logging.debug("Coverage creation: %s", end) logging.debug("Coverage creation takes: %s", end - start) return self.covjson - coords = {} - mars_metadata = {} - range_dict = {} - fields = {} - fields["lat"] = 0 - fields["param"] = 0 - fields["number"] = [0] - fields["step"] = [0] - fields["dates"] = [] - fields["levels"] = [0] - fields["times"] = [] - - start = time.time() - logging.debug("Tree walking starts at: %s", start) # noqa: E501 - self.walk_tree_step(result, fields, coords, mars_metadata, range_dict) - end = time.time() - delta = end - start - logging.debug("Tree walking ends at: %s", end) # noqa: E501 - logging.debug("Tree walking takes: %s", delta) # noqa: E501 - - start = time.time() - logging.debug("Coords creation: %s", start) # noqa: E501 - - self.add_reference( - { - "coordinates": ["x", "y", "z"], - "system": { - "type": "GeographicCRS", - "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", - }, - } - ) - - coordinates = {} - - levels = fields["levels"] - if fields["param"] == 0: - raise ValueError("No data was returned.") - for para in fields["param"]: - self.add_parameter(para) - - logging.debug("The parameters added were: %s", self.parameters) # noqa: E501 - - # points = len(coords[fields["dates"][0]]["composite"]) - - coordinates = {} - - for date in coords.keys(): - for dt in coords[date]["t"]: - coordinates[dt] = {} - coordinates[dt]["composite"] = [] - coordinates[dt]["t"] = [dt] - coord = coords[date]["composite"] - for level in levels: - for cor in coord: - coordinates[dt]["composite"].append([cor[0], cor[1], level]) - - end = time.time() - delta = end - start - logging.debug("Coords creation: %s", end) # noqa: E501 - logging.debug("Coords creation: %s", delta) # noqa: E501 - - start = time.time() - logging.debug("Coverage creation: %s", start) # noqa: E501 - - val_dict = {} - for i, t in enumerate(fields["times"]): - val_dict[t] = {} - for level in fields["levels"]: - for num in fields["number"]: - for date in fields["dates"]: - for para in fields["param"]: - val_dict[t][para] = [] - key = (date, level, num, para) - vals = int(len(range_dict[key]) / len(fields["times"])) - # for val in range_dict[key]: - # vals.append(val[i]) - val_dict[t][para].extend(range_dict[key][i * vals : (i + 1) * vals]) - for para in fields["param"]: - val_dict[t][para] = [item for sublist in val_dict[t][para] for item in sublist] - mm = mars_metadata.copy() - mm["number"] = num - # mm["Forecast date"] = date - datetime = pd.Timestamp(date) + t - self.add_coverage(mm, coordinates[str(datetime).split("+")[0] + "Z"], val_dict[t]) - - end = time.time() - delta = end - start - logging.debug("Coverage creation: %s", end) # noqa: E501 - logging.debug("Coverage creation: %s", delta) # noqa: E501 - - return self.covjson diff --git a/tests/test_encoder_grid_from_polytope.py b/tests/test_encoder_grid_from_polytope.py index 146f3ec..a0c362e 100644 --- a/tests/test_encoder_grid_from_polytope.py +++ b/tests/test_encoder_grid_from_polytope.py @@ -14,17 +14,16 @@ GRID_2X2_AXES = { "t": {"values": [0]}, - "latitude": {"values": [48.0, 50.0]}, - "longitude": {"values": [11.0, 12.0]}, - "levelist": {"values": [0]}, + "y": {"values": [48.0, 50.0]}, + "x": {"values": [11.0, 12.0]}, } GRID_2X2_RANGES = { "2t": { "type": "NdArray", "dataType": "float", - "shape": [1, 1, 2, 2], - "axisNames": ["t", "levelist", "latitude", "longitude"], + "shape": [1, 2, 2], + "axisNames": ["t", "y", "x"], "values": [264.9, 265.1, 266.3, 267.5], } } @@ -65,15 +64,19 @@ def test_2x2_grid(self): assert covjson["type"] == "CoverageCollection" assert covjson["domainType"] == "Grid" - # Collection-level referencing + # Collection-level referencing (split GeographicCRS / TemporalRS) assert covjson["referencing"] == [ { - "coordinates": ["latitude", "longitude", "levelist"], + "coordinates": ["x", "y"], "system": { "type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", }, - } + }, + { + "coordinates": ["t"], + "system": {"type": "TemporalRS", "calendar": "Gregorian"}, + }, ] # Collection-level parameters @@ -124,16 +127,15 @@ def test_1x1_grid(self): assert cov["domain"]["axes"] == { "t": {"values": [0]}, - "latitude": {"values": [48.0]}, - "longitude": {"values": [11.0]}, - "levelist": {"values": [0]}, + "y": {"values": [48.0]}, + "x": {"values": [11.0]}, } assert cov["ranges"] == { "2t": { "type": "NdArray", "dataType": "float", - "shape": [1, 1, 1, 1], - "axisNames": ["t", "levelist", "latitude", "longitude"], + "shape": [1, 1, 1], + "axisNames": ["t", "y", "x"], "values": [264.9], } } diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index 5e3cc93..4b55c76 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -17,7 +17,7 @@ import numpy as np import pytest -from conftest import assert_valid_covjson, chain, make_leaf, make_point, node +from conftest import assert_valid_covjson, chain, make_leaf, make_point, node, tip from polytope_feature.datacube.tensor_index_tree import TensorIndexTree from covjsonkit.api import Covjsonkit @@ -495,3 +495,133 @@ def test_path_legacy_and_new_geojson_equivalent(self, tree_factory): feat["geometry"]["coordinates"] = feat["geometry"]["coordinates"][:2] assert legacy_gj == new_gj + + +def _grid_forecast_tree(): + """Single-date surface 2x2 grid (analysis) for Grid spec-compliance tests.""" + tree = chain( + TensorIndexTree(), + node("class", ("od",)), + node("date", (np.datetime64("2025-01-01T00:00:00"),)), + node("domain", ("g",)), + node("expver", ("0001",)), + node("levtype", ("sfc",)), + node("param", ("167",)), + node("step", (0,)), + node("stream", ("oper",)), + node("type", ("an",)), + ) + parent = tip(tree) + for lat, lon, vals in [ + (48.0, 11.0, [264.9]), + (48.0, 12.0, [265.1]), + (50.0, 11.0, [266.3]), + (50.0, 12.0, [267.5]), + ]: + parent.add_child(make_point(lat, lon, vals)) + return tree + + +def _grid_level_tree(): + """Single-date pressure-level 2x2 grid (analysis) for Grid spec tests.""" + tree = chain( + TensorIndexTree(), + node("class", ("od",)), + node("date", (np.datetime64("2025-01-01T00:00:00"),)), + node("domain", ("g",)), + node("expver", ("0001",)), + node("levtype", ("pl",)), + node("levelist", (500,)), + node("param", ("129",)), + node("step", (0,)), + node("stream", ("oper",)), + node("type", ("an",)), + ) + parent = tip(tree) + for lat, lon, vals in [ + (48.0, 11.0, [264.9]), + (48.0, 12.0, [265.1]), + (50.0, 11.0, [266.3]), + (50.0, 12.0, [267.5]), + ]: + parent.add_child(make_point(lat, lon, vals)) + return tree + + +def _to_legacy_grid(covjson): + """Convert a spec-compliant Grid covjson into the legacy form: named axes + latitude/longitude/levelist (always including a z axis) with range axisNames + ["t", "levelist", "latitude", "longitude"] and a single combined + GeographicCRS. Used to assert the decoder reads both.""" + legacy = copy.deepcopy(covjson) + legacy["referencing"] = [ + { + "coordinates": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ] + for coverage in legacy["coverages"]: + axes = coverage["domain"]["axes"] + has_z = "z" in axes + new_axes = {"t": axes["t"]} + new_axes["levelist"] = axes["z"] if has_z else {"values": [0]} + new_axes["latitude"] = axes["y"] + new_axes["longitude"] = axes["x"] + coverage["domain"]["axes"] = new_axes + for prange in coverage["ranges"].values(): + prange["axisNames"] = ["t", "levelist", "latitude", "longitude"] + if not has_z: + # insert the length-1 level dim into the shape + t, y, x = prange["shape"] + prange["shape"] = [t, 1, y, x] + return legacy + + +class TestGridSpecCompliance: + """Grid encoder output validates and decodes with back-compat. + + Spec Grid uses named axes x (lon), y (lat), t and an optional z (level), + dropping z at the surface, with split referencing. The legacy form uses + latitude/longitude/levelist axes (always with z) and a single combined + GeographicCRS. The decoder reads both.""" + + @pytest.mark.parametrize("tree_factory", [_grid_forecast_tree, _grid_level_tree]) + def test_grid_output_validates(self, tree_factory): + covjson = Covjsonkit().encode("CoverageCollection", "Grid").from_polytope(tree_factory()) + assert_valid_covjson(covjson) + + def test_grid_surface_drops_z(self): + covjson = Covjsonkit().encode("CoverageCollection", "Grid").from_polytope(_grid_forecast_tree()) + axes = covjson["coverages"][0]["domain"]["axes"] + assert set(axes) == {"t", "y", "x"} + assert covjson["coverages"][0]["ranges"]["2t"]["axisNames"] == ["t", "y", "x"] + + def test_grid_level_includes_z(self): + covjson = Covjsonkit().encode("CoverageCollection", "Grid").from_polytope(_grid_level_tree()) + axes = covjson["coverages"][0]["domain"]["axes"] + assert set(axes) == {"t", "z", "y", "x"} + assert axes["z"]["values"] == [500] + param_range = next(iter(covjson["coverages"][0]["ranges"].values())) + assert param_range["axisNames"] == ["t", "z", "y", "x"] + + def test_grid_axes_lonlat_values(self): + # x holds longitude, y holds latitude. + covjson = Covjsonkit().encode("CoverageCollection", "Grid").from_polytope(_grid_forecast_tree()) + axes = covjson["coverages"][0]["domain"]["axes"] + assert axes["x"]["values"] == [11.0, 12.0] + assert axes["y"]["values"] == [48.0, 50.0] + + @pytest.mark.parametrize("tree_factory", [_grid_forecast_tree, _grid_level_tree]) + def test_grid_legacy_and_new_xarray_equivalent(self, tree_factory): + new_covjson = Covjsonkit().encode("CoverageCollection", "Grid").from_polytope(tree_factory()) + legacy_covjson = _to_legacy_grid(new_covjson) + + new_ds = Covjsonkit().decode(new_covjson).to_xarray() + legacy_ds = Covjsonkit().decode(legacy_covjson).to_xarray() + + assert new_ds.identical(legacy_ds) + assert "latitude" in new_ds.coords + assert "longitude" in new_ds.coords From 9358566753b9c238cb85143f4cbe8c1d5017eaab Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 4 Sep 2026 15:28:07 +0200 Subject: [PATCH 05/10] Fix bug when converting grid feature to_xarray --- covjsonkit/encoder/Grid.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index 0d9a821..a0e1f11 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -98,17 +98,16 @@ def from_xarray(self, dataset): data_var = self.convert_param_to_param_id(data_var) self.add_parameter(data_var) - # Prepare coordinates (spec-named axes: x=lon, y=lat, z=level). - coords = { - "t": [str(x) for x in dataset["steps"].values], - "y": dataset["latitude"].values.tolist(), - "x": dataset["longitude"].values.tolist(), - } + # Prepare shared spatial coordinates (spec-named axes: x=lon, y=lat, + # z=level). Each coverage carries a single time step, so the ``t`` axis + # and range shape use a t-dimension of 1 (set per coverage below). + y_vals = dataset["latitude"].values.tolist() + x_vals = dataset["longitude"].values.tolist() + z_vals = dataset["levelist"].values.tolist() if include_z else None if include_z: - coords["z"] = dataset["levelist"].values.tolist() - self.shp = [len(coords["t"]), len(coords["z"]), len(coords["y"]), len(coords["x"])] + self.shp = [1, len(z_vals), len(y_vals), len(x_vals)] else: - self.shp = [len(coords["t"]), len(coords["y"]), len(coords["x"])] + self.shp = [1, len(y_vals), len(x_vals)] for datetime in dataset["datetimes"].values: for num in dataset["number"].values: @@ -119,11 +118,15 @@ def from_xarray(self, dataset): mars_metadata["step"] = normalize_step_value(step) mars_metadata["Forecast date"] = str(datetime) for dv in dataset.data_vars: - nested_list = dataset[dv].sel(datetimes=datetime, number=num, steps=step).values.tolist() - flattened_list = [item for sublist in nested_list for item in sublist] - flattened_list = [item for sublist in flattened_list for item in sublist] - dv_dict[dv] = flattened_list - + # Flatten in axis order z, y, x (or y, x at the surface). + # ravel() handles both the 3D (level) and 2D (surface) + # cases without assuming a level dimension. + arr = dataset[dv].sel(datetimes=datetime, number=num, steps=step).values + dv_dict[dv] = arr.ravel().tolist() + + coords = {"t": [str(step)], "y": y_vals, "x": x_vals} + if include_z: + coords["z"] = z_vals self.add_coverage(mars_metadata, coords, dv_dict, include_z) # Return the generated CoverageJSON From 46f98f64b7d5b914554b85f2fb42037eef2daf36 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Mon, 7 Sep 2026 13:44:28 +0200 Subject: [PATCH 06/10] Convert all t values in coverages to be valid_time, remove forecast_date from metadata for stream efcl --- covjsonkit/encoder/BoundingBox.py | 17 +++++++-- covjsonkit/encoder/Circle.py | 14 +++++-- covjsonkit/encoder/Frame.py | 14 +++++-- covjsonkit/encoder/Grid.py | 18 +++++---- covjsonkit/encoder/Path.py | 25 +++++++----- covjsonkit/encoder/Position.py | 11 ++++-- covjsonkit/encoder/Shapefile.py | 14 +++++-- covjsonkit/encoder/VerticalProfile.py | 13 +++++-- covjsonkit/encoder/Wkt.py | 14 +++++-- covjsonkit/encoder/encoder.py | 37 ++++++++++++++++++ tests/conftest.py | 1 - ...test_encoder_bounding_box_from_polytope.py | 38 +++++++------------ tests/test_encoder_circle_from_polytope.py | 4 +- tests/test_encoder_frame_from_polytope.py | 10 +---- tests/test_encoder_grid_from_polytope.py | 31 +++++++-------- tests/test_encoder_path_from_polytope.py | 30 ++++++++------- tests/test_encoder_position_from_polytope.py | 3 +- tests/test_encoder_shapefile_from_polytope.py | 4 +- ..._encoder_vertical_profile_from_polytope.py | 11 +----- tests/test_encoder_wkt_from_polytope.py | 7 +--- tests/test_spec_compliance.py | 5 +++ 21 files changed, 197 insertions(+), 124 deletions(-) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index 7fc169b..4381c9b 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -1,7 +1,7 @@ import logging import time -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class BoundingBox(Encoder): @@ -188,6 +188,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: else: coords[date]["composite"].append([cor[1], cor[0]]) + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -200,8 +202,17 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step], include_z) + # ``t`` is the valid-time (forecast date + step offset). + cov_coords = dict(coords[date]) + cov_coords["t"] = [valid_time(date, step)] + if reanalysis: + # Reanalysis (class=ce, stream=efcl): expose only the + # valid-time; drop the scalar forecast-date metadata. + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + self.add_coverage(mm, cov_coords, val_dict[step], include_z) return self.covjson diff --git a/covjsonkit/encoder/Circle.py b/covjsonkit/encoder/Circle.py index 67c7700..eee67f2 100644 --- a/covjsonkit/encoder/Circle.py +++ b/covjsonkit/encoder/Circle.py @@ -1,6 +1,6 @@ import logging -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Circle(Encoder): @@ -176,6 +176,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: else: coords[date]["composite"].append([cor[1], cor[0]]) + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -188,8 +190,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step], include_z) + cov_coords = dict(coords[date]) + cov_coords["t"] = [valid_time(date, step)] + if reanalysis: + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + self.add_coverage(mm, cov_coords, val_dict[step], include_z) return self.covjson diff --git a/covjsonkit/encoder/Frame.py b/covjsonkit/encoder/Frame.py index 309d41f..67628ac 100644 --- a/covjsonkit/encoder/Frame.py +++ b/covjsonkit/encoder/Frame.py @@ -1,6 +1,6 @@ import logging -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Frame(Encoder): @@ -173,6 +173,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: else: coords[date]["composite"].append([cor[1], cor[0]]) + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -185,8 +187,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step], include_z) + cov_coords = dict(coords[date]) + cov_coords["t"] = [valid_time(date, step)] + if reanalysis: + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + self.add_coverage(mm, cov_coords, val_dict[step], include_z) return self.covjson diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index a0e1f11..b368434 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -1,8 +1,7 @@ import logging import time -from datetime import timedelta -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Grid(Encoder): @@ -188,9 +187,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: coordinates = {} for date in coords.keys(): coordinates[date] = {} - coordinates[date]["t"] = [ - int(s.total_seconds() // 3600) if isinstance(s, timedelta) else s for s in fields["step"] - ] + # ``t`` axis holds the valid-time (forecast date + step) for each step. + coordinates[date]["t"] = [valid_time(date, s) for s in fields["step"]] if include_z: coordinates[date]["z"] = list(fields["levels"]) coordinates[date]["y"] = [] @@ -214,6 +212,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: len(coordinates[first_date]["x"]), ] + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -223,8 +223,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: val_dict[para].extend(combined_dict[date][num][para][step]) mm = mars_metadata.copy() mm["number"] = num - mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date + if reanalysis: + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["step"] = normalize_step_value(step) + mm["Forecast date"] = date self.add_coverage(mm, coordinates[date], val_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/Path.py b/covjsonkit/encoder/Path.py index fe91f3f..a42c9bc 100644 --- a/covjsonkit/encoder/Path.py +++ b/covjsonkit/encoder/Path.py @@ -1,6 +1,6 @@ import logging -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Path(Encoder): @@ -164,22 +164,21 @@ def from_polytope(self, result, date_key: str = "date") -> dict: if (date, level, fields["number"][0], fields["param"][0], s) in range_dict: cor_len = len(range_dict[(date, level, fields["number"][0], fields["param"][0], s)]) end = start + cor_len - # Normalize the step value before it is placed into the - # composite (t) coordinate so that timedelta/Timedelta - # values are JSON serialisable. The original ``s`` is kept - # for the range_dict lookups above. - s_norm = normalize_step_value(s) + # The composite ``t`` coordinate is the valid-time + # (forecast date + step offset), as an ISO-8601 string. + # The original ``s`` is kept for the range_dict lookups. + t_valid = valid_time(date, s) if include_z: if len(fields["levels"]) != 1: for lev in fields["levels"]: for cor in coord[int(start) : int(end)]: - coords[date]["composite"].append([s_norm, cor[1], cor[0], lev]) + coords[date]["composite"].append([t_valid, cor[1], cor[0], lev]) else: for cor in coord[int(start) : int(end)]: - coords[date]["composite"].append([s_norm, cor[1], cor[0], level]) + coords[date]["composite"].append([t_valid, cor[1], cor[0], level]) else: for cor in coord[int(start) : int(end)]: - coords[date]["composite"].append([s_norm, cor[1], cor[0]]) + coords[date]["composite"].append([t_valid, cor[1], cor[0]]) start = end logging.debug("The coordinates returned from walking tree: %s", coords) # noqa: E501 @@ -227,9 +226,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: val_dict[para].extend(combined_dict[date][num][para][step]) mm = mars_metadata.copy() mm["number"] = num - mm["Forecast date"] = date if "levelist" in mm: del mm["levelist"] + if is_reanalysis(mars_metadata, date_key): + # Reanalysis (class=ce, stream=efcl): expose only the + # valid-time; drop the scalar forecast-date metadata. + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date self.add_coverage(mm, coords[date], val_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index 2068ac7..371aab2 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, is_reanalysis class Position(Encoder): @@ -244,10 +244,15 @@ def from_polytope(self, result, date_key: str = "date") -> dict: ) mm = mars_metadata.copy() mm["number"] = num - mm["Forecast date"] = date mm["levelist"] = level coordinates[date][i]["levelist"] = [level] - del mm["step"] + mm.pop("step", None) + if is_reanalysis(mars_metadata, date_key): + # Reanalysis (class=ce, stream=efcl): expose only the + # valid-time; drop the scalar forecast-date metadata. + mm.pop("Forecast date", None) + else: + mm["Forecast date"] = date self.add_coverage(mm, coordinates[date][i], val_dict, include_z) end = time.time() diff --git a/covjsonkit/encoder/Shapefile.py b/covjsonkit/encoder/Shapefile.py index 26b0a50..6be6dbe 100644 --- a/covjsonkit/encoder/Shapefile.py +++ b/covjsonkit/encoder/Shapefile.py @@ -1,6 +1,6 @@ import logging -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Shapefile(Encoder): @@ -173,6 +173,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: else: coords[date]["composite"].append([cor[1], cor[0]]) + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -185,8 +187,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step], include_z) + cov_coords = dict(coords[date]) + cov_coords["t"] = [valid_time(date, step)] + if reanalysis: + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + self.add_coverage(mm, cov_coords, val_dict[step], include_z) return self.covjson diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 2029570..d3d5825 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -4,7 +4,7 @@ import pandas as pd -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value class VerticalProfile(Encoder): @@ -221,9 +221,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: ) mm = mars_metadata.copy() mm["number"] = num - mm["Forecast date"] = date - mm["step"] = normalize_step_value(step) - # del mm["step"] + if is_reanalysis(mars_metadata, date_key): + # Reanalysis (class=ce, stream=efcl): expose only the + # valid-time; drop the scalar forecast-date metadata. + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + mm["step"] = normalize_step_value(step) self.add_coverage(mm, coordinates[date][i][s], val_dict[step]) end = time.time() diff --git a/covjsonkit/encoder/Wkt.py b/covjsonkit/encoder/Wkt.py index 86f84c3..cd552cb 100644 --- a/covjsonkit/encoder/Wkt.py +++ b/covjsonkit/encoder/Wkt.py @@ -4,7 +4,7 @@ import pandas as pd -from .encoder import Encoder, normalize_step_value +from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time class Wkt(Encoder): @@ -180,6 +180,8 @@ def from_polytope(self, result, date_key: str = "date") -> dict: else: coords[date]["composite"].append([cor[1], cor[0]]) + reanalysis = is_reanalysis(mars_metadata, date_key) + for date in combined_dict.keys(): for num in combined_dict[date].keys(): val_dict = {} @@ -192,8 +194,14 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm = mars_metadata.copy() mm["number"] = num mm["step"] = normalize_step_value(step) - mm["Forecast date"] = date - self.add_coverage(mm, coords[date], val_dict[step], include_z) + cov_coords = dict(coords[date]) + cov_coords["t"] = [valid_time(date, step)] + if reanalysis: + mm.pop("Forecast date", None) + mm.pop("step", None) + else: + mm["Forecast date"] = date + self.add_coverage(mm, cov_coords, val_dict[step], include_z) return self.covjson diff --git a/covjsonkit/encoder/encoder.py b/covjsonkit/encoder/encoder.py index 52f06ef..6330ef3 100644 --- a/covjsonkit/encoder/encoder.py +++ b/covjsonkit/encoder/encoder.py @@ -196,6 +196,43 @@ def normalize_step_value(step): return str(step) +def valid_time(date, step): + """Return the ISO-8601 valid-time string for a forecast ``date`` and ``step``. + + The valid time is ``date + step``. ``date`` may be any value accepted by + ``pandas.Timestamp`` (e.g. a MARS ``"YYYYMMDDTHHMMSS"`` string, optionally + suffixed with ``"Z"``). ``step`` may be a ``timedelta``/``numpy.timedelta64`` + or an hour count (int, float, numeric string, or a single-element list). + + The result is an ISO-8601 string with a trailing ``"Z"``, matching the + valid-time strings produced by the PointSeries encoders. + """ + date_format = "%Y%m%dT%H%M%S" + new_date = pd.Timestamp(str(date).rstrip("Z")).strftime(date_format) + start_time = pd.Timestamp(new_date).to_pydatetime() + if isinstance(step, (timedelta, pd.Timedelta, np.timedelta64)): + stamp = start_time + pd.Timedelta(step).to_pytimedelta() + else: + if isinstance(step, (list, tuple)): + step = step[0] + try: + hours = int(step) + except (TypeError, ValueError): + hours = int(str(step).split("h")[0]) + stamp = start_time + timedelta(hours=hours) + return stamp.isoformat() + "Z" + + +def is_reanalysis(mars_metadata: dict, date_key: str = "date") -> bool: + """True for the collapsed reanalysis path (``class=ce``, ``stream=efcl``). + + When walking the ``hdate`` axis for ``class=ce``/``stream=efcl`` requests the + encoders omit the scalar ``"Forecast date"`` metadata and expose only the + valid-time on the ``t`` axis, mirroring the PointSeries behaviour. + """ + return date_key == "hdate" and mars_metadata.get("class") == "ce" and mars_metadata.get("stream") == "efcl" + + class Encoder(ABC): def __init__(self, type, domaintype): """ diff --git a/tests/conftest.py b/tests/conftest.py index 1d8b4cb..5f0c043 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -48,7 +48,6 @@ def assert_valid_covjson(covjson): "domain": "g", "expver": "4321", "levtype": "sfc", - "step": 0, "stream": "efcl", "type": "sfo", "number": 0, diff --git a/tests/test_encoder_bounding_box_from_polytope.py b/tests/test_encoder_bounding_box_from_polytope.py index 1df0548..ef207f5 100644 --- a/tests/test_encoder_bounding_box_from_polytope.py +++ b/tests/test_encoder_bounding_box_from_polytope.py @@ -107,15 +107,15 @@ def test_two_dates_two_steps_two_points(self): } expected = [ - ("2025-01-01T00:00:00Z", 0, [264.9, 265.1]), - ("2025-01-01T00:00:00Z", 6, [270.1, 271.3]), - ("2025-01-02T00:00:00Z", 0, [266.0, 267.0]), - ("2025-01-02T00:00:00Z", 6, [272.0, 273.0]), + ("2025-01-01T00:00:00Z", "2025-01-01T00:00:00Z", 0, [264.9, 265.1]), + ("2025-01-01T00:00:00Z", "2025-01-01T06:00:00Z", 6, [270.1, 271.3]), + ("2025-01-02T00:00:00Z", "2025-01-02T00:00:00Z", 0, [266.0, 267.0]), + ("2025-01-02T00:00:00Z", "2025-01-02T06:00:00Z", 6, [272.0, 273.0]), ] assert len(covjson["coverages"]) == len(expected) - for cov, (date, step, vals) in zip(covjson["coverages"], expected): + for cov, (date, valid, step, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "t": {"values": [date]}, + "t": {"values": [valid]}, "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == vals @@ -150,11 +150,7 @@ def test_reforecast_single_hdate_two_points(self): } } - assert cov["mars:metadata"] == { - **EXPECTED_REFORECAST_METADATA, - "Forecast date": "2025-07-14T06:00:00Z", - "step": 0, - } + assert cov["mars:metadata"] == EXPECTED_REFORECAST_METADATA def test_reforecast_two_hdates_two_points(self): tree = reforecast_tree( @@ -177,11 +173,7 @@ def test_reforecast_two_hdates_two_points(self): "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == vals - assert cov["mars:metadata"] == { - **EXPECTED_REFORECAST_METADATA, - "Forecast date": fc_date, - "step": 0, - } + assert cov["mars:metadata"] == EXPECTED_REFORECAST_METADATA def test_reforecast_single_hdate_two_steps_two_points(self): tree = reforecast_tree( @@ -197,18 +189,14 @@ def test_reforecast_single_hdate_two_steps_two_points(self): covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope_reforecast(tree) expected = [ - (0, [264.9, 265.1]), - (6, [270.1, 271.3]), + ("2025-07-14T06:00:00Z", [264.9, 265.1]), + ("2025-07-14T12:00:00Z", [270.1, 271.3]), ] assert len(covjson["coverages"]) == len(expected) - for cov, (step, vals) in zip(covjson["coverages"], expected): + for cov, (valid, vals) in zip(covjson["coverages"], expected): assert cov["domain"]["axes"] == { - "t": {"values": ["2025-07-14T06:00:00Z"]}, + "t": {"values": [valid]}, "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == vals - assert cov["mars:metadata"] == { - **EXPECTED_REFORECAST_METADATA, - "Forecast date": "2025-07-14T06:00:00Z", - "step": step, - } + assert cov["mars:metadata"] == EXPECTED_REFORECAST_METADATA diff --git a/tests/test_encoder_circle_from_polytope.py b/tests/test_encoder_circle_from_polytope.py index 4e94d5c..f6fbb5a 100644 --- a/tests/test_encoder_circle_from_polytope.py +++ b/tests/test_encoder_circle_from_polytope.py @@ -85,7 +85,7 @@ def test_reforecast_single_hdate_three_points(self): } } - assert cov["mars:metadata"] == {"Forecast date": "2025-07-14T06:00:00Z", **REFORECAST_METADATA_BASE} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_three_points(self): tree = reforecast_tree( @@ -111,4 +111,4 @@ def test_reforecast_two_hdates_three_points(self): "composite": THREE_POINTS_COMPOSITE, } assert cov["ranges"]["2t"]["values"] == range_vals - assert cov["mars:metadata"] == {"Forecast date": fc_date, **REFORECAST_METADATA_BASE} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_encoder_frame_from_polytope.py b/tests/test_encoder_frame_from_polytope.py index 743ab3e..2d7810f 100644 --- a/tests/test_encoder_frame_from_polytope.py +++ b/tests/test_encoder_frame_from_polytope.py @@ -75,10 +75,7 @@ def test_reforecast_single_hdate_two_points(self): } } - assert cov["mars:metadata"] == { - **REFORECAST_METADATA_BASE, - "Forecast date": "2025-07-14T06:00:00Z", - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_two_points(self): points = [(48.0, 11.0, [264.9]), (50.0, 12.0, [265.1])] @@ -102,7 +99,4 @@ def test_reforecast_two_hdates_two_points(self): "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == vals - assert cov["mars:metadata"] == { - **REFORECAST_METADATA_BASE, - "Forecast date": fc_date, - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_encoder_grid_from_polytope.py b/tests/test_encoder_grid_from_polytope.py index a0c362e..0eaa1d1 100644 --- a/tests/test_encoder_grid_from_polytope.py +++ b/tests/test_encoder_grid_from_polytope.py @@ -12,11 +12,14 @@ from covjsonkit.api import Covjsonkit -GRID_2X2_AXES = { - "t": {"values": [0]}, - "y": {"values": [48.0, 50.0]}, - "x": {"values": [11.0, 12.0]}, -} + +def grid_2x2_axes(t_val): + return { + "t": {"values": [t_val]}, + "y": {"values": [48.0, 50.0]}, + "x": {"values": [11.0, 12.0]}, + } + GRID_2X2_RANGES = { "2t": { @@ -90,7 +93,7 @@ def test_2x2_grid(self): assert len(covjson["coverages"]) == 1 cov = covjson["coverages"][0] - assert cov["domain"]["axes"] == GRID_2X2_AXES + assert cov["domain"]["axes"] == grid_2x2_axes("2025-01-01T00:00:00Z") assert cov["ranges"] == GRID_2X2_RANGES assert cov["mars:metadata"] == { "class": "od", @@ -126,7 +129,7 @@ def test_1x1_grid(self): cov = covjson["coverages"][0] assert cov["domain"]["axes"] == { - "t": {"values": [0]}, + "t": {"values": ["2025-01-01T00:00:00Z"]}, "y": {"values": [48.0]}, "x": {"values": [11.0]}, } @@ -160,12 +163,9 @@ def test_reforecast_single_hdate_2x2_grid(self): cov = covjson["coverages"][0] - assert cov["domain"]["axes"] == GRID_2X2_AXES + assert cov["domain"]["axes"] == grid_2x2_axes("2025-07-14T06:00:00Z") assert cov["ranges"] == GRID_2X2_RANGES - assert cov["mars:metadata"] == { - **REFORECAST_METADATA_BASE, - "Forecast date": "2025-07-14T06:00:00Z", - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_2x2_grid(self): """Two hdates each with 2x2 grid -> 2 Grid coverages.""" @@ -184,9 +184,6 @@ def test_reforecast_two_hdates_2x2_grid(self): ] assert len(covjson["coverages"]) == len(expected) for cov, fc_date in zip(covjson["coverages"], expected): - assert cov["domain"]["axes"] == GRID_2X2_AXES + assert cov["domain"]["axes"] == grid_2x2_axes(fc_date) assert cov["ranges"] == GRID_2X2_RANGES - assert cov["mars:metadata"] == { - **REFORECAST_METADATA_BASE, - "Forecast date": fc_date, - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_encoder_path_from_polytope.py b/tests/test_encoder_path_from_polytope.py index b77de66..e194bab 100644 --- a/tests/test_encoder_path_from_polytope.py +++ b/tests/test_encoder_path_from_polytope.py @@ -30,7 +30,10 @@ def test_two_points_single_step(self): "composite": { "dataType": "tuple", "coordinates": ["t", "x", "y"], - "values": [[0, 11.0, 48.0], [0, 12.0, 49.0]], + "values": [ + ["2025-01-01T00:00:00Z", 11.0, 48.0], + ["2025-01-01T00:00:00Z", 12.0, 49.0], + ], } } @@ -80,13 +83,15 @@ def test_two_points_single_step(self): class TestPathFromPolytopeReforecast: """Tests for Path (Trajectory) encoder's from_polytope_reforecast method.""" - EXPECTED_AXES = { - "composite": { - "dataType": "tuple", - "coordinates": ["t", "x", "y"], - "values": [[0, 11.0, 48.0], [0, 12.0, 50.0]], + @staticmethod + def axes_for(t_val): + return { + "composite": { + "dataType": "tuple", + "coordinates": ["t", "x", "y"], + "values": [[t_val, 11.0, 48.0], [t_val, 12.0, 50.0]], + } } - } EXPECTED_RANGES = { "2t": { @@ -113,12 +118,9 @@ def test_reforecast_single_hdate_two_points(self): cov = covjson["coverages"][0] - assert cov["domain"]["axes"] == self.EXPECTED_AXES + assert cov["domain"]["axes"] == self.axes_for("2025-07-14T06:00:00Z") assert cov["ranges"] == self.EXPECTED_RANGES - assert cov["mars:metadata"] == { - **REFORECAST_METADATA_BASE, - "Forecast date": "2025-07-14T06:00:00Z", - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_two_points(self): """Two hdates each with 2 path points -> 2 Trajectory coverages.""" @@ -137,6 +139,6 @@ def test_reforecast_two_hdates_two_points(self): ] assert len(covjson["coverages"]) == len(expected) for cov, fc_date in zip(covjson["coverages"], expected): - assert cov["domain"]["axes"] == self.EXPECTED_AXES + assert cov["domain"]["axes"] == self.axes_for(fc_date) assert cov["ranges"] == self.EXPECTED_RANGES - assert cov["mars:metadata"] == {**REFORECAST_METADATA_BASE, "Forecast date": fc_date} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_encoder_position_from_polytope.py b/tests/test_encoder_position_from_polytope.py index c21d829..6f2568d 100644 --- a/tests/test_encoder_position_from_polytope.py +++ b/tests/test_encoder_position_from_polytope.py @@ -150,7 +150,6 @@ def test_reforecast_single_hdate_two_points(self): shared_metadata = { "class": "ce", "date": "2024-03-01", - "Forecast date": "2025-07-14T06:00:00Z", "domain": "g", "expver": "4321", "levtype": "sfc", @@ -212,4 +211,4 @@ def test_reforecast_two_hdates_two_points(self): "t": {"values": t}, } assert cov["ranges"]["2t"]["values"] == vals - assert cov["mars:metadata"] == {**shared_metadata, "Forecast date": fc_date} + assert cov["mars:metadata"] == shared_metadata diff --git a/tests/test_encoder_shapefile_from_polytope.py b/tests/test_encoder_shapefile_from_polytope.py index 24289e0..5eb75bd 100644 --- a/tests/test_encoder_shapefile_from_polytope.py +++ b/tests/test_encoder_shapefile_from_polytope.py @@ -75,7 +75,7 @@ def test_reforecast_single_hdate_two_points(self): } } - assert cov["mars:metadata"] == {"Forecast date": "2025-07-14T06:00:00Z", **REFORECAST_METADATA_BASE} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_two_points(self): points = [(48.0, 11.0, [264.9]), (50.0, 12.0, [265.1])] @@ -99,4 +99,4 @@ def test_reforecast_two_hdates_two_points(self): "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == range_vals - assert cov["mars:metadata"] == {"Forecast date": fc_date, **REFORECAST_METADATA_BASE} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_encoder_vertical_profile_from_polytope.py b/tests/test_encoder_vertical_profile_from_polytope.py index a443af0..c9d4f04 100644 --- a/tests/test_encoder_vertical_profile_from_polytope.py +++ b/tests/test_encoder_vertical_profile_from_polytope.py @@ -222,7 +222,6 @@ class TestVerticalProfileFromPolytopeReforecast: "domain": "g", "expver": "4321", "levtype": "pl", - "step": 6, "stream": "efcl", "type": "sfo", "levelist": 1000, @@ -267,10 +266,7 @@ def test_reforecast_single_hdate_three_levels(self): } } - assert cov["mars:metadata"] == { - "Forecast date": "2025-07-14T06:00:00Z", - **self.EXPECTED_REFORECAST_METADATA, - } + assert cov["mars:metadata"] == self.EXPECTED_REFORECAST_METADATA def test_reforecast_two_hdates_three_levels(self): """2 hdates, 3 levels, 1 point → 2 coverages (one per hdate).""" @@ -318,7 +314,4 @@ def test_reforecast_two_hdates_three_levels(self): "values": vals, } } - assert cov["mars:metadata"] == { - "Forecast date": fc_date, - **self.EXPECTED_REFORECAST_METADATA, - } + assert cov["mars:metadata"] == self.EXPECTED_REFORECAST_METADATA diff --git a/tests/test_encoder_wkt_from_polytope.py b/tests/test_encoder_wkt_from_polytope.py index 1c4721e..f0db244 100644 --- a/tests/test_encoder_wkt_from_polytope.py +++ b/tests/test_encoder_wkt_from_polytope.py @@ -75,10 +75,7 @@ def test_reforecast_single_hdate_two_points(self): } } - assert cov["mars:metadata"] == { - "Forecast date": "2025-07-14T06:00:00Z", - **REFORECAST_METADATA_BASE, - } + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE def test_reforecast_two_hdates_two_points(self): points = [(48.0, 11.0, [264.9]), (50.0, 12.0, [265.1])] @@ -102,4 +99,4 @@ def test_reforecast_two_hdates_two_points(self): "composite": COMPOSITE_TWO_POINTS, } assert cov["ranges"]["2t"]["values"] == range_vals - assert cov["mars:metadata"] == {"Forecast date": fc_date, **REFORECAST_METADATA_BASE} + assert cov["mars:metadata"] == REFORECAST_METADATA_BASE diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index 4b55c76..ebf752b 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -364,6 +364,11 @@ def test_bounding_box_geojson_lonlat_z_order(self): gj = Covjsonkit().decode(covjson).to_geojson() assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0, 500] + @pytest.mark.xfail( + reason="Decoder round-trip does not yet support valid-time t axis (encoder-only change); " + "decoder update tracked separately.", + strict=False, + ) @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) def test_bounding_box_legacy_and_new_xarray_equivalent(self, tree_factory): new_covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(tree_factory()) From 2b4be3e6701853a84ea6e93b06bd081bb6a2e6e2 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Wed, 7 Oct 2026 13:42:49 +0200 Subject: [PATCH 07/10] Align position and vertical profile with timeseries --- covjsonkit/encoder/Position.py | 17 ++--- covjsonkit/encoder/VerticalProfile.py | 16 ++--- tests/conftest.py | 6 +- ...est_encoder_separate_datetime_pointwise.py | 64 +++++++++++++++---- 4 files changed, 65 insertions(+), 38 deletions(-) diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index d11bd98..525a8ba 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -274,15 +274,8 @@ def from_polytope_reforecast(self, result) -> dict: 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", - }, - } - ) + include_z = self._tree_has_axis(result, "levelist") + self._set_references(include_z) exclude_meta = { "latitude", @@ -327,7 +320,9 @@ def from_polytope_reforecast(self, result) -> dict: continue meta[name] = self._reforecast_stringify(rec[name]) meta["number"] = number - meta["Forecast date"] = ref.isoformat() + "Z" + if not is_reanalysis(meta, "hdate"): + # Reanalysis (class=ce, stream=efcl) exposes only the valid-time. + meta["Forecast date"] = ref.isoformat() + "Z" coverages[key] = { "lat": lat, "lon": lon, @@ -362,7 +357,7 @@ def from_polytope_reforecast(self, result) -> dict: 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) + self.add_coverage(cov["meta"], coords, val_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 21e1a9f..6d51062 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -250,15 +250,7 @@ def from_polytope_reforecast(self, result) -> dict: 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", - }, - } - ) + self._set_references() exclude_meta = { "latitude", @@ -303,8 +295,10 @@ def from_polytope_reforecast(self, result) -> dict: continue meta[name] = self._reforecast_stringify(rec[name]) meta["number"] = number - meta["step"] = step - meta["Forecast date"] = ref.isoformat() + "Z" + if not is_reanalysis(meta, "hdate"): + # Reanalysis (class=ce, stream=efcl) exposes only the valid-time. + meta["step"] = step + meta["Forecast date"] = ref.isoformat() + "Z" coverages[key] = { "lat": lat, "lon": lon, diff --git a/tests/conftest.py b/tests/conftest.py index 899cbd1..66ae065 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -215,6 +215,7 @@ def reforecast_separate_datetime_tree( step=(0,), param="167", point_factory=make_point, + stream="efcl", ): """Build a separate-datetime reforecast tree (class=ce). @@ -236,7 +237,7 @@ def reforecast_separate_datetime_tree( node("levtype", ("sfc",)), node("param", (param,)), node("step", step), - node("stream", ("efcl",)), + node("stream", (stream,)), node("time", times), node("type", ("sfo",)), ) @@ -310,6 +311,7 @@ def reforecast_separate_datetime_vertical_tree( step=(0,), param="130", point_factory=make_point, + stream="efcl", ): """Build a separate-datetime reforecast tree (class=ce) with a ``levelist`` axis. @@ -328,7 +330,7 @@ def reforecast_separate_datetime_vertical_tree( node("param", (param,)), node("levelist", levels), node("step", step), - node("stream", ("efcl",)), + node("stream", (stream,)), node("time", times), node("type", ("sfo",)), ) diff --git a/tests/test_encoder_separate_datetime_pointwise.py b/tests/test_encoder_separate_datetime_pointwise.py index a03f2dc..69b9ed0 100644 --- a/tests/test_encoder_separate_datetime_pointwise.py +++ b/tests/test_encoder_separate_datetime_pointwise.py @@ -25,7 +25,8 @@ (50.0, 12.0, [10.0, 20.0, 30.0, 40.0]), ] -# (point lat/lon, reference datetime, value) +# (point lat/lon, valid datetime, value). With a single step the valid time +# equals the reference datetime hdate + time. POSITION_EXPECTED = [ ((48.0, 11.0), "2025-07-14T00:00:00Z", 1.0), ((48.0, 11.0), "2025-07-14T12:00:00Z", 2.0), @@ -37,31 +38,51 @@ ((50.0, 12.0), "2025-07-15T12:00:00Z", 40.0), ] +POINT_SERIES_REFERENCING = [ + { + "coordinates": ["x", "y"], + "system": {"type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84"}, + }, + {"coordinates": ["t"], "system": {"type": "TemporalRS", "calendar": "Gregorian"}}, +] + 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 covjson["referencing"] == POINT_SERIES_REFERENCING 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"]] + # Surface data: no levelist axis in the tree, so no z axis. + assert set(axes) == {"x", "y", "t"} + lon = axes["x"]["values"][0] + lat = axes["y"]["values"][0] + (valid_time,) = axes["t"]["values"] (param_name,) = cov["ranges"].keys() (value,) = cov["ranges"][param_name]["values"] + # Reanalysis (efcl): only the valid-time is exposed. + assert "Forecast date" not in cov["mars:metadata"] assert "step" not in cov["mars:metadata"] assert cov["mars:metadata"]["class"] == "ce" - seen.add(((lat, lon), cov["mars:metadata"]["Forecast date"], value)) + seen.add(((lat, lon), valid_time, value)) assert seen == set(POSITION_EXPECTED) +def test_position_separate_datetime_non_efcl_keeps_forecast_date(): + tree = reforecast_separate_datetime_tree(POINTS, HDATES, TIMES, stream="enfh") + covjson = Covjsonkit().encode("CoverageCollection", "position").from_polytope_reforecast(tree) + + for cov in covjson["coverages"]: + # Single step -> the valid time equals the reference datetime. + assert cov["domain"]["axes"]["t"]["values"] == [cov["mars:metadata"]["Forecast date"]] + + 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 @@ -72,7 +93,8 @@ def test_position_separate_datetime_multi_step_collapses_into_t(): # 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"]} + # Key by the first valid time (reference + step 0). + by_ref = {c["domain"]["axes"]["t"]["values"][0]: 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"] == [ @@ -93,19 +115,33 @@ def test_verticalprofile_separate_datetime(): covjson = Covjsonkit().encode("CoverageCollection", "verticalprofile").from_polytope_reforecast(tree) assert covjson["domainType"] == "VerticalProfile" + assert [r["coordinates"] for r in covjson["referencing"]] == [["x", "y"], ["z"], ["t"]] # 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"] + by_time = {c["domain"]["axes"]["t"]["values"][0]: c for c in covjson["coverages"]} + cov = by_time["2025-07-14T00:00:00Z"] axes = cov["domain"]["axes"] - assert axes["levelist"]["values"] == [500, 850] - assert axes["t"]["values"] == ["2025-07-14T00:00:00Z"] + assert axes["x"]["values"] == [11.0] + assert axes["y"]["values"] == [48.0] + assert axes["z"]["values"] == [500, 850] (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 + assert cov["ranges"][param_name]["axisNames"] == ["z"] + # Reanalysis (efcl): only the valid-time is exposed. + assert "Forecast date" not in cov["mars:metadata"] + assert "step" not in cov["mars:metadata"] + + +def test_verticalprofile_separate_datetime_non_efcl_keeps_forecast_date(): + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + tree = reforecast_separate_datetime_vertical_tree(points, HDATES, TIMES, (500, 850), stream="enfh") + covjson = Covjsonkit().encode("CoverageCollection", "verticalprofile").from_polytope_reforecast(tree) + + for cov in covjson["coverages"]: + assert cov["domain"]["axes"]["t"]["values"] == [cov["mars:metadata"]["Forecast date"]] + assert cov["mars:metadata"]["step"] == 0 def test_pointwise_reforecast_backward_compat_delegates(): From 4d6017a5918581b0b6eae18e58010ee29b031751 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 8 Oct 2026 09:03:15 +0200 Subject: [PATCH 08/10] Update encoders to calculate valid time correctly --- covjsonkit/encoder/BoundingBox.py | 17 +- covjsonkit/encoder/Circle.py | 17 +- covjsonkit/encoder/Frame.py | 17 +- covjsonkit/encoder/Grid.py | 53 +++- covjsonkit/encoder/Path.py | 45 ++- covjsonkit/encoder/Position.py | 115 ++------ covjsonkit/encoder/Shapefile.py | 17 +- covjsonkit/encoder/TimeSeries.py | 92 ++---- covjsonkit/encoder/VerticalProfile.py | 126 ++------ covjsonkit/encoder/Wkt.py | 24 +- covjsonkit/encoder/encoder.py | 398 +++++++++++++++----------- 11 files changed, 452 insertions(+), 469 deletions(-) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index 4381c9b..35fbc48 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -1,7 +1,13 @@ import logging import time -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class BoundingBox(Encoder): @@ -41,8 +47,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -118,7 +125,9 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict, include_z) + cov_coords = dict(coords) + cov_coords["t"] = [valid_time(datetime, step)] + self.add_coverage(mars_metadata, cov_coords, dv_dict, include_z) # Return the generated CoverageJSON return self.covjson diff --git a/covjsonkit/encoder/Circle.py b/covjsonkit/encoder/Circle.py index eee67f2..add05ef 100644 --- a/covjsonkit/encoder/Circle.py +++ b/covjsonkit/encoder/Circle.py @@ -1,6 +1,12 @@ import logging -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class Circle(Encoder): @@ -36,8 +42,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -109,7 +116,9 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict, include_z) + cov_coords = dict(coords) + cov_coords["t"] = [valid_time(datetime, step)] + self.add_coverage(mars_metadata, cov_coords, dv_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/Frame.py b/covjsonkit/encoder/Frame.py index 67628ac..4e2ecfd 100644 --- a/covjsonkit/encoder/Frame.py +++ b/covjsonkit/encoder/Frame.py @@ -1,6 +1,12 @@ import logging -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class Frame(Encoder): @@ -36,8 +42,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -109,7 +116,9 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict, include_z) + cov_coords = dict(coords) + cov_coords["t"] = [valid_time(datetime, step)] + self.add_coverage(mars_metadata, cov_coords, dv_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index b368434..84a6168 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -123,7 +123,7 @@ def from_xarray(self, dataset): arr = dataset[dv].sel(datetimes=datetime, number=num, steps=step).values dv_dict[dv] = arr.ravel().tolist() - coords = {"t": [str(step)], "y": y_vals, "x": x_vals} + coords = {"t": [valid_time(datetime, step)], "y": y_vals, "x": x_vals} if include_z: coords["z"] = z_vals self.add_coverage(mars_metadata, coords, dv_dict, include_z) @@ -227,12 +227,61 @@ def from_polytope(self, result, date_key: str = "date") -> dict: mm.pop("Forecast date", None) mm.pop("step", None) else: - mm["step"] = normalize_step_value(step) + if len(fields["step"]) == 1: + mm["step"] = normalize_step_value(fields["step"][0]) + else: + # ``t`` carries the valid time of every step. + mm.pop("step", None) mm["Forecast date"] = date self.add_coverage(mm, coordinates[date], val_dict, include_z) return self.covjson + def from_polytope_reforecast(self, result) -> dict: + """Encode reforecast/reanalysis data into a Grid collection. + + Legacy merged trees (no independent ``time`` axis) are delegated to + :meth:`from_polytope` with ``date_key="hdate"``. For separate-datetime + (``class=ce``) trees one coverage is produced per + ``(reference = hdate + time, number)`` whose ``t`` axis holds the valid + time of every step, with values ordered ``[t, (z,) y, x]``. + """ + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + include_z = self._tree_has_axis(result, "levelist") + self._set_references(include_z) + + coverages, param_order = self._reforecast_coverages(result, by_step=False) + for para in param_order: + self.add_parameter(para) + + for cov in coverages: + entries = cov["entries"] + times = sorted({valid for valid, _lat, _lon, _level in entries}) + levels = list(dict.fromkeys(level for _valid, _lat, _lon, level in entries)) + lats = list(dict.fromkeys(lat for _valid, lat, _lon, _level in entries)) + lons = list(dict.fromkeys(lon for _valid, _lat, lon, _level in entries)) + + coords = {"t": times, "y": lats, "x": lons} + if include_z: + coords["z"] = levels + self.shp = [len(times), len(levels), len(lats), len(lons)] + else: + self.shp = [len(times), len(lats), len(lons)] + + val_dict = {para: [] for para in param_order} + for valid in times: + for level in levels: + for lat in lats: + for lon in lons: + values = entries.get((valid, lat, lon, level), {}) + for para in param_order: + val_dict[para].append(values.get(para)) + self.add_coverage(dict(cov["meta"]), coords, val_dict, include_z) + + return self.covjson + def from_polytope_month(self, result): coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Path.py b/covjsonkit/encoder/Path.py index a42c9bc..b90b316 100644 --- a/covjsonkit/encoder/Path.py +++ b/covjsonkit/encoder/Path.py @@ -1,6 +1,12 @@ import logging -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class Path(Encoder): @@ -40,8 +46,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -239,6 +246,38 @@ 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 Trajectory collection. + + Legacy merged trees (no independent ``time`` axis) are delegated to + :meth:`from_polytope` with ``date_key="hdate"``. For separate-datetime + (``class=ce``) trees one coverage is produced per + ``(reference = hdate + time, number)`` whose composite holds + ``[t, x, y(, z)]`` for every point, ``t`` being its valid time. + """ + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + include_z = self._tree_has_axis(result, "levelist") + self._set_references(include_z) + + coverages, param_order = self._reforecast_coverages(result, by_step=False) + for para in param_order: + self.add_parameter(para) + + for cov in coverages: + # Stable sort by valid time keeps the path's point order within a step. + entries = sorted(cov["entries"].items(), key=lambda item: item[0][0]) + composite = [ + [valid, lon, lat, level] if include_z else [valid, lon, lat] for (valid, lat, lon, level), _ in entries + ] + val_dict = {para: [values.get(para) for _, values in entries] for para in param_order} + meta = dict(cov["meta"]) + meta.pop("levelist", None) + self.add_coverage(meta, {"composite": composite}, val_dict, include_z) + + return self.covjson + def from_polytope_month(self, result): coords = {} mars_metadata = {} diff --git a/covjsonkit/encoder/Position.py b/covjsonkit/encoder/Position.py index 525a8ba..6505735 100644 --- a/covjsonkit/encoder/Position.py +++ b/covjsonkit/encoder/Position.py @@ -1,10 +1,7 @@ import logging import time -from datetime import datetime, timedelta -import pandas as pd - -from .encoder import Encoder, is_reanalysis, normalize_step_value +from .encoder import Encoder, is_reanalysis, range_shape, valid_time class Position(Encoder): @@ -42,8 +39,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["t"] + shape, axis_names = range_shape(values[parameter], "t") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[ parameter ] # [values[parameter][val][0] for val in values[parameter].keys()] @@ -190,19 +188,7 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for num in fields["number"]: for para in fields["param"]: for step in fields["step"]: - date_format = "%Y%m%dT%H%M%S" - new_date = pd.Timestamp(date).strftime(date_format) - start_time = datetime.strptime(new_date, date_format) - # add current date to list by converting it to iso format - if isinstance(step, timedelta): - stamp = start_time + step - else: - try: - int(step) - except ValueError: - step = step[0] - stamp = start_time + timedelta(hours=int(step)) - coordinates[date][i]["t"].append(stamp.isoformat() + "Z") + coordinates[date][i]["t"].append(valid_time(date, step)) break break break @@ -277,87 +263,19 @@ def from_polytope_reforecast(self, result) -> dict: include_z = self._tree_has_axis(result, "levelist") self._set_references(include_z) - 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 - if not is_reanalysis(meta, "hdate"): - # Reanalysis (class=ce, stream=efcl) exposes only the valid-time. - 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.") - + coverages, param_order = self._reforecast_coverages(result, by_step=False) 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, include_z) + for cov in coverages: + series_by_point = {} + for (valid, lat, lon, level), values in cov["entries"].items(): + series_by_point.setdefault((lat, lon, level), {})[valid] = values + for (lat, lon, level), series in series_by_point.items(): + times = sorted(series) + coords = {"latitude": [lat], "longitude": [lon], "levelist": [level], "t": times} + val_dict = {para: [series[t].get(para) for t in times] for para in param_order} + self.add_coverage(dict(cov["meta"]), coords, val_dict, include_z) return self.covjson @@ -506,9 +424,8 @@ def from_polytope_step(self, result): for para in fields["param"]: for date in fields["dates"]: for times in fields["times"]: - datetime = pd.Timestamp(date) + times coordinates[fields["dates"][0]][(i * len(fields["levels"]) + j)]["t"].append( - str(datetime).split("+")[0] + "Z" + self._step_path_valid_time(date, times, fields) ) break break diff --git a/covjsonkit/encoder/Shapefile.py b/covjsonkit/encoder/Shapefile.py index 6be6dbe..d06b1b1 100644 --- a/covjsonkit/encoder/Shapefile.py +++ b/covjsonkit/encoder/Shapefile.py @@ -1,6 +1,12 @@ import logging -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class Shapefile(Encoder): @@ -36,8 +42,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -109,7 +116,9 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict, include_z) + cov_coords = dict(coords) + cov_coords["t"] = [valid_time(datetime, step)] + self.add_coverage(mars_metadata, cov_coords, dv_dict, include_z) return self.covjson diff --git a/covjsonkit/encoder/TimeSeries.py b/covjsonkit/encoder/TimeSeries.py index ada6471..0ec8553 100644 --- a/covjsonkit/encoder/TimeSeries.py +++ b/covjsonkit/encoder/TimeSeries.py @@ -4,7 +4,15 @@ import pandas as pd -from .encoder import Encoder +from .encoder import ( + Encoder, + iso_utc, + range_shape, + reference_datetime, + set_forecast_date, + valid_datetime, + valid_time, +) class TimeSeries(Encoder): @@ -42,14 +50,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - # A single-valued t axis is emitted as a 0D array: the spec allows - # single-valued axes to be omitted and covjson-validator rejects ["t"]. - if len(values[parameter]) == 1: - coverage["ranges"][param]["shape"] = [] - coverage["ranges"][param]["axisNames"] = [] - else: - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["t"] + shape, axis_names = range_shape(values[parameter], "t") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[ parameter ] # [values[parameter][val][0] for val in values[parameter].keys()] @@ -80,33 +83,6 @@ def _set_references(self, include_z): ) self.covjson["referencing"] = refs - @staticmethod - 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: - int(step) - except ValueError: - step = step[0] - return start_time + timedelta(hours=int(step)) - def _collapse_reanalysis(self, fields, coords, mars_metadata, range_dict): """Collapse all hdates for each point into a single PointSeries coverage. @@ -130,12 +106,12 @@ def _collapse_reanalysis(self, fields, coords, mars_metadata, range_dict): stamp_order = [] for date in fields["dates"]: for step in fields["step"]: - stamp = self._hdate_step_timestamp(date, step, time_offset) + stamp = valid_datetime(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]) - t_values = [stamp.isoformat() + "Z" for stamp, _, _ in stamp_order] + t_values = [iso_utc(stamp) for stamp, _, _ in stamp_order] for i in range(points): lat = coords[first_date]["composite"][i][0] @@ -297,19 +273,7 @@ def from_polytope(self, result, date_key: str = "date", reforecast: bool = False for num in fields["number"]: for para in fields["param"]: for step in fields["step"]: - date_format = "%Y%m%dT%H%M%S" - new_date = pd.Timestamp(date).strftime(date_format) - start_time = datetime.strptime(new_date, date_format) - # add current date to list by converting it to iso format - if isinstance(step, timedelta): - stamp = start_time + step - else: - try: - int(step) - except ValueError: - step = step[0] - stamp = start_time + timedelta(hours=int(step)) - coordinates[date][i]["t"].append(stamp.isoformat() + "Z") + coordinates[date][i]["t"].append(valid_time(date, step)) break break break @@ -448,9 +412,7 @@ def emit(full_path, flat_result): 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) + stamp = valid_datetime(hdate, step, time_off) # Stream determines the coverage grouping for class=ce: # * efcl (reforecast/reanalysis): collapse all hdates into a @@ -468,7 +430,7 @@ def emit(full_path, flat_result): key = (lat, lon, level, number) elif forecast: # Reference datetime of the forecast run = date + time. - reference = self._hdate_step_timestamp(hdate, 0, time_off) + reference = reference_datetime(hdate, time_off) key = (lat, lon, level, number, reference) else: key = (lat, lon, level, number, stringify(hdate)) @@ -480,13 +442,9 @@ def emit(full_path, flat_result): 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" + # efas: Forecast date is the run (date + time); efcl: none; + # otherwise one coverage per hdate. + set_forecast_date(meta, reference if forecast else hdate) coverages[key] = { "lat": lat, "lon": lon, @@ -545,7 +503,7 @@ def recurse(node, path): 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] + t_values = [iso_utc(reference[i][0]) for i in order] val_dict = {para: [params[para][i][1] for i in order] for para in paras} coord_entry = { "latitude": [cov["lat"]], @@ -717,14 +675,8 @@ def from_polytope_step(self, result): for para in fields["param"]: for date in fields["dates"]: for times in fields["times"]: - # date_format = "%Y%m%dT%H%M%S" - # new_date = pd.Timestamp(date).strftime(date_format) - # start_time = datetime.strptime(new_date, date_format) - # add current date to list by converting it to iso format - # stamp = start_time + timedelta(hours=int(step)) - datetime = pd.Timestamp(date) + times coordinates[fields["dates"][0]][(i * len(fields["levels"]) + j)]["t"].append( - str(datetime).split("+")[0] + "Z" + self._step_path_valid_time(date, times, fields) ) break break diff --git a/covjsonkit/encoder/VerticalProfile.py b/covjsonkit/encoder/VerticalProfile.py index 6d51062..fdb3cde 100644 --- a/covjsonkit/encoder/VerticalProfile.py +++ b/covjsonkit/encoder/VerticalProfile.py @@ -1,10 +1,13 @@ import logging import time -from datetime import datetime, timedelta -import pandas as pd - -from .encoder import Encoder, is_reanalysis, normalize_step_value +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class VerticalProfile(Encoder): @@ -41,8 +44,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["z"] + shape, axis_names = range_shape(values[parameter], "z") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -114,7 +118,7 @@ def from_xarray(self, datasets): "latitude": [float(x) for x in dataset["latitude"].values], "longitude": [float(x) for x in dataset["longitude"].values], "levelist": [float(x) for x in dataset["levelist"].values], - "t": [str(x) for x in dataset["datetime"].values], + "t": [valid_time(x, step) for x in dataset["datetime"].values], }, dv_dict, ) @@ -164,24 +168,12 @@ def from_polytope(self, result, date_key: str = "date") -> dict: for i, point in enumerate(range(points)): coordinates[date][i] = [] for step in fields["step"]: - date_format = "%Y%m%dT%H%M%S" - new_date = pd.Timestamp(date).strftime(date_format) - start_time = datetime.strptime(new_date, date_format) - # add current date to list by converting it to iso format - if isinstance(step, timedelta): - stamp = start_time + step - else: - try: - int(step) - except ValueError: - step = step[0] - stamp = start_time + timedelta(hours=int(step)) coordinates[date][i].append( { "latitude": [coords[date]["composite"][i][0]], "longitude": [coords[date]["composite"][i][1]], "levelist": list(levels), - "t": [stamp.isoformat() + "Z"], + "t": [valid_time(date, step)], } ) @@ -245,97 +237,27 @@ def from_polytope_reforecast(self, result) -> dict: (``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``. + its ``z`` 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._set_references() - 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 - if not is_reanalysis(meta, "hdate"): - # Reanalysis (class=ce, stream=efcl) exposes only the valid-time. - 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.") - + coverages, param_order = self._reforecast_coverages(result, by_step=True) 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) + for cov in coverages: + profile_by_point = {} + for (valid, lat, lon, level), values in cov["entries"].items(): + profile = profile_by_point.setdefault((lat, lon), {"t": valid, "levels": {}}) + profile["levels"][level] = values + for (lat, lon), profile in profile_by_point.items(): + levels = sorted(profile["levels"]) + coords = {"latitude": [lat], "longitude": [lon], "levelist": levels, "t": [profile["t"]]} + val_dict = {para: [profile["levels"][lev].get(para) for lev in levels] for para in param_order} + self.add_coverage(dict(cov["meta"]), coords, val_dict) return self.covjson diff --git a/covjsonkit/encoder/Wkt.py b/covjsonkit/encoder/Wkt.py index cd552cb..645647b 100644 --- a/covjsonkit/encoder/Wkt.py +++ b/covjsonkit/encoder/Wkt.py @@ -2,9 +2,13 @@ import logging import time -import pandas as pd - -from .encoder import Encoder, is_reanalysis, normalize_step_value, valid_time +from .encoder import ( + Encoder, + is_reanalysis, + normalize_step_value, + range_shape, + valid_time, +) class Wkt(Encoder): @@ -40,8 +44,9 @@ def add_range(self, coverage, values): coverage["ranges"][param] = {} coverage["ranges"][param]["type"] = "NdArray" coverage["ranges"][param]["dataType"] = "float" - coverage["ranges"][param]["shape"] = [len(values[parameter])] - coverage["ranges"][param]["axisNames"] = ["composite"] + shape, axis_names = range_shape(values[parameter], "composite") + coverage["ranges"][param]["shape"] = shape + coverage["ranges"][param]["axisNames"] = axis_names coverage["ranges"][param]["values"] = values[parameter] def add_mars_metadata(self, coverage, metadata): @@ -113,7 +118,9 @@ def from_xarray(self, dataset): for dv in dataset.data_vars: dv_dict[dv] = dataset[dv].sel(number=num, steps=step, datetimes=datetime).values.tolist() - self.add_coverage(mars_metadata, coords, dv_dict, include_z) + cov_coords = dict(coords) + cov_coords["t"] = [valid_time(datetime, step)] + self.add_coverage(mars_metadata, cov_coords, dv_dict, include_z) return self.covjson @@ -278,9 +285,8 @@ def from_polytope_step(self, result): val_dict[para].extend(vals) mm = mars_metadata.copy() mm["number"] = num - # mm["Forecast date"] = date - datetime = pd.Timestamp(date) + t - self.add_coverage(mm, coordinates[str(datetime).split("+")[0] + "Z"], val_dict, include_z) + valid = self._step_path_valid_time(date, t, fields) + self.add_coverage(mm, coordinates[valid], val_dict, include_z) end = time.time() delta = end - start diff --git a/covjsonkit/encoder/encoder.py b/covjsonkit/encoder/encoder.py index cd3e4c7..cf75764 100644 --- a/covjsonkit/encoder/encoder.py +++ b/covjsonkit/encoder/encoder.py @@ -1,5 +1,6 @@ from __future__ import annotations +import re from abc import ABC, abstractmethod from datetime import datetime, timedelta from typing import Any @@ -32,6 +33,17 @@ def is_merged_node(node) -> bool: return hasattr(node, "axes") and getattr(node, "axes", None) is not None +def range_shape(values, axis_name): + """Return ``(shape, axisNames)`` for a one-dimensional NdArray range. + + A single value is emitted as a 0D array (``[]``, ``[]``): the spec allows + single-valued axes to be omitted and covjson-validator rejects them. + """ + if len(values) == 1: + return [], [] + return [len(values)], [axis_name] + + def timedelta_to_step_string(td: timedelta) -> str: """ Convert a timedelta object to a step string in the format 'XhYm'. @@ -64,38 +76,110 @@ def timedelta_to_step_string(td: timedelta) -> str: def parse_step_string(step_str: str) -> float: """ - Parse a step string in format 'XhYm' and return total hours as float. + Parse a step value and return total hours as float. - Args: - step_str: Step string in format 'Xh', 'Ym', or 'XhYm' - - Returns: - Total hours as a float value + Accepts every format understood by :func:`step_to_timedelta`. Examples: parse_step_string("13h") -> 13.0 parse_step_string("13h30m") -> 13.5 parse_step_string("30m") -> 0.5 - parse_step_string("0h") -> 0.0 + parse_step_string("6") -> 6.0 + parse_step_string("0-6") -> 6.0 + """ + return step_to_timedelta(step_str).total_seconds() / 3600.0 + + +_STEP_UNITS_RE = re.compile(r"(?:(\d+(?:\.\d+)?)h)?(?:(\d+(?:\.\d+)?)m)?") + + +def step_to_timedelta(step) -> timedelta: + """Convert a MARS/polytope ``step`` value to a :class:`datetime.timedelta`. + + Accepts timedelta-likes, hour counts (int, float or numeric string), unit + strings (``"6h"``, ``"30m"``, ``"13h30m"``) and single-element sequences. + A step range ``"a-b"`` (accumulations, means) is valid at the end of its + period, so ``b`` is used. """ - if isinstance(step_str, (int, float)): - return float(step_str) + if step is None: + return timedelta(0) + if isinstance(step, (list, tuple, np.ndarray)): + if len(step) != 1: + raise ValueError(f"Expected a single step value, got {step!r}") + return step_to_timedelta(step[0]) + if isinstance(step, (timedelta, pd.Timedelta, np.timedelta64)): + return pd.Timedelta(step).to_pytimedelta() + if isinstance(step, (int, float, np.integer, np.floating)): + return timedelta(hours=float(step)) + + text = str(step).strip() + if "-" in text.lstrip("-"): + text = text.rsplit("-", 1)[1] + try: + return timedelta(hours=float(text)) + except ValueError: + pass + match = _STEP_UNITS_RE.fullmatch(text) + if not match or not any(match.groups()): + raise ValueError(f"Unrecognised step value: {step!r}") + hours, minutes = match.groups() + return timedelta(hours=float(hours or 0), minutes=float(minutes or 0)) - step_str = str(step_str) - hours = 0.0 - minutes = 0.0 - # Parse hours - if "h" in step_str: - parts = step_str.split("h") - hours = float(parts[0]) - step_str = parts[1] if len(parts) > 1 else "" +def time_to_timedelta(value) -> timedelta: + """Convert a MARS ``time`` (time-of-day) value to a :class:`datetime.timedelta`. - # Parse minutes - if "m" in step_str and step_str: - minutes = float(step_str.replace("m", "")) + Accepts timedelta-likes, ``"HH:MM[:SS]"`` strings and ``HHMM`` numbers or + strings (``600``, ``"0600"``, ``"1200"``). Numbers below 100 are hours. + """ + if value is None: + return timedelta(0) + if isinstance(value, (list, tuple, np.ndarray)): + if len(value) != 1: + raise ValueError(f"Expected a single time value, got {value!r}") + return time_to_timedelta(value[0]) + if isinstance(value, (timedelta, pd.Timedelta, np.timedelta64)): + return pd.Timedelta(value).to_pytimedelta() + if isinstance(value, str) and ":" in value: + return pd.to_timedelta(value).to_pytimedelta() + try: + number = int(float(str(value))) + except ValueError: + return pd.to_timedelta(value).to_pytimedelta() + if number < 100: + return timedelta(hours=number) + return timedelta(hours=number // 100, minutes=number % 100) + + +def to_naive_utc(value) -> pd.Timestamp: + """Return ``value`` as a timezone-naive UTC :class:`pandas.Timestamp`. + + Accepts anything ``pandas.Timestamp`` does, including ISO strings with a + trailing ``"Z"`` and timezone-aware values (converted to UTC). + """ + if isinstance(value, str): + value = value.strip() + if value.endswith("Z"): + value = value[:-1] + stamp = pd.Timestamp(value) + if stamp.tzinfo is not None: + stamp = stamp.tz_convert("UTC").tz_localize(None) + return stamp - return hours + (minutes / 60.0) + +def iso_utc(value) -> str: + """Format a datetime-like as an ISO-8601 UTC string with a trailing ``"Z"``.""" + return to_naive_utc(value).isoformat() + "Z" + + +def reference_datetime(date, time=None) -> pd.Timestamp: + """Reference datetime of a forecast/reforecast: ``date`` (or ``hdate``) + ``time``.""" + return to_naive_utc(date) + time_to_timedelta(time) + + +def valid_datetime(date, step=0, time=None) -> pd.Timestamp: + """Valid datetime: ``date + time + step`` (see :func:`step_to_timedelta`).""" + return reference_datetime(date, time) + step_to_timedelta(step) def sort_step_values(steps: list) -> list: @@ -196,31 +280,14 @@ def normalize_step_value(step): return str(step) -def valid_time(date, step): - """Return the ISO-8601 valid-time string for a forecast ``date`` and ``step``. - - The valid time is ``date + step``. ``date`` may be any value accepted by - ``pandas.Timestamp`` (e.g. a MARS ``"YYYYMMDDTHHMMSS"`` string, optionally - suffixed with ``"Z"``). ``step`` may be a ``timedelta``/``numpy.timedelta64`` - or an hour count (int, float, numeric string, or a single-element list). +def valid_time(date, step=0, time=None) -> str: + """Return the ISO-8601 valid-time string ``date + time + step`` with a trailing ``"Z"``. - The result is an ISO-8601 string with a trailing ``"Z"``, matching the - valid-time strings produced by the PointSeries encoders. + ``date`` may be any value accepted by :func:`to_naive_utc`, ``time`` any + value accepted by :func:`time_to_timedelta` and ``step`` any value accepted + by :func:`step_to_timedelta`. """ - date_format = "%Y%m%dT%H%M%S" - new_date = pd.Timestamp(str(date).rstrip("Z")).strftime(date_format) - start_time = pd.Timestamp(new_date).to_pydatetime() - if isinstance(step, (timedelta, pd.Timedelta, np.timedelta64)): - stamp = start_time + pd.Timedelta(step).to_pytimedelta() - else: - if isinstance(step, (list, tuple)): - step = step[0] - try: - hours = int(step) - except (TypeError, ValueError): - hours = int(str(step).split("h")[0]) - stamp = start_time + timedelta(hours=hours) - return stamp.isoformat() + "Z" + return iso_utc(valid_datetime(date, step, time)) def is_reanalysis(mars_metadata: dict, date_key: str = "date") -> bool: @@ -233,6 +300,25 @@ def is_reanalysis(mars_metadata: dict, date_key: str = "date") -> bool: return date_key == "hdate" and mars_metadata.get("class") == "ce" and mars_metadata.get("stream") == "efcl" +def set_forecast_date(meta: dict, reference) -> dict: + """Apply the ``"Forecast date"`` rules to ``meta`` in place and return it. + + * ``class=ce, stream=efcl`` (reanalysis): no forecast date; ``t`` carries + the valid time. + * ``class=ce, stream=efas`` (forecast): ``date + time`` is the run + reference, so the raw ``date`` is dropped in favour of ``"Forecast date"``. + * otherwise ``"Forecast date"`` is ``reference``. + """ + is_ce = meta.get("class") == "ce" + if is_ce and meta.get("stream") == "efcl": + meta.pop("Forecast date", None) + return meta + if is_ce and meta.get("stream") == "efas": + meta.pop("date", None) + meta["Forecast date"] = iso_utc(reference) + return meta + + class Encoder(ABC): def __init__(self, type, domaintype): """ @@ -415,7 +501,7 @@ def handle_specific_axes(child): return child.values if child.axis.name == "param": return child.values - if child.axis.name in [date_key, "time"]: + 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: @@ -423,12 +509,40 @@ def handle_specific_axes(child): coords[date]["composite"] = [] coords[date]["t"] = [date] return dates + if child.axis.name == "time": + fold_time_into_date(child.values) + return None if child.axis.name == "number": return child.values if child.axis.name == "step": return child.values return None + def fold_time_into_date(times): + """Shift the current date key by a single time-of-day offset. + + Independent ``date``/``time`` axes with several times are encoded by + ``from_polytope_step``/``from_polytope_reforecast``; here only a + single time can be folded into the date dimension. + """ + if len(times) != 1: + raise ValueError( + "A multi-valued 'time' axis is not supported by from_polytope; " + "use from_polytope_step or from_polytope_reforecast." + ) + base_dates = fields.setdefault("_base_dates", {}) + current = fields["dates"][-1] + base = base_dates.get(current, current) + shifted = iso_utc(reference_datetime(base, times[0])) + if shifted == current: + return + base_dates[shifted] = base + entry = coords.pop(current) + entry["t"] = [shifted] + coords.setdefault(shifted, entry) + fields["dates"][-1] = shifted + mars_metadata["Forecast date"] = shifted + 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) @@ -497,7 +611,7 @@ def emit_leaf(lat, lon_values, result): fields["l"].extend(result) elif child.axis.name == "param": fields["param"] = result - elif child.axis.name in [date_key, "time"]: + elif child.axis.name == date_key: fields["dates"].extend(result) elif child.axis.name == "number": fields["number"] = result @@ -562,7 +676,7 @@ def handle_specific_axes(child): # 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"] = self._reforecast_timedelta(child.values[0]) + fields["time_offset"] = time_to_timedelta(child.values[0]) return None if child.axis.name == "number": return child.values @@ -684,10 +798,7 @@ def handle_specific_axes_step(child): for date in fields["dates"]: coords[date] = {} coords[date]["composite"] = [] - coords[date]["t"] = [] - for time in child.values: - datetime = pd.Timestamp(date) + time - coords[date]["t"].append(str(datetime).split("+")[0] + "Z") + coords[date]["t"] = [self._step_path_valid_time(date, time, fields) for time in child.values] return child.values return None @@ -979,35 +1090,16 @@ def _reforecast_stringify(value): 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) + def _step_path_valid_time(self, date, time, fields) -> str: + """Valid time for the date/time ``*_step`` paths: ``date + time + step``. - @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 + Their values are indexed by date and time only, so a step is applied + when it is single-valued. + """ + steps = fields.get("step") + if isinstance(steps, (list, tuple, np.ndarray)) and len(steps) == 1: + return valid_time(date, steps[0], time) + return valid_time(date, 0, time) @staticmethod def _tree_has_axis(tree, axis_name) -> bool: @@ -1063,59 +1155,24 @@ def recurse(node, path): recurse(result, []) return records - def from_polytope_reforecast(self, result) -> dict: - """Encode reforecast/reanalysis data that uses ``"hdate"`` as the time axis. + def _reforecast_coverages(self, result, by_step=True): + """Group separate-datetime (``class=ce``) records into coverages. - Two representations are supported: + Coverages are keyed by ``(reference, number)`` -- plus ``step`` when + ``by_step`` -- where the reference datetime is ``hdate`` (or ``date``) + ``+ time``. Each coverage maps ``(valid, lat, lon, level)`` entries, in + record order, to ``{param: value}``. Coverages are returned in + reference -> step -> number order with ``"Forecast date"`` applied by + :func:`set_forecast_date`. - * **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. + Returns ``(coverages, param_order)``. """ - 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) + exclude_meta = {"latitude", "longitude", "hdate", "time", "step", "param", "number", "levelist"} 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) @@ -1128,65 +1185,70 @@ def to_timedelta(value): 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) + ref = reference_datetime(rec.get("hdate", rec.get("date")), rec.get("time")) + valid = iso_utc(ref + step_to_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) + key = (ref, step_to_timedelta(step), number) if by_step else (ref, 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[name] = self._reforecast_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 by_step: + meta["step"] = step + set_forecast_date(meta, ref) + if "Forecast date" not in meta: + # Reanalysis (efcl): the valid time on ``t`` replaces date and step. + meta.pop("step", None) + coverages[key] = {"meta": meta, "entries": {}} + if para not in param_order: param_order.append(para) - cov["values"].setdefault(para, {})[point] = float(value) + entry = (valid, float(rec["latitude"]), float(rec["longitude"]), level) + coverages[key]["entries"].setdefault(entry, {})[para] = float(rec["__value__"]) if not coverages: raise ValueError("No data was returned.") + ordered = [coverages[key] for key in sorted(coverages, key=lambda k: k)] + return ordered, param_order + + def from_polytope_reforecast(self, result) -> dict: + """Encode reforecast/reanalysis data that uses ``"hdate"`` as the time axis. + + 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 MultiPoint coverage is produced per + ``(reference-datetime, step, number)`` combination, holding all spatial + points in its ``composite`` axis and the valid time on ``t``. + + Encoders whose domain is not MultiPoint (Grid, Path, PointSeries, + VerticalProfile) override this method. + """ + if not self._tree_has_axis(result, "time"): + return self.from_polytope(result, date_key="hdate") + + include_z = self._tree_has_axis(result, "levelist") + self._set_references(include_z) + + coverages, param_order = self._reforecast_coverages(result, by_step=True) 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) + for cov in coverages: + entries = cov["entries"] + valid = next(iter(entries))[0] + composite = [] + for _valid, lat, lon, level in entries: + composite.append([lon, lat, level] if include_z else [lon, lat]) + val_dict = {para: [values.get(para) for values in entries.values()] for para in param_order} + self.add_coverage(cov["meta"], {"composite": composite, "t": [valid]}, val_dict, include_z) return self.covjson From 54a2c9fb975b61b4413968380a7a9582a24f92c0 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 8 Oct 2026 10:44:00 +0200 Subject: [PATCH 09/10] Align decoders with new outputs --- covjsonkit/decoder/BoundingBox.py | 28 +++++++++++++-------------- covjsonkit/decoder/Circle.py | 28 +++++++++++++-------------- covjsonkit/decoder/Grid.py | 9 ++++----- covjsonkit/decoder/Path.py | 26 +++++++++++-------------- covjsonkit/decoder/Position.py | 4 ++++ covjsonkit/decoder/VerticalProfile.py | 15 +++++--------- covjsonkit/decoder/decoder.py | 19 ++++++++++++++++++ 7 files changed, 69 insertions(+), 60 deletions(-) diff --git a/covjsonkit/decoder/BoundingBox.py b/covjsonkit/decoder/BoundingBox.py index 0e91430..e1e59f4 100644 --- a/covjsonkit/decoder/BoundingBox.py +++ b/covjsonkit/decoder/BoundingBox.py @@ -183,21 +183,16 @@ def to_xarray(self): if "step" not in coverage["mars:metadata"]: coverage["mars:metadata"]["step"] = 0 steps.append(coverage["mars:metadata"]["step"]) - datetimes.append(coverage["domain"]["axes"]["t"]["values"][0]) + # ``t`` is the valid time (date + step); the datetimes dimension is + # the forecast reference so that datetimes x steps stays dense. + datetime = self.forecast_reference(coverage) + datetimes.append(datetime) + number = coverage["mars:metadata"]["number"] + step = coverage["mars:metadata"]["step"] for parameter in self.parameters: - # values[parameter].append(coverage["ranges"][parameter]["values"]) - if coverage["domain"]["axes"]["t"]["values"][0] not in values[parameter]: - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]] = {} - if ( - coverage["mars:metadata"]["number"] - not in values[parameter][coverage["domain"]["axes"]["t"]["values"][0]] - ): - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]][ - coverage["mars:metadata"]["number"] - ] = {} - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]][coverage["mars:metadata"]["number"]][ - coverage["mars:metadata"]["step"] - ] = coverage["ranges"][parameter]["values"] + values[parameter].setdefault(datetime, {}).setdefault(number, {})[step] = coverage["ranges"][parameter][ + "values" + ] datetimes = sorted(list(set(datetimes))) numbers = sorted(list(set(numbers))) @@ -211,7 +206,10 @@ def to_xarray(self): for j, number in enumerate(numbers): new_values[parameter][i].append([]) for k, step in enumerate(steps): - new_values[parameter][i][j].append(values[parameter][datetime][number][step]) + cell = values[parameter].get(datetime, {}).get(number, {}).get(step) + if cell is None: + cell = [None] * len(longitude) + new_values[parameter][i][j].append(cell) for parameter in self.parameters: dataarray = xr.DataArray(new_values[parameter], dims=dims) diff --git a/covjsonkit/decoder/Circle.py b/covjsonkit/decoder/Circle.py index adddfca..34863cf 100644 --- a/covjsonkit/decoder/Circle.py +++ b/covjsonkit/decoder/Circle.py @@ -124,21 +124,16 @@ def to_xarray(self): if "step" not in coverage["mars:metadata"]: coverage["mars:metadata"]["step"] = 0 steps.append(coverage["mars:metadata"]["step"]) - datetimes.append(coverage["domain"]["axes"]["t"]["values"][0]) + # ``t`` is the valid time (date + step); the datetimes dimension is + # the forecast reference so that datetimes x steps stays dense. + datetime = self.forecast_reference(coverage) + datetimes.append(datetime) + number = coverage["mars:metadata"]["number"] + step = coverage["mars:metadata"]["step"] for parameter in self.parameters: - # values[parameter].append(coverage["ranges"][parameter]["values"]) - if coverage["domain"]["axes"]["t"]["values"][0] not in values[parameter]: - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]] = {} - if ( - coverage["mars:metadata"]["number"] - not in values[parameter][coverage["domain"]["axes"]["t"]["values"][0]] - ): - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]][ - coverage["mars:metadata"]["number"] - ] = {} - values[parameter][coverage["domain"]["axes"]["t"]["values"][0]][coverage["mars:metadata"]["number"]][ - coverage["mars:metadata"]["step"] - ] = coverage["ranges"][parameter]["values"] + values[parameter].setdefault(datetime, {}).setdefault(number, {})[step] = coverage["ranges"][parameter][ + "values" + ] datetimes = sorted(list(set(datetimes))) numbers = sorted(list(set(numbers))) @@ -152,7 +147,10 @@ def to_xarray(self): for j, number in enumerate(numbers): new_values[parameter][i].append([]) for k, step in enumerate(steps): - new_values[parameter][i][j].append(values[parameter][datetime][number][step]) + cell = values[parameter].get(datetime, {}).get(number, {}).get(step) + if cell is None: + cell = [None] * len(longitude) + new_values[parameter][i][j].append(cell) for parameter in self.parameters: dataarray = xr.DataArray(new_values[parameter], dims=dims) diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index a0cf722..44a4064 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -174,7 +174,7 @@ def to_xarray(self): if "mars:metadata" not in self.covjson["coverages"][0]: times = [0] else: - times = sorted({cov["mars:metadata"]["Forecast date"] for cov in self.covjson["coverages"]}) + times = sorted({self.forecast_reference(cov) for cov in self.covjson["coverages"]}) if "mars:metadata" not in self.covjson["coverages"][0]: numbers = [0] else: @@ -213,11 +213,10 @@ def to_xarray(self): # Fill arrays for coverage in self.covjson["coverages"]: if "mars:metadata" not in coverage: - md = {"Forecast date": 0, "number": 0} + t_idx, n_idx = 0, 0 else: - md = coverage["mars:metadata"] - t_idx = times.index(md["Forecast date"]) - n_idx = numbers.index(md.get("number", 0)) + t_idx = times.index(self.forecast_reference(coverage)) + n_idx = numbers.index(coverage["mars:metadata"].get("number", 0)) for pname, prange in coverage["ranges"].items(): arr = np.array(prange["values"]).reshape(prange["shape"]) diff --git a/covjsonkit/decoder/Path.py b/covjsonkit/decoder/Path.py index 6365aa7..7a72e92 100644 --- a/covjsonkit/decoder/Path.py +++ b/covjsonkit/decoder/Path.py @@ -145,21 +145,14 @@ def to_xarray(self): if "step" not in coverage["mars:metadata"]: coverage["mars:metadata"]["step"] = 0 steps.append(coverage["mars:metadata"]["step"]) - datetimes.append(coverage["mars:metadata"]["Forecast date"]) + datetime = self.forecast_reference(coverage) + datetimes.append(datetime) + number = coverage["mars:metadata"]["number"] + step = coverage["mars:metadata"]["step"] for parameter in self.parameters: - # values[parameter].append(coverage["ranges"][parameter]["values"]) - if coverage["mars:metadata"]["Forecast date"] not in values[parameter]: - values[parameter][coverage["mars:metadata"]["Forecast date"]] = {} - if ( - coverage["mars:metadata"]["number"] - not in values[parameter][coverage["mars:metadata"]["Forecast date"]] - ): - values[parameter][coverage["mars:metadata"]["Forecast date"]][ - coverage["mars:metadata"]["number"] - ] = {} - values[parameter][coverage["mars:metadata"]["Forecast date"]][coverage["mars:metadata"]["number"]][ - coverage["mars:metadata"]["step"] - ] = coverage["ranges"][parameter]["values"] + values[parameter].setdefault(datetime, {}).setdefault(number, {})[step] = coverage["ranges"][parameter][ + "values" + ] datetimes = list(set(datetimes)) numbers = list(set(numbers)) @@ -173,7 +166,10 @@ def to_xarray(self): for j, number in enumerate(numbers): new_values[parameter][i].append([]) for k, step in enumerate(steps): - new_values[parameter][i][j].append(values[parameter][datetime][number][step]) + cell = values[parameter].get(datetime, {}).get(number, {}).get(step) + if cell is None: + cell = [None] * len(longitude) + new_values[parameter][i][j].append(cell) for parameter in self.parameters: dataarray = xr.DataArray(new_values[parameter], dims=dims) diff --git a/covjsonkit/decoder/Position.py b/covjsonkit/decoder/Position.py index fa7cc34..3666e8f 100644 --- a/covjsonkit/decoder/Position.py +++ b/covjsonkit/decoder/Position.py @@ -2,6 +2,7 @@ import xarray as xr from .decoder import Decoder +from .TimeSeries import TimeSeries class Position(Decoder): @@ -108,6 +109,9 @@ def to_geojson(self): # function to convert covjson to xarray dataset def to_xarray(self): + if not all("Forecast date" in cov.get("mars:metadata", {}) for cov in self.covjson["coverages"]): + # Reanalysis (efcl) coverages carry only valid times on ``t``. + return TimeSeries._to_xarray_no_forecast_date(self) if self.z_name: dims = ["latitude", "longitude", "levelist", "number", "datetime", "t"] else: diff --git a/covjsonkit/decoder/VerticalProfile.py b/covjsonkit/decoder/VerticalProfile.py index 591a588..e417fe9 100644 --- a/covjsonkit/decoder/VerticalProfile.py +++ b/covjsonkit/decoder/VerticalProfile.py @@ -1,9 +1,6 @@ -from datetime import datetime as dt -from datetime import timedelta - import xarray as xr -from ..encoder.encoder import parse_step_string, sort_step_values +from ..encoder.encoder import sort_step_values, valid_time from .decoder import Decoder @@ -148,8 +145,8 @@ def to_xarray(self): steps = [] for coverage in self.covjson["coverages"]: num.append(coverage["mars:metadata"]["number"]) - datetime.append(coverage["mars:metadata"]["Forecast date"]) - steps.append(coverage["mars:metadata"]["step"]) + datetime.append(self.forecast_reference(coverage)) + steps.append(coverage["mars:metadata"].get("step", 0)) nums = list(set(num)) datetime = list(set(datetime)) @@ -171,12 +168,10 @@ def to_xarray(self): if len(param_values[parameter][domain_idx][i][j]) <= k: param_values[parameter][domain_idx][i][j].append([]) for coverage in self.covjson["coverages"]: - new_step = ( - dt.fromisoformat(date.replace("Z", "")) + timedelta(hours=parse_step_string(step)) - ).isoformat() + "Z" + new_step = valid_time(date, step) if ( coverage["mars:metadata"]["number"] == num - and coverage["mars:metadata"]["Forecast date"] == date + and self.forecast_reference(coverage) == date and coverage["domain"]["axes"][self.x_name]["values"][0] == longitude[0] and coverage["domain"]["axes"][self.y_name]["values"][0] == latitude[0] and coverage["domain"]["axes"]["t"]["values"][0] == new_step diff --git a/covjsonkit/decoder/decoder.py b/covjsonkit/decoder/decoder.py index 6922e46..90c132d 100644 --- a/covjsonkit/decoder/decoder.py +++ b/covjsonkit/decoder/decoder.py @@ -91,6 +91,25 @@ def get_mars_metadata(self): mars_metadata.append(coverage["mars:metadata"]) return mars_metadata + @staticmethod + def forecast_reference(coverage): + """Return a coverage's ``"Forecast date"``, or its first valid time when absent. + + Reanalysis coverages (``class=ce``, ``stream=efcl``) carry no forecast + date; their valid time (from ``t`` or the composite ``t`` coordinate) + identifies them instead. + """ + mars_metadata = coverage.get("mars:metadata", {}) + if "Forecast date" in mars_metadata: + return mars_metadata["Forecast date"] + axes = coverage["domain"]["axes"] + if "t" in axes: + return axes["t"]["values"][0] + composite = axes.get("composite", {}) + if "t" in composite.get("coordinates", []): + return composite["values"][0][composite["coordinates"].index("t")] + return 0 + @abstractmethod def get_ranges(self): pass From 4b5bebfcadc18a5eab5d970dc9cf9cf5e79d5d2b Mon Sep 17 00:00:00 2001 From: awarde96 Date: Thu, 8 Oct 2026 11:01:32 +0200 Subject: [PATCH 10/10] Update tests --- tests/test_encoder_position_from_polytope.py | 4 +- ...st_encoder_separate_datetime_reforecast.py | 55 ++-- ..._encoder_vertical_profile_from_polytope.py | 4 +- tests/test_spec_compliance.py | 5 - tests/test_valid_time.py | 255 ++++++++++++++++++ 5 files changed, 296 insertions(+), 27 deletions(-) create mode 100644 tests/test_valid_time.py diff --git a/tests/test_encoder_position_from_polytope.py b/tests/test_encoder_position_from_polytope.py index 6f2568d..1fcf9a4 100644 --- a/tests/test_encoder_position_from_polytope.py +++ b/tests/test_encoder_position_from_polytope.py @@ -125,8 +125,8 @@ def test_single_step(self): "2t": { "type": "NdArray", "dataType": "float", - "shape": [1], - "axisNames": ["t"], + "shape": [], + "axisNames": [], "values": [264.9], } } diff --git a/tests/test_encoder_separate_datetime_reforecast.py b/tests/test_encoder_separate_datetime_reforecast.py index 2679e84..995a1df 100644 --- a/tests/test_encoder_separate_datetime_reforecast.py +++ b/tests/test_encoder_separate_datetime_reforecast.py @@ -8,7 +8,11 @@ import numpy as np import pytest -from conftest import forecast_separate_datetime_tree, reforecast_separate_datetime_tree +from conftest import ( + forecast_separate_datetime_tree, + reforecast_separate_datetime_tree, + reforecast_separate_datetime_vertical_tree, +) from covjsonkit.api import Covjsonkit @@ -34,23 +38,41 @@ 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 [r["coordinates"] for r in covjson["referencing"]] == [["x", "y"], ["t"]] 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], - ] + composite = cov["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["x", "y"] + # [x, y] == [longitude, latitude] + assert composite["values"] == [[11.0, 48.0], [12.0, 50.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 "step" not in cov["mars:metadata"] assert cov["mars:metadata"]["class"] == "ce" assert cov["mars:metadata"]["stream"] == "efcl" +@pytest.mark.parametrize("feature", ["BoundingBox", "Frame", "Circle", "Shapefile", "polygon"]) +def test_separate_datetime_reforecast_composite_encoders_with_levels(feature): + # Layout product(hdate, levelist, step, time) for a single point. + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + tree = reforecast_separate_datetime_vertical_tree(points, HDATES, TIMES, (500, 850)) + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope_reforecast(tree) + + assert [r["coordinates"] for r in covjson["referencing"]] == [["x", "y"], ["t"], ["z"]] + cov = covjson["coverages"][0] + assert cov["domain"]["axes"]["t"]["values"] == ["2025-07-14T00:00:00Z"] + composite = cov["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["x", "y", "z"] + assert composite["values"] == [[11.0, 48.0, 500], [11.0, 48.0, 850]] + # L0 t0 -> idx0=1, L1 t0 -> idx2=3 + assert cov["ranges"]["t"]["values"] == [1.0, 3.0] + + 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, ... @@ -60,19 +82,16 @@ def test_separate_datetime_reforecast_single_point_two_steps(): 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"] - ] + seen = [(c["domain"]["axes"]["t"]["values"][0], 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), + ("2025-07-14T00:00:00Z", 1.0), + ("2025-07-14T06:00:00Z", 3.0), + ("2025-07-14T12:00:00Z", 2.0), + ("2025-07-14T18:00:00Z", 4.0), + ("2025-07-15T00:00:00Z", 5.0), + ("2025-07-15T06:00:00Z", 7.0), + ("2025-07-15T12:00:00Z", 6.0), + ("2025-07-15T18:00:00Z", 8.0), ] diff --git a/tests/test_encoder_vertical_profile_from_polytope.py b/tests/test_encoder_vertical_profile_from_polytope.py index c9d4f04..c6b11af 100644 --- a/tests/test_encoder_vertical_profile_from_polytope.py +++ b/tests/test_encoder_vertical_profile_from_polytope.py @@ -183,8 +183,8 @@ def test_step_offset(self): "t": { "type": "NdArray", "dataType": "float", - "shape": [1], - "axisNames": ["z"], + "shape": [], + "axisNames": [], "values": [290.0], } } diff --git a/tests/test_spec_compliance.py b/tests/test_spec_compliance.py index ebf752b..4b55c76 100644 --- a/tests/test_spec_compliance.py +++ b/tests/test_spec_compliance.py @@ -364,11 +364,6 @@ def test_bounding_box_geojson_lonlat_z_order(self): gj = Covjsonkit().decode(covjson).to_geojson() assert gj["features"][0]["geometry"]["coordinates"] == [11.0, 48.0, 500] - @pytest.mark.xfail( - reason="Decoder round-trip does not yet support valid-time t axis (encoder-only change); " - "decoder update tracked separately.", - strict=False, - ) @pytest.mark.parametrize("tree_factory", [_surface_forecast_tree, _level_forecast_tree]) def test_bounding_box_legacy_and_new_xarray_equivalent(self, tree_factory): new_covjson = Covjsonkit().encode("CoverageCollection", "BoundingBox").from_polytope(tree_factory()) diff --git a/tests/test_valid_time.py b/tests/test_valid_time.py new file mode 100644 index 0000000..159614b --- /dev/null +++ b/tests/test_valid_time.py @@ -0,0 +1,255 @@ +"""Tests for the shared valid-time helpers and their use across encoders/decoders.""" + +from datetime import timedelta + +import numpy as np +import pandas as pd +import pytest +from conftest import ( + chain, + forecast_tree, + make_point, + node, + reforecast_separate_datetime_tree, + reforecast_separate_datetime_vertical_tree, + tip, +) +from polytope_feature.datacube.tensor_index_tree import TensorIndexTree + +from covjsonkit.api import Covjsonkit +from covjsonkit.encoder.encoder import ( + iso_utc, + range_shape, + reference_datetime, + set_forecast_date, + step_to_timedelta, + time_to_timedelta, + valid_time, +) + +HDATES = (np.datetime64("2025-07-14T00:00:00"), np.datetime64("2025-07-15T00:00:00")) +TIMES = (np.timedelta64(0, "h"), np.timedelta64(12, "h")) + + +@pytest.mark.parametrize( + "step,expected", + [ + (6, timedelta(hours=6)), + (6.5, timedelta(hours=6, minutes=30)), + ("6", timedelta(hours=6)), + ("6h", timedelta(hours=6)), + ("30m", timedelta(minutes=30)), + ("13h30m", timedelta(hours=13, minutes=30)), + ("0-6", timedelta(hours=6)), + ("12-18", timedelta(hours=18)), + (timedelta(hours=3), timedelta(hours=3)), + (np.timedelta64(90, "m"), timedelta(minutes=90)), + (pd.Timedelta(hours=2), timedelta(hours=2)), + ([6], timedelta(hours=6)), + (None, timedelta(0)), + ], +) +def test_step_to_timedelta(step, expected): + assert step_to_timedelta(step) == expected + + +@pytest.mark.parametrize("step", ["abc", "6x", [0, 6]]) +def test_step_to_timedelta_rejects_unknown(step): + with pytest.raises(ValueError): + step_to_timedelta(step) + + +@pytest.mark.parametrize( + "value,expected", + [ + (np.timedelta64(6, "h"), timedelta(hours=6)), + (timedelta(hours=12), timedelta(hours=12)), + ("06:00:00", timedelta(hours=6)), + ("0600", timedelta(hours=6)), + (1230, timedelta(hours=12, minutes=30)), + (6, timedelta(hours=6)), + (None, timedelta(0)), + ], +) +def test_time_to_timedelta(value, expected): + assert time_to_timedelta(value) == expected + + +def test_valid_time_formats_iso_utc(): + assert valid_time("2025-01-01T00:00:00Z", "30m") == "2025-01-01T00:30:00Z" + assert valid_time(np.datetime64("2025-01-01"), 6, np.timedelta64(12, "h")) == "2025-01-01T18:00:00Z" + # Timezone-aware input is converted to UTC; no "+00:00Z" suffix. + assert valid_time(pd.Timestamp("2025-01-01T01:00:00+01:00"), 0) == "2025-01-01T00:00:00Z" + assert iso_utc(pd.Timestamp("2025-01-01", tz="UTC")) == "2025-01-01T00:00:00Z" + assert reference_datetime("2025-01-01Z", "1200") == pd.Timestamp("2025-01-01T12:00:00") + + +def test_set_forecast_date_rules(): + ref = pd.Timestamp("2025-01-01T12:00:00") + assert set_forecast_date({"class": "od"}, ref)["Forecast date"] == "2025-01-01T12:00:00Z" + efcl = set_forecast_date({"class": "ce", "stream": "efcl", "Forecast date": "x"}, ref) + assert "Forecast date" not in efcl + efas = set_forecast_date({"class": "ce", "stream": "efas", "date": "2025-01-01"}, ref) + assert efas["Forecast date"] == "2025-01-01T12:00:00Z" + assert "date" not in efas + + +def test_range_shape(): + assert range_shape([1.0], "t") == ([], []) + assert range_shape([1.0, 2.0], "t") == ([2], ["t"]) + + +@pytest.mark.parametrize("feature", ["timeseries", "position", "boundingbox", "verticalprofile", "grid", "path"]) +@pytest.mark.parametrize("step,expected", [("30m", "2025-01-01T00:30:00Z"), ("6", "2025-01-01T06:00:00Z")]) +def test_string_steps_from_polytope(feature, step, expected): + covjson = getattr(Covjsonkit().encode("CoverageCollection", feature), "from_polytope")( + forecast_tree([(48.0, 11.0, [1.0])], step=(step,)) + ) + axes = covjson["coverages"][0]["domain"]["axes"] + t = axes["t"]["values"] if "t" in axes else [axes["composite"]["values"][0][0]] + assert t == [expected] + + +@pytest.mark.parametrize("feature", ["timeseries", "position", "verticalprofile"]) +@pytest.mark.parametrize("step,expected", [("30m", "2025-07-14T00:30:00Z"), ("6", "2025-07-14T06:00:00Z")]) +def test_string_steps_from_polytope_reforecast(feature, step, expected): + tree = reforecast_separate_datetime_tree([(48.0, 11.0, [1.0])], HDATES[:1], TIMES[:1], step=(step,)) + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope_reforecast(tree) + assert covjson["coverages"][0]["domain"]["axes"]["t"]["values"] == [expected] + + +def _forecast_tree_with_time(times): + tree = chain( + TensorIndexTree(), + node("class", ("od",)), + node("date", (np.datetime64("2025-01-01T00:00:00"),)), + node("domain", ("g",)), + node("expver", ("0001",)), + node("levtype", ("sfc",)), + node("param", ("167",)), + node("step", (0, 6)), + node("stream", ("oper",)), + node("time", times), + node("type", ("fc",)), + ) + tip(tree).add_child(make_point(48.0, 11.0, [1.0, 2.0])) + return tree + + +def test_walk_tree_folds_single_time_into_date(): + covjson = ( + Covjsonkit() + .encode("CoverageCollection", "boundingbox") + .from_polytope(_forecast_tree_with_time((np.timedelta64(12, "h"),))) + ) + seen = [(c["mars:metadata"]["Forecast date"], c["domain"]["axes"]["t"]["values"]) for c in covjson["coverages"]] + assert seen == [ + ("2025-01-01T12:00:00Z", ["2025-01-01T12:00:00Z"]), + ("2025-01-01T12:00:00Z", ["2025-01-01T18:00:00Z"]), + ] + + +def test_walk_tree_rejects_multiple_times(): + with pytest.raises(ValueError, match="from_polytope_step or from_polytope_reforecast"): + Covjsonkit().encode("CoverageCollection", "boundingbox").from_polytope( + _forecast_tree_with_time((np.timedelta64(0, "h"), np.timedelta64(12, "h"))) + ) + + +def _step_tree(step): + return chain( + TensorIndexTree(), + node("class", ("d1",)), + node("date", (np.datetime64("2025-01-01T00:00:00"),)), + node("expver", ("0001",)), + node("levtype", ("sfc",)), + node("param", ("167",)), + node("step", (step,)), + node("stream", ("clte",)), + node("time", (timedelta(0), timedelta(hours=6))), + make_point(48.0, 11.0, [264.0, 263.0]), + ) + + +@pytest.mark.parametrize("feature", ["timeseries", "position"]) +def test_step_path_valid_time_is_iso_and_includes_step(feature): + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope_step(_step_tree(1)) + assert covjson["coverages"][0]["domain"]["axes"]["t"]["values"] == [ + "2025-01-01T01:00:00Z", + "2025-01-01T07:00:00Z", + ] + + +def test_grid_reforecast_separate_datetime(): + points = [(48.0, 11.0, [1.0, 2.0, 3.0, 4.0]), (48.0, 12.0, [5.0, 6.0, 7.0, 8.0])] + tree = reforecast_separate_datetime_tree(points, HDATES, TIMES) + covjson = Covjsonkit().encode("CoverageCollection", "grid").from_polytope_reforecast(tree) + + assert [r["coordinates"] for r in covjson["referencing"]] == [["x", "y"], ["t"]] + assert len(covjson["coverages"]) == 4 + cov = covjson["coverages"][0] + assert cov["domain"]["axes"] == { + "t": {"values": ["2025-07-14T00:00:00Z"]}, + "y": {"values": [48.0]}, + "x": {"values": [11.0, 12.0]}, + } + assert cov["ranges"]["2t"]["shape"] == [1, 1, 2] + assert cov["ranges"]["2t"]["axisNames"] == ["t", "y", "x"] + assert cov["ranges"]["2t"]["values"] == [1.0, 5.0] + assert "Forecast date" not in cov["mars:metadata"] + + +def test_grid_reforecast_steps_share_a_coverage(): + 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", "grid").from_polytope_reforecast(tree) + cov = covjson["coverages"][0] + # h0 t0: step 0 -> idx0=1, step 6 -> idx2=3 + assert cov["domain"]["axes"]["t"]["values"] == ["2025-07-14T00:00:00Z", "2025-07-14T06:00:00Z"] + assert cov["ranges"]["2t"]["values"] == [1.0, 3.0] + + +def test_path_reforecast_separate_datetime(): + points = [(48.0, 11.0, [1.0, 2.0, 3.0, 4.0]), (50.0, 12.0, [5.0, 6.0, 7.0, 8.0])] + tree = reforecast_separate_datetime_tree(points, HDATES, TIMES) + covjson = Covjsonkit().encode("CoverageCollection", "path").from_polytope_reforecast(tree) + + assert [r["coordinates"] for r in covjson["referencing"]] == [["x", "y"], ["t"]] + cov = covjson["coverages"][0] + composite = cov["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["t", "x", "y"] + assert composite["values"] == [["2025-07-14T00:00:00Z", 11.0, 48.0], ["2025-07-14T00:00:00Z", 12.0, 50.0]] + assert cov["ranges"]["2t"]["values"] == [1.0, 5.0] + + +def test_path_reforecast_with_levels(): + points = [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])] + tree = reforecast_separate_datetime_vertical_tree(points, HDATES, TIMES, (500, 850)) + covjson = Covjsonkit().encode("CoverageCollection", "path").from_polytope_reforecast(tree) + composite = covjson["coverages"][0]["domain"]["axes"]["composite"] + assert composite["coordinates"] == ["t", "x", "y", "z"] + assert composite["values"] == [["2025-07-14T00:00:00Z", 11.0, 48.0, 500], ["2025-07-14T00:00:00Z", 11.0, 48.0, 850]] + + +def test_multipoint_multi_step_round_trips_to_xarray(): + tree = forecast_tree([(48.0, 11.0, [1.0, 2.0]), (50.0, 12.0, [3.0, 4.0])], step=(0, 6)) + covjson = Covjsonkit().encode("CoverageCollection", "boundingbox").from_polytope(tree) + ds = Covjsonkit().decode(covjson).to_xarray() + assert list(ds["steps"].values) == [0, 6] + assert list(ds["datetimes"].values) == ["2025-01-01T00:00:00Z"] + assert ds["2t"].values.tolist() == [[[[1.0, 3.0], [2.0, 4.0]]]] + + +@pytest.mark.parametrize("feature", ["boundingbox", "grid", "path", "verticalprofile", "position"]) +def test_efcl_reforecast_round_trips_to_xarray(feature): + if feature == "verticalprofile": + tree = reforecast_separate_datetime_vertical_tree( + [(48.0, 11.0, [1, 2, 3, 4, 5, 6, 7, 8])], HDATES, TIMES, (500, 850) + ) + else: + tree = reforecast_separate_datetime_tree([(48.0, 11.0, [1.0, 2.0, 3.0, 4.0])], HDATES, TIMES) + covjson = Covjsonkit().encode("CoverageCollection", feature).from_polytope_reforecast(tree) + assert all("Forecast date" not in c["mars:metadata"] for c in covjson["coverages"]) + decoder = Covjsonkit().decode(covjson) + decoder.to_xarray() + decoder.to_geojson()