Add query_proba_results: implement a function to query PROBA result…

src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py

@@ -1,11 +1,9 @@
-import os, sys
 import pandas as pd
 import numpy as np
 import xarray as xr
 import re
-from typing import Dict, List, Tuple, Union
+from typing import Dict, List
 from openeo.udf import inspect
-from datetime import datetime
 from openeo.metadata import CubeMetadata
 
 def apply_metadata(metadata: CubeMetadata, context:Dict) -> CubeMetadata:
@@ -14,7 +12,10 @@ def apply_metadata(metadata: CubeMetadata, context:Dict) -> CubeMetadata:
     :param context: Context of the UDF
     :return: renamed labels
     """
-    return metadata.rename_labels(dimension="bands", target=[f'{context.get('typology','EUNIS')} habitat level3'])
+    band_name = context.get('typology', 'EUNIS')
+    final_band_name = f"{band_name} habitat level3"
+
+    return metadata.rename_labels(dimension="bands", target=[final_band_name])
 
 def _select_highest_prob_class(cube: xr.DataArray, raster_codes) -> xr.DataArray:
     """ Select per model the highest probability of occurrence class
@@ -150,8 +151,9 @@ def apply_datacube(cube: xr.DataArray, context:Dict) -> xr.DataArray:
 
     ### get the list of classes as output from inference run
     # use returned metadata to build up the class dictionary
-    inspect(message=cube.indexes["bands"].values)
-    df = parse_prob_classes_fromStac(cube.indexes["bands"].values)
+    input_band_names = cube.indexes["bands"].values
+    inspect(message=f"input cube band names ({len(input_band_names)}): {input_band_names}")
+    df = parse_prob_classes_fromStac(input_band_names)
 
     inspect(message=f"## context parameters")
     inspect(message=f"{df}")

src/eo_processing/utils/stac_helper.py

@@ -1,10 +1,12 @@
 import pystac_client
+import os
 import logging
 from urllib.request import urlopen
 from io import BytesIO
 from typing import Optional, List
 import geopandas as gpd
 import pandas as pd
