154-s1orbitselection

src/openeo_gfmap/utils/__init__.py

@@ -1,5 +1,7 @@
 """This sub-module contains utilitary function and tools for OpenEO-GFMap"""
 
+import logging
+
 from openeo_gfmap.utils.build_df import load_json
 from openeo_gfmap.utils.intervals import quintad_intervals
 from openeo_gfmap.utils.netcdf import update_nc_attributes
@@ -12,6 +14,18 @@
     select_sar_bands,
 )
 
+_log = logging.getLogger(__name__)
+_log.setLevel(logging.INFO)
+
+ch = logging.StreamHandler()
+ch.setLevel(logging.INFO)
+
+formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
+ch.setFormatter(formatter)
+
+_log.addHandler(ch)
+
+
 __all__ = [
     "load_json",
     "normalize_array",

src/openeo_gfmap/utils/catalogue.py

@@ -3,6 +3,7 @@
 from typing import Optional
 
 import geojson
+import pandas as pd
 import requests
 from pyproj.crs import CRS
 from rasterio.warp import transform_bounds
@@ -17,6 +18,7 @@
     SpatialContext,
     TemporalContext,
 )
+from openeo_gfmap.utils import _log
 
 request_sessions: Optional[requests.Session] = None
 
@@ -41,13 +43,15 @@ class UncoveredS1Exception(Exception):
 
 
 def _parse_cdse_products(response: dict):
-    """Parses the geometry of products from the CDSE catalogue."""
-    geoemetries = []
+    """Parses the geometry and timestamps of products from the CDSE catalogue."""
+    geometries = []
+    timestamps = []
     products = response["features"]
 
     for product in products:
-        geoemetries.append(shape(product["geometry"]))
-    return geoemetries
+        geometries.append(shape(product["geometry"]))
+        timestamps.append(pd.to_datetime(product["properties"]["startDate"]))
+    return geometries, timestamps
 
 
 def _query_cdse_catalogue(
@@ -139,8 +143,8 @@ def s1_area_per_orbitstate_vvvh(
     temporal_extent: TemporalContext,
 ) -> dict:
     """
-    Evaluates for both the ascending and descending state orbits the area of interesection for the
-    available products with a VV&VH polarisation.
+    Evaluates for both the ascending and descending state orbits the area of interesection and
+    maximum temporal gap for the available products with a VV&VH polarisation.
 
     Parameters
     ----------
@@ -155,8 +159,8 @@ def s1_area_per_orbitstate_vvvh(
     Returns
     ------
     dict
-        Keys containing the orbit state and values containing the total area of intersection in
-        km^2
+        Keys containing the orbit state and values containing the total area of intersection and
+        in km^2 and maximum temporal gap in days.
     """
     if isinstance(spatial_extent, geojson.FeatureCollection):
         # Transform geojson into shapely geometry and compute bounds
@@ -211,26 +215,46 @@ def s1_area_per_orbitstate_vvvh(
     spatial_extent = box(*bounds)
 
     # Computes if there is the full overlap for each of those states
-    union_ascending = unary_union(ascending_products)
-    union_descending = unary_union(descending_products)
+    union_ascending = unary_union(ascending_products[0])
+    union_descending = unary_union(descending_products[0])
 
     ascending_covers = union_ascending.contains(spatial_extent)
     descending_covers = union_descending.contains(spatial_extent)
 
+    # Computes max temporal gap. Include requested start and end date so we dont miss
+    # any start or end gap before first/last observation
+    ascending_timestamps = pd.DatetimeIndex(
+        sorted(
+            [pd.to_datetime(temporal_extent.start_date, utc=True)]
+            + ascending_products[1]
+            + [pd.to_datetime(temporal_extent.end_date, utc=True)]
+        )
+    )
+
+    descending_timestamps = pd.DatetimeIndex(
+        sorted(
+            [pd.to_datetime(temporal_extent.start_date, utc=True)]
+            + descending_products[1]
+            + [pd.to_datetime(temporal_extent.end_date, utc=True)]
+        )
+    )
+
     # Computes the area of intersection
     return {
         "ASCENDING": {
             "full_overlap": ascending_covers,
+            "max_temporal_gap": ascending_timestamps.diff().max().days,
             "area": sum(
                 product.intersection(spatial_extent).area
-                for product in ascending_products
+                for product in ascending_products[0]
             ),
         },
         "DESCENDING": {
             "full_overlap": descending_covers,
+            "max_temporal_gap": descending_timestamps.diff().max().days,
             "area": sum(
                 product.intersection(spatial_extent).area
-                for product in descending_products
+                for product in descending_products[0]
             ),
         },
     }
@@ -240,9 +264,15 @@ def select_s1_orbitstate_vvvh(
     backend: BackendContext,
     spatial_extent: SpatialContext,
     temporal_extent: TemporalContext,
+    max_temporal_gap: int = 30,
 ) -> str:
-    """Selects the orbit state that covers the most area of intersection for the
-    available products with a VV&VH polarisation.
+    """Selects the orbit state based on some predefined rules that
+    are checked in sequential order:
+    1. prefer an orbit with full coverage over the requested bounds
+    2. prefer an orbit with a maximum temporal gap under a
+        predefined threshold
+    3. prefer the orbit that covers the most area of intersection
+        for the available products
 
     Parameters
     ----------
@@ -253,6 +283,8 @@ def select_s1_orbitstate_vvvh(
         The spatial extent to be checked, it will check within its bounding box.
     temporal_extent : TemporalContext
         The temporal period to be checked.
+    max_temporal_gap: int, optional, default: 30
+        The maximum temporal gap in days to be considered for the orbit state.
 
     Returns
     ------
@@ -265,21 +297,55 @@ def select_s1_orbitstate_vvvh(
 
     ascending_overlap = areas["ASCENDING"]["full_overlap"]
     descending_overlap = areas["DESCENDING"]["full_overlap"]
+    ascending_gap_too_large = areas["ASCENDING"]["max_temporal_gap"] > max_temporal_gap
+    descending_gap_too_large = (
+        areas["DESCENDING"]["max_temporal_gap"] > max_temporal_gap
+    )
+
+    orbit_choice = None
 
+    # Rule 1: Prefer an orbit with full coverage over the requested bounds
     if ascending_overlap and not descending_overlap:
-        return "ASCENDING"
+        orbit_choice = "ASCENDING"
+        reason = "Only orbit fully covering the requested area."
     elif descending_overlap and not ascending_overlap:
-        return "DESCENDING"
+        orbit_choice = "DESCENDING"
+        reason = "Only orbit fully covering the requested area."
+    elif ascending_gap_too_large and not descending_gap_too_large:
+        orbit_choice = "DESCENDING"
+        reason = (
+            "Only orbit with temporal gap under the threshold. "
+            f"{areas['DESCENDING']['max_temporal_gap']} days < {max_temporal_gap} days"
+        )
+    elif descending_gap_too_large and not ascending_gap_too_large:
+        orbit_choice = "ASCENDING"
+        reason = (
+            "Only orbit with temporal gap under the threshold. "
+            f"{areas['ASCENDING']['max_temporal_gap']} < {max_temporal_gap}"
+        )
+    # Rule 3: Prefer the orbit that covers the most area of intersection
+    # for the available products
     elif ascending_overlap and descending_overlap:
         ascending_cover_area = areas["ASCENDING"]["area"]
         descending_cover_area = areas["DESCENDING"]["area"]
 
         # Selects the orbit state that covers the most area
         if ascending_cover_area > descending_cover_area:
-            return "ASCENDING"
+            orbit_choice = "ASCENDING"
+            reason = (
+                "Orbit has more cumulative intersected area. "
+                f"{ascending_cover_area} > {descending_cover_area}"
+            )
         else:
-            return "DESCENDING"
+            reason = (
+                "Orbit has more cumulative intersected area. "
+                f"{descending_cover_area} > {ascending_cover_area}"
+            )
+            orbit_choice = "DESCENDING"
 
+    if orbit_choice is not None:
+        _log.info(f"Selected orbit state: {orbit_choice}. Reason: {reason}")
+        return orbit_choice
     raise UncoveredS1Exception(
         "No product available to fully cover the given spatio-temporal context."
     )

tests/test_openeo_gfmap/test_utils.py

@@ -1,6 +1,7 @@
 import os
 from pathlib import Path
 
+import geojson
 import pystac
 import pytest
 from netCDF4 import Dataset
@@ -44,6 +45,9 @@ def test_query_cdse_catalogue():
     assert response["ASCENDING"]["full_overlap"] is True
     assert response["DESCENDING"]["full_overlap"] is True
 
+    assert response["ASCENDING"]["max_temporal_gap"] > 0.0
+    assert response["DESCENDING"]["max_temporal_gap"] > 0.0
+
     # Testing the decision maker, it should return DESCENDING
     decision = select_s1_orbitstate_vvvh(
         backend=backend_context,
@@ -54,6 +58,29 @@ def test_query_cdse_catalogue():
     assert decision == "DESCENDING"
 
 
+def test_query_cdse_catalogue_with_s1_gap():
+    """This example has a large S1 gap in ASCENDING,
+    so the decision should be DESCENDING
+    """
+    backend_context = BackendContext(Backend.CDSE)
+
+    spatial_extent = geojson.loads(
+        (
+            '{"features": [{"geometry": {"coordinates": [[[35.85799, 49.705688], [35.85799, 49.797363], [36.039682, 49.797363], '
+            '[36.039682, 49.705688], [35.85799, 49.705688]]], "type": "Polygon"}, "id": "0", "properties": '
+            '{"GT_available": true, "extract": 1, "index": 12, "sample_id": "ukraine_sunflower", "tile": '
+            '"36UYA", "valid_time": "2019-05-01", "year": 2019}, "type": "Feature"}], "type": "FeatureCollection"}'
+        )
+    )
+    temporal_extent = TemporalContext("2019-01-30", "2019-08-31")
+
+    decision = select_s1_orbitstate_vvvh(
+        backend_context, spatial_extent, temporal_extent
+    )
+
+    assert decision == "DESCENDING"
+
+
 @pytest.fixture
 def temp_nc_file():
     temp_file = Path("temp_test.nc")
@@ -176,3 +203,4 @@ def test_split_collection_by_epsg(tmp_path):
         collection.add_item(missing_epsg_item)
         collection.normalize_and_save(input_dir)
         split_collection_by_epsg(path=input_dir, output_dir=output_dir)
+        split_collection_by_epsg(path=input_dir, output_dir=output_dir)