diff --git a/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py b/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py index 550a57c..3348772 100644 --- a/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py +++ b/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}") diff --git a/src/eo_processing/utils/stac_helper.py b/src/eo_processing/utils/stac_helper.py index 046a87c..4f9ac1e 100644 --- a/src/eo_processing/utils/stac_helper.py +++ b/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',