+import re
 
 def get_stac_collection_url(collection_id: str, catalog_url: str = "https://catalogue.weed.apex.esa.int/") -> str:
     """
@@ -96,6 +98,132 @@ def query_modelID_output_bands(model_id: str,
 
     return item_collection.items[0].properties["output_band_names"]
 
+def query_proba_results(df_AOI: gpd.GeoDataFrame, collection_id:str, processing_year:int,
+                         stac_url:str = 'https://catalogue.weed.apex.esa.int',
+                         info_debug: bool = True, postprocess: bool = True) -> gpd.GeoDataFrame:
+    """
+    Queries and retrieves PROBA results intersecting with a given Area of Interest (AOI) from a STAC catalog.
+
+    This function searches for PROBA results within the bounding box of the provided AOI and retrieves metadata
+    and asset information from the specified STAC catalog. The retrieved data is reformatted and returned
+    as a GeoDataFrame containing details about the intersecting PROBA tiles.
+
+    Arguments:
+    :param df_AOI: A GeoDataFrame representing the Area of Interest (AOI). The GeoDataFrame should contain geometry
+        information and a coordinate reference system. If the coordinate reference system is not EPSG:4326,
+        the function will reproject it to EPSG:4326.
+    :param collection_id: A string specifying the collection identifier to search within the STAC catalog.
+    :param processing_year: An integer specifying the year for which to retrieve PROBA results.
+    :param stac_url: A string specifying the URL of the STAC catalog to query. Defaults to 'https://catalogue.weed.apex.esa.int'.
+    :param info_debug: A boolean flag to enable or disable debug logging information. Defaults to True.
+    :param postprocess: A boolean flag to determine whether to postprocess the retrieved data. Defaults to True.
+
+    Returns:
+    A GeoDataFrame containing metadata and details about the PROBA results that intersect with the AOI. The GeoDataFrame
+    includes additional columns extracted from the metadata and asset information of the intersecting tiles, such as
+    datetime, bounding box, tile ID, and others.
+
+    Raises:
+    ValueError: If no intersecting PROBA tiles are found in the STAC catalog for the specified collection.
+    """
+    if info_debug: print(f"get_modelID_asset_geometry_from_STAC")
+    # convert AOI into BBOX in 4326
+    if df_AOI.crs != "EPSG:4326":
+        df_AOI_4326 = df_AOI.to_crs("EPSG:4326")
+    else:
+        df_AOI_4326 = df_AOI.copy()
+
+    bbox_4326 = df_AOI_4326.total_bounds
+
+    if info_debug: print(f"- searching for PROBA results in {bbox_4326}")
+    if info_debug: print(f"- using STAC url: {stac_url}")
+    if info_debug: print(f"- using collection id: {collection_id}")
+
+    client = pystac_client.Client.open(stac_url)
+
+    search = client.search(
+        collections=[collection_id],
+        bbox=bbox_4326,
+        fields=["properties", "assets.openEO.href"],
+    )
+
+    results = []
+    for item in search.items_as_dicts():
+        results.append(
+            [item['properties']['datetime'], item['properties']['proj:bbox'], item['properties']['proj:shape'],
+             item['properties']['proj:code'], item['assets']['openEO']['href']])
+
+    # build dataframe
+    df_result = pd.DataFrame(results, columns=['datetime', 'file_bbox', 'file_shape', 'file_epsg', 'file_url'])
+    if info_debug: print(f"- found {len(df_result)} intersecting PROBA tiles")
+    # check if there are any results
+    if df_result.empty:
+        ValueError(f"No intersecting PROBA tiles found in the STAC ({collection_id}).")
+
+    if postprocess:
+        if info_debug: print(f"- postprocessing PROBA tiles")
+        # split out from file_url important parts (file_name, tile_id, etc)
+        df_result['basename'] = df_result['file_url'].apply(lambda x: os.path.basename(x))
+        df_result[['project_typology', 'type', 'processing_year', 'tileID', 'model_short', 'inference_run_version',
+                   'procesisng_start']] = df_result['basename'].str.split('_', expand=True)
+        df_result['processing_year'] = df_result['processing_year'].str[-4:].astype(int)
+
+        # first limit results to processing year
+        df_result = df_result[df_result['processing_year'] == processing_year]
+
+        # check if we have tiles smaller than our standard 20x20km grid - yes then make sure tile name is correct
+        def extract_real_tileid(tileid_variant):
+            """
+            Extract the real tileID by removing trailing letter suffixes.
+            Ensures the tileID ends with a number.
+            Example: '48πXH34a' -> '48πXH34'
+            """
+            # Remove any trailing letters after the last digit
+            return re.sub(r'[a-zA-Z]+$', '', tileid_variant)
+
+        for idx, row in df_result.iterrows():
+            if row.file_shape != [2000, 2000]:
+                df_result.at[idx, 'tileID'] = extract_real_tileid(row.tileID)
+
+        # now we can filter out spatial duplicates for same used modelID_short name
+        # NOTE: that assumes that NEVER different inference runs of smae modelID were saved in same STAC catalog
+        df_result = df_result.drop_duplicates(subset=['tileID', 'model_short'], keep='first')
+
+        # last step. we have to prepare the output file_name.
+        # Step 1: Check if we have duplicate tileIDs with different model_short values
+        duplicate_tiles = df_result.groupby('tileID')['model_short'].apply(lambda x: list(x.unique())).to_dict()
+        tiles_with_multiple_models = {k: v for k, v in duplicate_tiles.items() if len(v) > 1}
+
+        if tiles_with_multiple_models:
+            if info_debug: print(f" -- Found {len(tiles_with_multiple_models)} tiles with multiple model_short values")
+
+            # Step 2: For tiles with multiple models, condense model_short names
+            # Create a condensed model_short by combining unique values
+            for tile, models in tiles_with_multiple_models.items():
+                # Sort models to ensure consistent naming
+                strata = [x.split('-')[0] for x in models]
+                condensed_name = '-'.join(sorted(strata)) + '-' + '-'.join(models[0].split('-')[1:])
+                # Update all rows for this tileID with the condensed name
+                df_result.loc[df_result['tileID'] == tile, 'model_short'] = condensed_name
+
+            # Step 3: Now remove duplicate tileIDs (keeping first occurrence)
+            df_result = df_result.drop_duplicates(subset=['tileID'], keep='first')
+        else:
+            if info_debug: print(" -- No duplicate tileIDs found with different model_short values")
+            # Still remove any exact duplicates
+            df_result = df_result.drop_duplicates(subset=['tileID'], keep='first')
+
+        # Step 4: Create the file_prefix column properly
+        df_result['file_prefix'] = df_result.apply(
+            lambda
+                row: f"{row['project_typology']}_mece-cube_year{row['processing_year']}_{row['tileID']}_{row['model_short']}_{row['inference_run_version']}",
+            axis=1
+        )
+        # filter to final needed
+        df_result = df_result[['tileID', 'file_prefix']]
+
+    return df_result
+
 def get_modelID_asset_geometry_from_STAC(df_AOI: gpd.GeoDataFrame, typology_schema: str = 'IUCNGET',
                                      model_version: Optional[str]=None,
                                      stac_url:str = 'https://catalogue.weed.apex.esa.int',