Feature/polytope indicies

covjsonkit/decoder/Grid.py

@@ -133,6 +133,103 @@ 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": ("indicies", indicies),
+            },
+        )
+
+        # Attach parameter metadata
+        for pname, param in parameters.items():
+            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):
+        """
+        Convert a Grid domainType CoverageJSON CoverageCollection into an xarray.Dataset.
+
         Dimensions:
         - time (Forecast date)
         - number (ensemble member)

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:

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:

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:

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):
@@ -60,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",
@@ -92,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
@@ -132,15 +123,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 +183,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 +299,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

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:

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: