From 9f354d4e09bbb3167d0dece7511935c6232e5ee1 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 17 Oct 2025 12:18:18 +0000 Subject: [PATCH 1/5] Add encoder for grid using indicies --- covjsonkit/encoder/Grid.py | 123 +++++++++++++++++++++++++++++----- covjsonkit/encoder/encoder.py | 3 + 2 files changed, 111 insertions(+), 15 deletions(-) diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index d5774f4..ec190e5 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -30,10 +30,8 @@ def add_domain(self, coverage, coords): 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"]["indicies"] = {} + coverage["domain"]["axes"]["indicies"]["values"] = coords["indicies"] coverage["domain"]["axes"]["levelist"] = {} coverage["domain"]["axes"]["levelist"]["values"] = coords["levelist"] @@ -44,7 +42,7 @@ def add_range(self, coverage, values): 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]["axisNames"] = ["t", "levelist", "indicies"] coverage["ranges"][param]["values"] = values[parameter] # [values[parameter]] def add_mars_metadata(self, coverage, metadata): @@ -132,15 +130,19 @@ def from_polytope(self, result): fields["step"] = [0] fields["dates"] = [] fields["levels"] = [0] + fields["indicies"] = [] self.walk_tree(result, fields, coords, mars_metadata, range_dict) + print(fields) + print(range_dict) + print(coords) 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"], + "coordinates": ["indicies", "levelist"], "system": { "type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", @@ -188,19 +190,12 @@ def from_polytope(self, result): coordinates[date] = {} coordinates[date]["t"] = list(fields["step"]) coordinates[date]["levelist"] = list(fields["levels"]) - coordinates[date]["latitude"] = [] - coordinates[date]["longitude"] = [] - 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"]) + coordinates[date]["indicies"] = fields["indicies"] 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"]), + len(coordinates[fields["dates"][0]]["indicies"]), ] for date in combined_dict.keys(): @@ -311,3 +306,101 @@ def from_polytope_step(self, result): logging.debug("Coverage creation: %s", delta) # noqa: E501 return self.covjson + + def from_polytope_old_latlon(self, result): + + coords = {} + mars_metadata = {} + range_dict = {} + fields = {} + fields["lat"] = 0 + fields["param"] = 0 + fields["number"] = [0] + fields["step"] = [0] + fields["dates"] = [] + fields["levels"] = [0] + + self.walk_tree(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": ["latitude", "longitude", "levelist"], + "system": { + "type": "GeographicCRS", + "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", + }, + } + ) + + combined_dict = {} + + for date in fields["dates"]: + if date not in combined_dict: + combined_dict[date] = {} + for level in fields["levels"]: + for num in fields["number"]: + if num not in combined_dict[date]: + combined_dict[date][num] = {} + 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 + combined_dict[date][num][para][s] += range_dict[key] + + 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 + + logging.debug("The fields retrieved were: %s", fields) # noqa: E501 + logging.debug("The range_dict created was: %s", range_dict) # noqa: E501 + + coordinates = {} + coordinates["t"] = list(fields["step"]) + + for date in coords.keys(): + coordinates[date] = {} + coordinates[date]["t"] = list(fields["step"]) + coordinates[date]["levelist"] = list(fields["levels"]) + coordinates[date]["latitude"] = [] + coordinates[date]["longitude"] = [] + 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"]), + ] + + for date in combined_dict.keys(): + for num in combined_dict[date].keys(): + val_dict = {} + for para in combined_dict[date][num].keys(): + val_dict[para] = [] + for step in combined_dict[date][num][para].keys(): + val_dict[para].extend(combined_dict[date][num][para][step]) + mm = mars_metadata.copy() + mm["number"] = num + mm["step"] = step + mm["Forecast date"] = date + self.add_coverage(mm, coordinates[date], val_dict) + + return self.covjson diff --git a/covjsonkit/encoder/encoder.py b/covjsonkit/encoder/encoder.py index 59c1011..d6f5781 100644 --- a/covjsonkit/encoder/encoder.py +++ b/covjsonkit/encoder/encoder.py @@ -221,6 +221,9 @@ def append_composite_coords(dates, tree_values, lat, coords): del range_dict[key] else: tree.result = [float(val) if val is not None else val for val in tree.result] + print(tree.indexes) + if "indicies" in fields: + fields["indicies"].extend(tree.indexes) level_len = len(tree.result) / len(fields["levels"]) num_len = level_len / len(fields["number"]) para_len = num_len / len(fields["param"]) From e3f5c7be660d82bcc987f355e588d796865c71a3 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 17 Oct 2025 12:35:15 +0000 Subject: [PATCH 2/5] Add to and from xarray for grid --- covjsonkit/decoder/Grid.py | 91 +++++++++++++++++++++++++++++++ covjsonkit/encoder/BoundingBox.py | 2 +- covjsonkit/encoder/Circle.py | 2 +- covjsonkit/encoder/Frame.py | 2 +- covjsonkit/encoder/Grid.py | 49 +++++++---------- covjsonkit/encoder/Shapefile.py | 2 +- covjsonkit/encoder/Wkt.py | 2 +- 7 files changed, 117 insertions(+), 33 deletions(-) diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index 1abf68f..96c633b 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -133,6 +133,97 @@ def to_xarray(self): """ Convert a Grid domainType CoverageJSON CoverageCollection into an xarray.Dataset. + Dimensions: + - time (Forecast date) + - number (ensemble member) + - step ('t' axis) + - level ('levelist') + - indicies + """ + if self.covjson["type"] != "CoverageCollection": + raise ValueError("Expected CoverageCollection as root object") + + parameters = self.covjson.get("parameters", {}) + + # Collect metadata for unique coords + 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"]}) + if "mars:metadata" not in self.covjson["coverages"][0]: + numbers = [0] + else: + numbers = sorted({cov["mars:metadata"].get("number", 0) for cov in self.covjson["coverages"]}) + + # Initialize coords from first coverage + first_cov = self.covjson["coverages"][0] + domain = first_cov["domain"]["axes"] + + if "indicies" in domain: + i_coords = "indicies" + if "levelist" in domain: + z_coords = "levelist" + else: + z_coords = "z" + + steps = np.array(domain.get("t", {}).get("values", [0])) + levels = np.array(domain.get(z_coords, {}).get("values", [0])) + indicies = np.array(domain[i_coords]["values"]) + + # Prepare arrays for each parameter + data_arrays = { + pname: np.full( + (len(times), len(numbers), len(steps), len(levels), len(indicies)), + np.nan, + dtype=float, + ) + for pname in first_cov["ranges"].keys() + } + + # Fill arrays + for coverage in self.covjson["coverages"]: + if "mars:metadata" not in coverage: + md = {"Forecast date": 0, "number": 0} + else: + md = coverage["mars:metadata"] + t_idx = times.index(md["Forecast date"]) + n_idx = numbers.index(md.get("number", 0)) + + for pname, prange in coverage["ranges"].items(): + arr = np.array(prange["values"]).reshape(prange["shape"]) + data_arrays[pname][t_idx, n_idx, :, :, :] = arr + + # Build xarray Dataset + xr_vars = { + pname: ( + ["datetimes", "number", "steps", "levelist", "indicies"], + arr, + ) + for pname, arr in data_arrays.items() + } + + ds = xr.Dataset( + xr_vars, + coords={ + "datetimes": ("datetimes", np.array(times)), + "number": ("number", np.array(numbers)), + "steps": ("steps", steps), + "levelist": ("levelist", levels), + "indicies": ("latitude", indicies), + }, + ) + + # Attach parameter metadata + for pname, param in parameters.items(): + for k, v in param.items(): + ds[pname].attrs[k] = v + + return ds + + def to_xarray_old(self): + """ + Convert a Grid domainType CoverageJSON CoverageCollection into an xarray.Dataset. + Dimensions: - time (Forecast date) - number (ensemble member) diff --git a/covjsonkit/encoder/BoundingBox.py b/covjsonkit/encoder/BoundingBox.py index a2e1502..d5521ec 100644 --- a/covjsonkit/encoder/BoundingBox.py +++ b/covjsonkit/encoder/BoundingBox.py @@ -56,7 +56,7 @@ def from_xarray(self, dataset): """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] if "latitude" in dataset.coords: diff --git a/covjsonkit/encoder/Circle.py b/covjsonkit/encoder/Circle.py index 8061351..c99b8f4 100644 --- a/covjsonkit/encoder/Circle.py +++ b/covjsonkit/encoder/Circle.py @@ -53,7 +53,7 @@ def from_xarray(self, dataset): """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] if "latitude" in dataset.coords: diff --git a/covjsonkit/encoder/Frame.py b/covjsonkit/encoder/Frame.py index 89e5af8..e582640 100644 --- a/covjsonkit/encoder/Frame.py +++ b/covjsonkit/encoder/Frame.py @@ -53,7 +53,7 @@ def from_xarray(self, dataset): """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] if "latitude" in dataset.coords: diff --git a/covjsonkit/encoder/Grid.py b/covjsonkit/encoder/Grid.py index ec190e5..8b480e2 100644 --- a/covjsonkit/encoder/Grid.py +++ b/covjsonkit/encoder/Grid.py @@ -58,24 +58,18 @@ def from_xarray(self, dataset): """ 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 "indicies" in dataset.coords: + i_coord = "indicies" if "levelist" in dataset.coords: z_coord = "levelist" # Add reference system self.add_reference( { - "coordinates": [x_coord, y_coord, z_coord], + "coordinates": [i_coord, z_coord], "system": { "type": "GeographicCRS", "id": "http://www.opengis.net/def/crs/OGC/1.3/CRS84", @@ -90,30 +84,29 @@ def from_xarray(self, dataset): # Prepare coordinates coords = { "t": [str(x) for x in dataset["steps"].values], - "latitude": dataset["latitude"].values.tolist(), - "longitude": dataset["longitude"].values.tolist(), + "indicies": dataset["indicies"].values.tolist(), "levelist": dataset["levelist"].values.tolist(), } - self.shp = [len(coords["t"]), len(coords["levelist"]), len(coords["latitude"]), len(coords["longitude"])] + self.shp = [len(coords["t"]), len(coords["levelist"]), len(coords["indicies"])] for datetime in dataset["datetimes"].values: for num in dataset["number"].values: - for step in dataset["steps"].values: - dv_dict = {} - mars_metadata = {metadata: dataset.attrs[metadata] for metadata in dataset.attrs} - mars_metadata["number"] = int(num) - mars_metadata["step"] = int(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() - 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) + # for step in dataset["steps"].values: + dv_dict = {} + mars_metadata = {metadata: dataset.attrs[metadata] for metadata in dataset.attrs} + mars_metadata["number"] = int(num) + # mars_metadata["step"] = int(step) + mars_metadata["Forecast date"] = str(datetime) + for dv in dataset.data_vars: + nested_list = dataset[dv].sel(datetimes=datetime, number=num).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) # Return the generated CoverageJSON return self.covjson diff --git a/covjsonkit/encoder/Shapefile.py b/covjsonkit/encoder/Shapefile.py index 2657f34..199ff72 100644 --- a/covjsonkit/encoder/Shapefile.py +++ b/covjsonkit/encoder/Shapefile.py @@ -53,7 +53,7 @@ def from_xarray(self, dataset): """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] if "latitude" in dataset.coords: diff --git a/covjsonkit/encoder/Wkt.py b/covjsonkit/encoder/Wkt.py index 4aecee9..075781f 100644 --- a/covjsonkit/encoder/Wkt.py +++ b/covjsonkit/encoder/Wkt.py @@ -57,7 +57,7 @@ def from_xarray(self, dataset): """ self.covjson["type"] = "CoverageCollection" - self.covjson["domainType"] = "PointSeries" + self.covjson["domainType"] = "MultiPoint" self.covjson["coverages"] = [] if "latitude" in dataset.coords: From 11ba815eced451d284b58bee335e8dd617279a27 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Fri, 17 Oct 2025 13:06:03 +0000 Subject: [PATCH 3/5] Update tests --- tests/test_encoder_grid.py | 3 +-- tests/test_geotiff.py | 10 +++++++++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/test_encoder_grid.py b/tests/test_encoder_grid.py index 07b109f..6aedf26 100644 --- a/tests/test_encoder_grid.py +++ b/tests/test_encoder_grid.py @@ -191,8 +191,7 @@ def test_add_coverage(self): } coords = {} coords["t"] = ["2017-01-01T00:00:00Z"] - coords["latitude"] = [20, 21, 17] - coords["longitude"] = [1, 2, 3] + coords["indicies"] = [20, 21, 17] coords["levelist"] = [1] value = {"2t": [111, 222, 333, 444, 555, 666, 777, 888, 999]} encoder.shp = [1, 1, 3, 3] diff --git a/tests/test_geotiff.py b/tests/test_geotiff.py index df0f4ea..c09e8b2 100644 --- a/tests/test_geotiff.py +++ b/tests/test_geotiff.py @@ -14,8 +14,16 @@ def setup_method(self): self.covjson = json.load(f) def test_geotiff_multipoint(self): - with pytest.raises(ImportError): + try: cov = Covjsonkit().decode(self.covjson) cov.to_geotiff() + assert os.path.exists("multipoint_2t.tif") + assert os.path.exists("multipoint_10u.tif") os.remove("multipoint_2t.tif") os.remove("multipoint_10u.tif") + except ImportError: + with pytest.raises(ImportError): + cov = Covjsonkit().decode(self.covjson) + cov.to_geotiff() + os.remove("multipoint_2t.tif") + os.remove("multipoint_10u.tif") From d17d051ff08a9d5ef85056ee13eea751f0420576 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Mon, 20 Oct 2025 13:19:17 +0000 Subject: [PATCH 4/5] Add metadata when converting to xarray --- covjsonkit/decoder/Grid.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index 96c633b..30ef888 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -218,6 +218,12 @@ def to_xarray(self): for k, v in param.items(): ds[pname].attrs[k] = v + for mars_metadata in self.mars_metadata[0]: + ds.attrs[mars_metadata] = self.mars_metadata[0][mars_metadata] + + # Add date attribute + ds.attrs["date"] = self.get_coordinates()["t"]["values"][0] + return ds def to_xarray_old(self): From c8f73580545cc204e6836a700ac3a45e5aa4ec31 Mon Sep 17 00:00:00 2001 From: awarde96 Date: Mon, 20 Oct 2025 13:57:37 +0000 Subject: [PATCH 5/5] Remove reference to latitude --- covjsonkit/decoder/Grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/covjsonkit/decoder/Grid.py b/covjsonkit/decoder/Grid.py index 30ef888..c5564b7 100644 --- a/covjsonkit/decoder/Grid.py +++ b/covjsonkit/decoder/Grid.py @@ -209,7 +209,7 @@ def to_xarray(self): "number": ("number", np.array(numbers)), "steps": ("steps", steps), "levelist": ("levelist", levels), - "indicies": ("latitude", indicies), + "indicies": ("indicies", indicies), }, )