Region code generation should match postgis logic

Author: jeremyh

cubedash/summary/_extents.py

@@ -36,7 +36,7 @@
 from datacube.drivers.postgres._fields import PgDocField, RangeDocField
 
 from datacube.index import Index
-from datacube.model import DatasetType, Field, MetadataType
+from datacube.model import DatasetType, Field, MetadataType, Dataset
 
 _LOG = structlog.get_logger()
 
@@ -527,7 +527,9 @@ def __init__(
     units_label: str = "regions"
 
     @classmethod
-    def for_product(cls, dt: DatasetType, region_shapes: Dict[str, GeometryCollection]):
+    def for_product(
+        cls, dt: DatasetType, region_shapes: Dict[str, GeometryCollection] = None
+    ):
         region_code_field: Field = dt.metadata_type.dataset_fields.get("region_code")
         grid_spec = dt.grid_spec
         # Ingested grids trump the "region_code" field because they've probably sliced it up smaller.
@@ -557,7 +559,23 @@ def geojson_extent(self, region_code):
             "properties": {"region_code": region_code},
         }
 
+    def dataset_region_code(self, dataset: Dataset) -> Optional[str]:
+        """
+        Get the region code for a dataset.
+
+        This should always give the same result as the alchemy_expression() function,
+        but is computed in pure python.
+
+        Classes that override alchemy_expression should override this to match.
+        """
+        return dataset.metadata.region_code
+
     def alchemy_expression(self):
+        """
+        Get an alchemy expression that computes dataset's region code
+
+        Classes that override this should also override dataset_region_code to match.
+        """
         dt = self.product
         region_code_field: Field = dt.metadata_type.dataset_fields.get("region_code")
         # `alchemy_expression` is part of the postgres driver (PgDocField),
@@ -625,6 +643,21 @@ def alchemy_expression(self):
             func.floor((func.ST_Y(center_point) - origin_y) / size_y).cast(String),
         )
 
+    def dataset_region_code(self, dataset: Dataset) -> Optional[str]:
+        tiles = [
+            tile
+            for tile, _ in dataset.type.grid_spec.tiles(
+                dataset.extent.centroid.boundingbox
+            )
+        ]
+        if not len(tiles) == 1:
+            raise ValueError(
+                "Tiled dataset should only have one tile? "
+                f"Got {tiles!r} for {dataset!r}"
+            )
+        x, y = tiles[0]
+        return f"{x}_{y}"
+
 
 def _from_xy_region_code(region_code: str):
     """
@@ -676,6 +709,19 @@ def alchemy_expression(self):
             else_=path_field.lower.alchemy_expression.cast(String),
         )
 
+    def dataset_region_code(self, dataset: Dataset) -> Optional[str]:
+        path_range = dataset.metadata.fields["sat_path"]
+        row_range = dataset.metadata.fields["sat_row"]
+        if row_range is None and path_range is None:
+            return None
+
+        # If it's just one scene? Include it specifically
+        if row_range[0] == row_range[1]:
+            return f"{path_range[0]}_{row_range[1]}"
+        # Otherwise it's a range of rows, so we say the whole path.
+        else:
+            return f"{path_range[0]}"
+
 
 def _region_code_field(dt: DatasetType):
     """

cubedash/warmup.py

@@ -7,6 +7,7 @@
 import click
 from click import secho
 
+from cubedash.summary import RegionInfo
 from datacube.index import Index
 from datacube.ui.click import config_option, environment_option, pass_index
 
@@ -48,11 +49,15 @@ def find_examples_of_all_public_urls(index: Index):
             # TODO: Do non-region_code regions too (such as ingested data)
             # TODO: Actually we have no EO3 in this test data, so it does nothing.
             #       Maybe add test data from test_eo3_support.py?
-            if "region_code" in dataset.metadata.fields:
-                yield f"/api/regions/{name}/{time:%Y/%m/%d}"
+            region_info = RegionInfo.for_product(dt)
 
-                yield f"/region/{name}/{dataset.metadata.region_code}"
-                yield f"/region/{name}/{dataset.metadata.region_code}/{time:%Y/%m/%d}"
+            if region_info is not None:
+                region_code = region_info.dataset_region_code(dataset)
+                if region_code is not None:
+                    yield f"/api/regions/{name}/{time:%Y/%m/%d}"
+
+                    yield f"/region/{name}/{region_code}"
+                    yield f"/region/{name}/{region_code}/{time:%Y/%m/%d}"
 
     for [dataset_id] in index.datasets.search_returning(("id",), limit=10):
         yield f"/dataset/{dataset_id}"

integration_tests/test_summarise_data.py

@@ -12,6 +12,7 @@
 
 from cubedash._utils import alchemy_engine
 from cubedash.summary import SummaryStore
+from cubedash.summary._extents import GridRegionInfo
 from cubedash.summary._schema import CUBEDASH_SCHEMA
 from datacube.index import Index
 from datacube.index.hl import Doc2Dataset
@@ -427,3 +428,28 @@ def _get_product_seq_value():
     assert (
         value_after_rerun == original_value
     ), "Product sequence was incremented without any new products being added."
+
+
+def test_computed_regions_match_those_summarised(summary_store: SummaryStore):
+    """
+    The region code for all datasets should be computed identically when
+    done in both SQL and Python.
+    """
+    summary_store.refresh_all_products()
+
+    # Loop through all datasets in the test data to check that the the DB and Python
+    # functions give identical region codes.
+    for product in summary_store.index.products.get_all():
+        region_info = GridRegionInfo.for_product(product, None)
+        for dataset in summary_store.index.datasets.search(product=product.name):
+            (
+                footprint,
+                alchemy_calculated_region_code,
+            ) = summary_store.get_dataset_footprint_region(dataset.id)
+
+            python_calculated_region_code = region_info.dataset_region_code(dataset)
+            assert python_calculated_region_code == alchemy_calculated_region_code, (
+                "Python and DB calculated region codes didn't product the same value. "
+                f"{python_calculated_region_code!r} != {alchemy_calculated_region_code!r}"
+                f"for product {dataset.type.name!r}, dataset {dataset!r}"
+            )