Correct extent in bbox, polygon and cube:dimensions

deep_code/tests/utils/test_dataset_stac_generator.py

@@ -36,8 +36,8 @@ def setUp(self):
         """Set up a mock dataset and generator."""
         self.mock_dataset = Dataset(
             coords={
-                "lon": ("lon", np.linspace(-180, 180, 10)),
-                "lat": ("lat", np.linspace(-90, 90, 5)),
+                "lon": ("lon", np.arange(-175, 176, 10)),
+                "lat": ("lat", np.arange(-87.5, 88, 5)),
                 "time": (
                     "time",
                     [
@@ -50,7 +50,7 @@ def setUp(self):
             data_vars={
                 "var1": (
                     ("time", "lat", "lon"),
-                    np.random.rand(2, 5, 10),
+                    np.random.rand(2, 36, 36),
                     {
                         "description": "dummy",
                         "standard_name": "var1",
@@ -59,7 +59,7 @@ def setUp(self):
                 ),
                 "var2": (
                     ("time", "lat", "lon"),
-                    np.random.rand(2, 5, 10),
+                    np.random.rand(2, 36, 36),
                     {
                         "description": "dummy",
                         "standard_name": "var2",
@@ -103,9 +103,19 @@ def test_open_dataset(self):
 
     def test_get_spatial_extent(self):
         """Test spatial extent extraction."""
+
+        # center
         extent = self.generator._get_spatial_extent(self.mock_dataset)
         self.assertEqual(extent.bboxes[0], [-180.0, -90.0, 180.0, 90.0])
 
+        # left
+        extent = self.generator._get_spatial_extent(self.mock_dataset, "left")
+        self.assertEqual(extent.bboxes[0], [-175.0, -87.5, 185.0, 92.5])
+
+        # right
+        extent = self.generator._get_spatial_extent(self.mock_dataset, "right")
+        self.assertEqual(extent.bboxes[0], [-185.0, -92.5, 175.0, 87.5])
+
     def test_get_temporal_extent(self):
         """Test temporal extent extraction."""
         extent = self.generator._get_temporal_extent(self.mock_dataset)
@@ -751,16 +761,16 @@ def test_spatial_extent_longitude_latitude(self, mock_open_ds):
         mock_open_ds.return_value = ds
         gen = self._make_generator(ds)
         extent = gen._get_spatial_extent(ds)
-        self.assertAlmostEqual(extent.bboxes[0][0], -10.0)
-        self.assertAlmostEqual(extent.bboxes[0][1], -5.0)
+        self.assertAlmostEqual(extent.bboxes[0][0], -15.0)
+        self.assertAlmostEqual(extent.bboxes[0][1], -10.0)
 
     @patch("deep_code.utils.dataset_stac_generator.open_dataset")
     def test_spatial_extent_x_y(self, mock_open_ds):
         ds = self._make_dataset("x_y")
         mock_open_ds.return_value = ds
         gen = self._make_generator(ds)
         extent = gen._get_spatial_extent(ds)
-        self.assertAlmostEqual(extent.bboxes[0][0], 0.0)
+        self.assertAlmostEqual(extent.bboxes[0][0], -25.0)
 
     @patch("deep_code.utils.dataset_stac_generator.open_dataset")
     def test_spatial_extent_unknown_coords_raises(self, mock_open_ds):
@@ -861,7 +871,7 @@ class TestPRRCollection(unittest.TestCase):
     def setUp(self):
         self.dataset = Dataset(
             coords={
-                "lon": ("lon", np.linspace(-20, 20, 4)),
+                "lon": ("lon", np.linspace(-20, 20, 3)),
                 "lat": ("lat", np.linspace(-10, 10, 3)),
                 "time": (
                     "time",
@@ -874,12 +884,12 @@ def setUp(self):
             data_vars={
                 "sst": (
                     ("time", "lat", "lon"),
-                    np.random.rand(2, 3, 4),
+                    np.random.rand(2, 3, 3),
                     {"units": "K", "long_name": "Sea surface temperature"},
                 ),
                 "chl": (
                     ("time", "lat", "lon"),
-                    np.random.rand(2, 3, 4),
+                    np.random.rand(2, 3, 3),
                     {"units": "mg m-3"},
                 ),
                 # A CRS var that must be excluded from cube:variables and drive EPSG.
@@ -938,15 +948,39 @@ def test_get_epsg_default_4326(self, mock_open_ds):
         )
         self.assertEqual(gen._get_crs(ds).to_epsg(), 4326)
 
-    def test_get_cube_dimensions(self):
+    def test_get_cube_dimensions_center(self):
         dims = self.gen._get_cube_dimensions(self.dataset)
         self.assertEqual(set(dims), {"lon", "lat", "time"})
         self.assertEqual(dims["lon"]["type"], "spatial")
         self.assertEqual(dims["lon"]["axis"], "x")
         self.assertEqual(dims["lon"]["reference_system"], 3035)
-        self.assertEqual(dims["lon"]["extent"], [-20.0, 20.0])
+        self.assertEqual(dims["lon"]["extent"], [-30.0, 30.0])
+        self.assertEqual(dims["lat"]["axis"], "y")
+        self.assertEqual(dims["lat"]["extent"], [-15.0, 15.0])
+        self.assertEqual(dims["time"]["type"], "temporal")
+        self.assertEqual(len(dims["time"]["extent"]), 2)
+
+    def test_get_cube_dimensions_left(self):
+        dims = self.gen._get_cube_dimensions(self.dataset, "left")
+        self.assertEqual(set(dims), {"lon", "lat", "time"})
+        self.assertEqual(dims["lon"]["type"], "spatial")
+        self.assertEqual(dims["lon"]["axis"], "x")
+        self.assertEqual(dims["lon"]["reference_system"], 3035)
+        self.assertEqual(dims["lon"]["extent"], [-20.0, 40.0])
+        self.assertEqual(dims["lat"]["axis"], "y")
+        self.assertEqual(dims["lat"]["extent"], [-10.0, 20.0])
+        self.assertEqual(dims["time"]["type"], "temporal")
+        self.assertEqual(len(dims["time"]["extent"]), 2)
+
+    def test_get_cube_dimensions_right(self):
+        dims = self.gen._get_cube_dimensions(self.dataset, "right")
+        self.assertEqual(set(dims), {"lon", "lat", "time"})
+        self.assertEqual(dims["lon"]["type"], "spatial")
+        self.assertEqual(dims["lon"]["axis"], "x")
+        self.assertEqual(dims["lon"]["reference_system"], 3035)
+        self.assertEqual(dims["lon"]["extent"], [-40.0, 20.0])
         self.assertEqual(dims["lat"]["axis"], "y")
-        self.assertEqual(dims["lat"]["extent"], [-10.0, 10.0])
+        self.assertEqual(dims["lat"]["extent"], [-20.0, 10.0])
         self.assertEqual(dims["time"]["type"], "temporal")
         self.assertEqual(len(dims["time"]["extent"]), 2)
 

deep_code/tools/new.py

@@ -86,6 +86,7 @@ def generate_dataset_template(output_path: str | None = None) -> str:
             "osc_project_title": "[OPTIONAL: display title of the OSC project as it appears in the catalog (e.g. DeepESDL). Defaults to a formatted version of osc_project if omitted]",
             "access_link_root": "[OPTIONAL: public S3 URL of the Zarr store — defaults to s3://deep-esdl-public]",
             "collection_title": "[OPTIONAL: title present in the collection and in the STAC browser]",
+            "coord_position": "[OPTIONAL: Position of the coordinates within each grid cell. 'center' assumes coordinates represent cell centers, 'left' assumes they represent the left/bottom edge, and 'right' assumes they represent the right/top edge.]",
             "cf_parameter": [
                 {"name": "[OPTIONAL: CF standard name]", "units": "[unit string]"}
             ],

deep_code/tools/prr.py

@@ -93,6 +93,7 @@ def generate_prr_collection(
         thumbnail_media_type=config.get("thumbnail_media_type"),
         sci_doi=config.get("sci_doi"),
         sci_citation=config.get("sci_citation"),
+        coord_position=config.get("coord_position") or "center",
     )
 
     out_dir = output_dir or config.get("prr_output_dir") or "prr"

deep_code/utils/dataset_stac_generator.py

@@ -7,7 +7,7 @@
 import logging
 from dataclasses import dataclass
 from datetime import datetime, timezone
-from typing import Any
+from typing import Any, Literal
 
 import numpy as np
 import pandas as pd
@@ -62,6 +62,10 @@ class OscDatasetStacGenerator:
         osc_missions: List of satellite missions associated with the dataset.
         cf_params: CF metadata parameters for the dataset.
         osc_project: OSC project identifier (default: "deep-earth-system-data-lab").
+        coord_position: Position of the coordinates within each grid cell.
+            ``"center"`` assumes coordinates represent cell centers,
+            ``"left"`` assumes they represent the left/bottom edge, and
+            ``"right"`` assumes they represent the right/top edge.
     """
 
     def __init__(
@@ -92,6 +96,7 @@ def __init__(
         thumbnail_media_type: str | None = None,
         sci_doi: str | None = None,
         sci_citation: str | None = None,
+        coord_position: Literal["left", "center", "right"] = "center",
     ):
         if " " in collection_id:
             raise ValueError(
@@ -128,10 +133,28 @@ def __init__(
         self.thumbnail_media_type = thumbnail_media_type
         self.sci_doi = sci_doi
         self.sci_citation = sci_citation
+        self.coord_position = coord_position
         self.logger = logging.getLogger(__name__)
 
-    def _get_spatial_extent(self, dataset: xr.Dataset) -> SpatialExtent:
-        """Extract the spatial extent and return it in EPSG:4326."""
+    def _get_spatial_extent(
+        self,
+        dataset: xr.Dataset,
+        coord_position: Literal["left", "center", "right"] = "center",
+    ) -> SpatialExtent:
+        """Extract the spatial extent and return it in EPSG:4326.
+
+        Args:
+            dataset: Input dataset.
+            coord_position: Position of the coordinates within each grid cell.
+                ``"center"`` assumes coordinates represent cell centers,
+                ``"left"`` assumes they represent the left/bottom edge, and
+                ``"right"`` assumes they represent the right/top edge.
+        """
+        if coord_position not in {"left", "center", "right"}:
+            raise ValueError(
+                f"Invalid coord_position: {coord_position!r}. "
+                "Must be 'left', 'center', or 'right'."
+            )
 
         if {"lon", "lat"}.issubset(dataset.coords):
             x_name, y_name = "lon", "lat"
@@ -145,10 +168,29 @@ def _get_spatial_extent(self, dataset: xr.Dataset) -> SpatialExtent:
                 "('lon', 'lat'), ('longitude', 'latitude'), or ('x', 'y')."
             )
 
-        x_min = float(dataset[x_name].min())
-        x_max = float(dataset[x_name].max())
-        y_min = float(dataset[y_name].min())
-        y_max = float(dataset[y_name].max())
+        x = dataset[x_name]
+        y = dataset[y_name]
+
+        x_min = float(x.min())
+        x_max = float(x.max())
+        y_min = float(y.min())
+        y_max = float(y.max())
+
+        # Calculate the grid resolution from adjacent coordinates.
+        x_res = float(abs(x.diff(x_name).median()))
+        y_res = float(abs(y.diff(y_name).median()))
+
+        if coord_position == "center":
+            x_min -= x_res / 2
+            x_max += x_res / 2
+            y_min -= y_res / 2
+            y_max += y_res / 2
+        elif coord_position == "left":
+            x_max += x_res
+            y_max += y_res
+        elif coord_position == "right":
+            x_min -= x_res
+            y_min -= y_res
 
         crs = self._get_crs(dataset)
 
@@ -623,7 +665,7 @@ def build_zarr_stac_item(
             f"for collection '{self.collection_id}'."
         )
         dataset = open_dataset(item_config.dataset_id, logger=self.logger)
-        spatial_extent = self._get_spatial_extent(dataset)
+        spatial_extent = self._get_spatial_extent(dataset, self.coord_position)
         temporal_extent = self._get_temporal_extent(dataset)
         general_metadata = self._get_general_metadata(dataset)
 
@@ -793,45 +835,77 @@ def _get_crs(dataset: xr.Dataset) -> pyproj.CRS:
                                 pass
         return pyproj.CRS.from_epsg(4326)
 
-    def _get_cube_dimensions(self, dataset: xr.Dataset) -> dict[str, dict]:
+    def _get_cube_dimensions(
+        self,
+        dataset: xr.Dataset,
+        coord_position: Literal["left", "center", "right"] = "center",
+    ) -> dict[str, dict]:
         """Build the ``cube:dimensions`` object from the dataset coordinates.
 
         Follows the datacube STAC extension: horizontal spatial dimensions are
         classified by axis (x/y), the ``time`` coordinate becomes a temporal
         dimension, and any remaining index coordinate is emitted as an
         additional dimension.
+
+        Args:
+            dataset: Dataset from which to extract the cube dimensions.
+            coord_position: Position of spatial coordinates within each grid cell.
+                ``"center"`` assumes cell-center coordinates,
+                ``"left"`` assumes left/bottom edge coordinates, and
+                ``"right"`` assumes right/top edge coordinates.
         """
+        if coord_position not in {"left", "center", "right"}:
+            raise ValueError(
+                f"Invalid coord_position: {coord_position!r}. "
+                "Must be 'left', 'center', or 'right'."
+            )
+
         crs = self._get_crs(dataset)
         x_names = {"lon", "longitude", "x"}
         y_names = {"lat", "latitude", "y"}
         dimensions: dict[str, dict] = {}
+
         for name, coord in dataset.coords.items():
             if name not in dataset.dims:
                 # Skip non-dimension coordinates (e.g. scalar or auxiliary coords).
                 continue
+
             name = str(name)
             lname = name.lower()
-            if lname in x_names:
-                dimensions[name] = {
-                    "type": "spatial",
-                    "axis": "x",
-                    "extent": [float(coord.min()), float(coord.max())],
-                    "reference_system": crs.to_epsg(),
-                }
-            elif lname in y_names:
+
+            if lname in x_names or lname in y_names:
+                coord_min = float(coord.min())
+                coord_max = float(coord.max())
+
+                if coord_position == "center":
+                    resolution = float(abs(coord[:2].diff(name))[0])
+                    extent_min = coord_min - resolution / 2
+                    resolution = float(abs(coord[-2:].diff(name))[0])
+                    extent_max = coord_max + resolution / 2
+                elif coord_position == "left":
+                    extent_min = coord_min
+                    resolution = float(abs(coord[-2:].diff(name))[0])
+                    extent_max = coord_max + resolution
+                else:  # coord_position == "right"
+                    resolution = float(abs(coord[:2].diff(name))[0])
+                    extent_min = coord_min - resolution
+                    extent_max = coord_max
+
                 dimensions[name] = {
                     "type": "spatial",
-                    "axis": "y",
-                    "extent": [float(coord.min()), float(coord.max())],
+                    "axis": "x" if lname in x_names else "y",
+                    "extent": [extent_min, extent_max],
                     "reference_system": crs.to_epsg(),
                 }
+
             elif lname == "time":
                 time_min = pd.to_datetime(coord.min().values).to_pydatetime()
                 time_max = pd.to_datetime(coord.max().values).to_pydatetime()
                 dimensions[name] = {
                     "type": "temporal",
                     "extent": [time_min.isoformat(), time_max.isoformat()],
                 }
+
             else:
                 try:
                     dimensions[name] = {
@@ -843,6 +917,7 @@ def _get_cube_dimensions(self, dataset: xr.Dataset) -> dict[str, dict]:
                         "type": lname,
                         "values": [str(v) for v in coord.values.tolist()],
                     }
+
         return dimensions
 
     @staticmethod
@@ -880,7 +955,7 @@ def build_prr_stac_item(self, item_config: ItemConfig) -> Item:
             f"for collection '{self.collection_id}'."
         )
         dataset = open_dataset(item_config.dataset_id, logger=self.logger)
-        spatial_extent = self._get_spatial_extent(dataset)
+        spatial_extent = self._get_spatial_extent(dataset, self.coord_position)
         temporal_extent = self._get_temporal_extent(dataset)
         general_metadata = self._get_general_metadata(dataset)
 
@@ -918,7 +993,9 @@ def build_prr_stac_item(self, item_config: ItemConfig) -> Item:
                 "description": general_metadata.get("description", ""),
                 "created": now_iso,
                 "updated": now_iso,
-                "cube:dimensions": self._get_cube_dimensions(dataset),
+                "cube:dimensions": self._get_cube_dimensions(
+                    dataset, self.coord_position
+                ),
                 "cube:variables": self._get_cube_variables(dataset),
             },
         )
@@ -1160,7 +1237,7 @@ def build_dataset_stac_collection(
         try:
             assert len(self.items_config) == 1
             dataset = open_dataset(self.items_config[0].dataset_id, logger=self.logger)
-            spatial_extent = self._get_spatial_extent(dataset)
+            spatial_extent = self._get_spatial_extent(dataset, self.coord_position)
             temporal_extent = self._get_temporal_extent(dataset)
             variables = self.get_variable_ids(dataset)
             general_metadata = self._get_general_metadata(dataset)