diff --git a/algorithm_catalog/terravision/s2_sam/benchmark_scenarios/s2_sam.json b/algorithm_catalog/terravision/s2_sam/benchmark_scenarios/s2_sam.json new file mode 100644 index 00000000..856a4cc6 --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/benchmark_scenarios/s2_sam.json @@ -0,0 +1,85 @@ +[ + { + "id": "s2_sam", + "type": "openeo", + "description": "Run Sentinel-2 based Spectral Angle Mapping (S2-SAM) on a monthly BAP composite and return the best matching USGS spectral library class index per pixel.", + "backend": "openeofed.dataspace.copernicus.eu", + "process_graph": { + "s2sam1": { + "arguments": { + "spatial_extent": { + "geometry": { + "coordinates": [ + [ + [ + -15.487987, + 27.813076 + ], + [ + -15.489753, + 27.814047 + ], + [ + -15.493856, + 27.809741 + ], + [ + -15.492555, + 27.808545 + ], + [ + -15.487809, + 27.806311 + ], + [ + -15.485954, + 27.809176 + ], + [ + -15.484348, + 27.812444 + ], + [ + -15.48635, + 27.81229 + ], + [ + -15.487987, + 27.813076 + ] + ] + ], + "type": "Polygon" + }, + "properties": { + "CNTR_ID": "ES", + "ID": "598" + }, + "type": "Feature" + }, + "temporal_extent": [ + "2024-05-01", + "2024-05-31" + ] + }, + "namespace": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/c0f8d549a866737a3c276f9838fa0a019e7975be/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json", + "process_id": "s2_sam" + }, + "saveresult1": { + "arguments": { + "data": { + "from_node": "s2sam1" + }, + "format": "NetCDF", + "options": {} + }, + "process_id": "save_result", + "result": true + } + }, + "reference_data": { + "job-results.json": "https://s3.waw3-1.cloudferro.com/apex-benchmarks-internal/HP120!tests_test_benchmarks.py__test_run_benchmark_s2_sam_!actual/job-results.json", + "openEO.nc": "https://s3.waw3-1.cloudferro.com/apex-benchmarks-internal/HP120!tests_test_benchmarks.py__test_run_benchmark_s2_sam_!actual/openEO.nc" + } + } +] \ No newline at end of file diff --git a/algorithm_catalog/terravision/s2_sam/notebook/s2_sam.ipynb b/algorithm_catalog/terravision/s2_sam/notebook/s2_sam.ipynb new file mode 100644 index 00000000..6344366f --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/notebook/s2_sam.ipynb @@ -0,0 +1,156 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "31162cad", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Authenticated using refresh token.\n" + ] + } + ], + "source": [ + "import openeo\n", + "\n", + "connection = openeo.connect(\"openeofed.dataspace.copernicus.eu\").authenticate_oidc()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "26bc37db", + "metadata": {}, + "outputs": [], + "source": [ + "spatial_extent = {\n", + " \"type\": \"Feature\",\n", + " \"properties\": {\n", + " \"ID\": \"598\",\n", + " \"CNTR_ID\": \"ES\"\n", + " },\n", + " \"geometry\": {\n", + " \"type\": \"Polygon\",\n", + " \"coordinates\": [[\n", + " [-15.487987, 27.813076],\n", + " [-15.489753, 27.814047],\n", + " [-15.493856, 27.809741],\n", + " [-15.492555, 27.808545],\n", + " [-15.487809, 27.806311],\n", + " [-15.485954, 27.809176],\n", + " [-15.484348, 27.812444],\n", + " [-15.486350, 27.812290],\n", + " [-15.487987, 27.813076]\n", + " ]]\n", + " }\n", + "}\n", + "temporal_extent = [\"2024-05-01\", \"2024-05-31\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "9fb09937", + "metadata": {}, + "outputs": [], + "source": [ + "public_url = \"https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/c0f8d549a866737a3c276f9838fa0a019e7975be/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json\"\n", + "sam_cube = connection.datacube_from_process(\n", + " \"s2_sam\", \n", + " namespace = public_url,\n", + " temporal_extent = temporal_extent,\n", + " spatial_extent= spatial_extent\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "4918618a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0:00:00 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': send 'start'\n", + "0:00:05 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:00:10 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:00:17 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:00:26 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:00:36 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:00:48 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': queued (progress 0%)\n", + "0:01:04 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': running (progress N/A)\n", + "0:01:23 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': running (progress N/A)\n", + "0:01:47 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': running (progress N/A)\n", + "0:02:17 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': running (progress N/A)\n", + "0:02:54 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': running (progress N/A)\n", + "0:03:41 Job 'cdse-j-2608111342194bf0ac97d86c30fb35e5': finished (progress 100%)\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " " + ], + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "job = sam_cube.create_job(\n", + " title=f'S2 SAM processing for {temporal_extent[0]} to {temporal_extent[1]}',\n", + " out_format=\"NetCDF\")\n", + "job.start_and_wait()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cdse", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/algorithm_catalog/terravision/s2_sam/openeo_udp/README.md b/algorithm_catalog/terravision/s2_sam/openeo_udp/README.md new file mode 100644 index 00000000..dda28d1a --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/openeo_udp/README.md @@ -0,0 +1,50 @@ +# Sentinel-2 based Spectral Angle Mapping (S2-SAM) + +This service applies Spectral Angle Mapping (SAM) to Sentinel-2 Level-2A reflectance data. +For each pixel, the algorithm compares the pixel spectrum to a spectral reference library and selects the closest match based on the smallest spectral angle. + +## What this service does + +SAM measures similarity in spectral *shape* (not absolute brightness): + +$$ +θ = \cos^{-1}\left(\frac{x \cdot y}{\|x\|\|y\|}\right) +$$ + +where: +- $x$ is the pixel spectrum. +- $y$ is a reference spectrum from the library. +- smaller $\theta$ means a better match. + +The output is a raster of class indices, where each pixel value is the index of the best-matching reference signature. + +## Sentinel-2 bands used + +The service uses only the Sentinel-2 L2A bands provided by the BAP composite input, not all L2A layers. +The band subset is (B01, B02, B03, B04, B05, B06, B07, B08, B8A, B09, B11, B12). + +Input compositing is based on the BAP process: +[https://algorithm-catalogue.apex.esa.int/apps/bap_composite](https://algorithm-catalogue.apex.esa.int/apps/bap_composite) + +## Reference spectral library + +The reference library is based on a Sentinel-2-resampled USGS mineral/material collection. +It includes a broad set of geology and mining-related targets (e.g. iron oxides, clays, sulfates, silicates, carbonates, and related materials). + +The CSV used by the service contains: +- 498 named reference signatures (material columns). +- 1 `cwvl` column with central wavelength values. + +Examples of included targets are minerals such as Acmite, Actinolite, Alunite, Jarosite, Kaolinite, Hematite, Goethite, Pyrite, Chalcopyrite, and many more variants. + +## Notes on interpretation + +- The product reports the **closest spectral match** in the selected library, not a definitive mineral identification. +- Similar materials can produce similar spectra at Sentinel-2 spectral resolution. +- Results are most useful when combined with geological context, ancillary data, and expert validation. + +## Missing data handling + +If some bands are missing (`NaN`) at a pixel, the implementation computes SAM using only the available bands, with a minimum valid-band threshold. +Pixels without enough valid information are marked as no match. + diff --git a/algorithm_catalog/terravision/s2_sam/openeo_udp/generate.py b/algorithm_catalog/terravision/s2_sam/openeo_udp/generate.py new file mode 100644 index 00000000..479b8390 --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/openeo_udp/generate.py @@ -0,0 +1,52 @@ +from datetime import date +import json +from pathlib import Path + +import openeo +from openeo.api.process import Parameter +from openeo.rest.udp import build_process_dict + +def generate(): + connection = openeo.connect("openeofed.dataspace.copernicus.eu").authenticate_oidc() + + spatial_extent = Parameter.spatial_extent( + name="spatial_extent", + description="Limits the data to process to the specified bounding box or polygons.\\n\\nFor raster data, the process loads the pixel into the data cube if the point at the pixel center intersects with the bounding box or any of the polygons (as defined in the Simple Features standard by the OGC).\\nFor vector data, the process loads the geometry into the data cube if the geometry is fully within the bounding box or any of the polygons (as defined in the Simple Features standard by the OGC). Empty geometries may only be in the data cube if no spatial extent has been provided.\\n\\nEmpty geometries are ignored.\\nSet this parameter to null to set no limit for the spatial extent." + ) + + temporal_extent = Parameter.temporal_interval( + name="temporal_extent", + description="Temporal extent specified as two-element array with start and end date/date-time." + ) + + + bap_cube = connection.datacube_from_process( + "bap_composite", + namespace = "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/refs/heads/main/algorithm_catalog/vito/bap_composite/openeo_udp/bap_composite.json", + temporal_extent = temporal_extent, + geometry= spatial_extent, + bands = ['B01','B02','B03','B04','B05','B06','B07','B08','B8A','B09','B11','B12'], + max_cloud_cover = 80 + ) + composite = bap_cube.aggregate_temporal_period("month","first") + + udf = openeo.UDF.from_file( + Path(__file__).parent / "sam_udf.py", + ) + processed_cube = composite.apply(process=udf) + + + return build_process_dict( + process_graph=processed_cube, + process_id="s2_sam", + description=(Path(__file__).parent / "README.md").read_text(), + parameters=[ + spatial_extent, + temporal_extent, + ] + ) + + +if __name__ == "__main__": + with open("s2_sam.json", "w") as f: + json.dump(generate(), f, indent=2) \ No newline at end of file diff --git a/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json b/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json new file mode 100644 index 00000000..1702820b --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json @@ -0,0 +1,196 @@ +{ + "process_graph": { + "bapcomposite1": { + "process_id": "bap_composite", + "arguments": { + "bands": [ + "B01", + "B02", + "B03", + "B04", + "B05", + "B06", + "B07", + "B08", + "B8A", + "B09", + "B11", + "B12" + ], + "geometry": { + "from_parameter": "spatial_extent" + }, + "max_cloud_cover": 80, + "temporal_extent": { + "from_parameter": "temporal_extent" + } + }, + "namespace": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/refs/heads/main/algorithm_catalog/vito/bap_composite/openeo_udp/bap_composite.json" + }, + "aggregatetemporalperiod1": { + "process_id": "aggregate_temporal_period", + "arguments": { + "data": { + "from_node": "bapcomposite1" + }, + "period": "month", + "reducer": { + "process_graph": { + "first1": { + "process_id": "first", + "arguments": { + "data": { + "from_parameter": "data" + } + }, + "result": true + } + } + } + } + }, + "apply1": { + "process_id": "apply", + "arguments": { + "data": { + "from_node": "aggregatetemporalperiod1" + }, + "process": { + "process_graph": { + "runudf1": { + "process_id": "run_udf", + "arguments": { + "data": { + "from_parameter": "x" + }, + "runtime": "Python", + "udf": "import xarray\nimport numpy as np\nimport pandas as pd\n\nfrom openeo.udf import inspect\n\n\"\"\"Spectral Angle Mapping UDF for Sentinel-2 based material matching.\n\nThis UDF compares each pixel spectrum against a USGS-derived reference library\nand returns the index of the best matching reference signature.\n\"\"\"\n\n\ndef compute_sam_w_nan(\n img: np.ndarray,\n usgs_df: pd.DataFrame,\n min_bands: int=12,\n) -> np.ndarray:\n \"\"\"Compute SAM\n\n Only overlapping finite bands are used for each pixel-reference pair.\n Pixels with fewer than min_bands valid overlaps are marked as no match.\n \"\"\"\n bands, height, width = img.shape\n # Flatten spatial dimensions so SAM can be computed in matrix form.\n pixel_spectra = img.reshape(bands, -1).T\n # Transpose library so each row corresponds to one reference signature.\n usgs_matrix = usgs_df.T.values\n\n n_pix = pixel_spectra.shape[0]\n n_mat = usgs_matrix.shape[0]\n\n # Pre-fill with +inf so non-computable matches remain easy to detect.\n angles = np.full((n_pix, n_mat), np.inf, dtype=np.float32)\n\n # Compute pixel finite mask once and reuse it for all reference signatures.\n Xfinite = np.isfinite(pixel_spectra)\n\n for j in range(n_mat):\n y = usgs_matrix[j]\n yfinite = np.isfinite(y)\n\n valid = Xfinite & yfinite[None, :]\n n_valid = valid.sum(axis=1)\n\n # Enforce a minimum number of overlapping bands per pixel-match pair.\n ok = n_valid >= min_bands\n if not np.any(ok):\n continue\n\n # Compute dot product and norms only over valid overlap.\n Xv = np.where(valid, pixel_spectra, 0.0)\n yv = np.where(yfinite, y, 0.0)\n\n dot = (Xv * yv[None, :]).sum(axis=1)\n nx = np.sqrt((Xv * Xv).sum(axis=1))\n ny = np.sqrt((yv * yv).sum())\n\n # Add a small epsilon to keep the denominator numerically safe.\n denom = nx * ny + 1e-12\n cosang = np.clip(dot / denom, -1.0, 1.0)\n\n # Convert cosine similarity to spectral angle in radians.\n ang = np.arccos(cosang).astype(np.float32)\n angles[ok, j] = ang[ok]\n\n # Select the reference with minimum angle for each pixel.\n best = np.argmin(angles, axis=1).reshape(height, width)\n # Pixels that never reached the min_bands threshold remain as no match.\n best[np.isinf(angles).all(axis=1).reshape(height, width)] = -1\n return best\n\ndef compute_sam_wo_nan(img: np.ndarray, usgs_df: pd.DataFrame) -> np.array:\n \"\"\"Compute SAM for fully valid arrays without missing values.\"\"\"\n bands, height, width = img.shape\n inspect(data=img.shape, message=\"SAM no-NaN path: input image shape\", level='debug')\n # Flatten to [pixels, bands] to vectorize pairwise comparisons.\n pixel_spectra = img.reshape(bands, -1).T\n\n # Normalize each pixel vector so comparison focuses on spectral shape.\n pixel_norms = np.linalg.norm(pixel_spectra, axis=1, keepdims=True)\n inspect(data=pixel_norms.shape, message=\"SAM no-NaN path: pixel norms shape\", level='debug')\n pixel_spectra_norm = pixel_spectra / (pixel_norms + 1e-10)\n inspect(data=pixel_spectra_norm.shape, message=\"SAM no-NaN path: normalized pixel spectra shape\", level='debug')\n\n # Normalize each reference signature with the same convention.\n inspect(data=usgs_df.shape, message=\"SAM no-NaN path: USGS dataframe shape\", level='debug')\n usgs_matrix = usgs_df.T.values\n inspect(data=usgs_matrix.shape, message=\"SAM no-NaN path: USGS matrix shape\", level='debug')\n usgs_norms = np.linalg.norm(usgs_matrix, axis=1, keepdims=True)\n inspect(data=usgs_norms.shape, message=\"SAM no-NaN path: USGS norms shape\", level='debug')\n usgs_matrix_norm = usgs_matrix / (usgs_norms + 1e-10)\n inspect(data=usgs_matrix_norm.shape, message=\"SAM no-NaN path: normalized USGS matrix shape\", level='debug')\n # Compute cosine similarity between each pixel and each reference spectrum.\n dot_product = np.dot(\n pixel_spectra_norm, usgs_matrix_norm.T\n )\n inspect(data=dot_product.shape, message=\"SAM no-NaN path: dot product shape\", level='debug')\n # Convert similarities to angles for SAM-based ranking.\n angles = np.arccos(np.clip(dot_product, -1.0, 1.0))\n inspect(data=angles.shape, message=\"SAM no-NaN path: angle matrix shape\", level='debug')\n # Lower angle means closer spectral match.\n best_match_indices = np.argmin(angles, axis=1)\n inspect(data=best_match_indices.shape, message=\"SAM no-NaN path: best match index vector shape\", level='debug')\n # Restore the original 2D spatial layout.\n match_index_map = best_match_indices.reshape(height, width)\n inspect(data=match_index_map.shape, message=\"SAM no-NaN path: output map shape\", level='debug')\n return match_index_map\n\n\n\n# UDF entry point\ndef apply_datacube(cube: xarray.DataArray, context: dict) -> xarray.DataArray:\n \"\"\"openEO UDF entry point.\n\n Expects an input cube with band-first layout and returns a 2D class-index map.\n \"\"\"\n inspect(data=cube.dims, message=\"SAM UDF: input cube dimensions\", level='debug')\n inspect(data=cube.shape, message=\"SAM UDF: input cube shape\", level='debug')\n # External reference library used for spectral matching.\n csv_url = \"https://artifactory.vgt.vito.be/artifactory/auxdata-public/openEO_terravision/S2_USGS_library498.csv\"\n usgs_df = pd.read_csv(csv_url, sep=\";\")\n inspect(data=usgs_df.shape, message=\"SAM UDF: loaded USGS dataframe shape\", level='debug')\n\n # Use the NaN-aware path only when needed.\n if np.isnan(cube.values).any():\n result_cube = compute_sam_w_nan(cube.values, usgs_df, 12)\n else:\n result_cube = compute_sam_wo_nan(cube.values, usgs_df)\n\n # Return only spatial dimensions, with class index per pixel.\n return xarray.DataArray(result_cube, dims=cube.dims[1:], coords={dim: cube.coords[dim] for dim in cube.dims[1:]})" + }, + "result": true + } + } + } + }, + "result": true + } + }, + "id": "s2_sam", + "description": "# Sentinel-2 based Spectral Angle Mapping (S2-SAM)\n\nThis service applies Spectral Angle Mapping (SAM) to Sentinel-2 Level-2A reflectance data.\nFor each pixel, the algorithm compares the pixel spectrum to a spectral reference library and selects the closest match based on the smallest spectral angle.\n\n## What this service does\n\nSAM measures similarity in spectral *shape* (not absolute brightness):\n\n$$\nθ = \\cos^{-1}\\left(\\frac{x \\cdot y}{\\|x\\|\\|y\\|}\\right)\n$$\n\nwhere:\n- $x$ is the pixel spectrum.\n- $y$ is a reference spectrum from the library.\n- smaller $\\theta$ means a better match.\n\nThe output is a raster of class indices, where each pixel value is the index of the best-matching reference signature.\n\n## Sentinel-2 bands used\n\nThe service uses only the Sentinel-2 L2A bands provided by the BAP composite input, not all L2A layers.\nThe band subset is (B01, B02, B03, B04, B05, B06, B07, B08, B8A, B09, B11, B12).\n\nInput compositing is based on the BAP process:\n[https://algorithm-catalogue.apex.esa.int/apps/bap_composite](https://algorithm-catalogue.apex.esa.int/apps/bap_composite)\n\n## Reference spectral library\n\nThe reference library is based on a Sentinel-2-resampled USGS mineral/material collection.\nIt includes a broad set of geology and mining-related targets (e.g. iron oxides, clays, sulfates, silicates, carbonates, and related materials).\n\nThe CSV used by the service contains:\n- 498 named reference signatures (material columns).\n- 1 `cwvl` column with central wavelength values.\n\nExamples of included targets are minerals such as Acmite, Actinolite, Alunite, Jarosite, Kaolinite, Hematite, Goethite, Pyrite, Chalcopyrite, and many more variants.\n\n## Notes on interpretation\n\n- The product reports the **closest spectral match** in the selected library, not a definitive mineral identification.\n- Similar materials can produce similar spectra at Sentinel-2 spectral resolution.\n- Results are most useful when combined with geological context, ancillary data, and expert validation.\n\n## Missing data handling\n\nIf some bands are missing (`NaN`) at a pixel, the implementation computes SAM using only the available bands, with a minimum valid-band threshold.\nPixels without enough valid information are marked as no match.\n\n", + "parameters": [ + { + "name": "spatial_extent", + "description": "Limits the data to process to the specified bounding box or polygons.\\n\\nFor raster data, the process loads the pixel into the data cube if the point at the pixel center intersects with the bounding box or any of the polygons (as defined in the Simple Features standard by the OGC).\\nFor vector data, the process loads the geometry into the data cube if the geometry is fully within the bounding box or any of the polygons (as defined in the Simple Features standard by the OGC). Empty geometries may only be in the data cube if no spatial extent has been provided.\\n\\nEmpty geometries are ignored.\\nSet this parameter to null to set no limit for the spatial extent.", + "schema": [ + { + "title": "Bounding Box", + "type": "object", + "subtype": "bounding-box", + "required": [ + "west", + "south", + "east", + "north" + ], + "properties": { + "west": { + "description": "West (lower left corner, coordinate axis 1).", + "type": "number" + }, + "south": { + "description": "South (lower left corner, coordinate axis 2).", + "type": "number" + }, + "east": { + "description": "East (upper right corner, coordinate axis 1).", + "type": "number" + }, + "north": { + "description": "North (upper right corner, coordinate axis 2).", + "type": "number" + }, + "base": { + "description": "Base (optional, lower left corner, coordinate axis 3).", + "type": [ + "number", + "null" + ], + "default": null + }, + "height": { + "description": "Height (optional, upper right corner, coordinate axis 3).", + "type": [ + "number", + "null" + ], + "default": null + }, + "crs": { + "description": "Coordinate reference system of the extent, specified as as [EPSG code](http://www.epsg-registry.org/) or [WKT2 CRS string](http://docs.opengeospatial.org/is/18-010r7/18-010r7.html). Defaults to `4326` (EPSG code 4326) unless the client explicitly requests a different coordinate reference system.", + "anyOf": [ + { + "title": "EPSG Code", + "type": "integer", + "subtype": "epsg-code", + "minimum": 1000, + "examples": [ + 3857 + ] + }, + { + "title": "WKT2", + "type": "string", + "subtype": "wkt2-definition" + } + ], + "default": 4326 + } + } + }, + { + "title": "Vector data cube", + "description": "Limits the data cube to the bounding box of the given geometries in the vector data cube. For raster data, all pixels inside the bounding box that do not intersect with any of the polygons will be set to no data (`null`). Empty geometries are ignored.", + "type": "object", + "subtype": "datacube", + "dimensions": [ + { + "type": "geometry" + } + ] + }, + { + "title": "No filter", + "description": "Don't filter spatially. All data is included in the data cube.", + "type": "null" + } + ] + }, + { + "name": "temporal_extent", + "description": "Temporal extent specified as two-element array with start and end date/date-time.", + "schema": { + "type": "array", + "subtype": "temporal-interval", + "uniqueItems": true, + "minItems": 2, + "maxItems": 2, + "items": { + "anyOf": [ + { + "type": "string", + "subtype": "date-time", + "format": "date-time" + }, + { + "type": "string", + "subtype": "date", + "format": "date" + }, + { + "type": "null" + } + ] + } + } + } + ] +} \ No newline at end of file diff --git a/algorithm_catalog/terravision/s2_sam/openeo_udp/sam_udf.py b/algorithm_catalog/terravision/s2_sam/openeo_udp/sam_udf.py new file mode 100644 index 00000000..131e4525 --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/openeo_udp/sam_udf.py @@ -0,0 +1,132 @@ +import xarray +import numpy as np +import pandas as pd + +from openeo.udf import inspect + +"""Spectral Angle Mapping UDF for Sentinel-2 based material matching. + +This UDF compares each pixel spectrum against a USGS-derived reference library +and returns the index of the best matching reference signature. +""" + + +def compute_sam_w_nan( + img: np.ndarray, + usgs_df: pd.DataFrame, + min_bands: int=12, +) -> np.ndarray: + """Compute SAM + + Only overlapping finite bands are used for each pixel-reference pair. + Pixels with fewer than min_bands valid overlaps are marked as no match. + """ + bands, height, width = img.shape + # Flatten spatial dimensions so SAM can be computed in matrix form. + pixel_spectra = img.reshape(bands, -1).T + # Transpose library so each row corresponds to one reference signature. + usgs_matrix = usgs_df.T.values + + n_pix = pixel_spectra.shape[0] + n_mat = usgs_matrix.shape[0] + + # Pre-fill with +inf so non-computable matches remain easy to detect. + angles = np.full((n_pix, n_mat), np.inf, dtype=np.float32) + + # Compute pixel finite mask once and reuse it for all reference signatures. + Xfinite = np.isfinite(pixel_spectra) + + for j in range(n_mat): + y = usgs_matrix[j] + yfinite = np.isfinite(y) + + valid = Xfinite & yfinite[None, :] + n_valid = valid.sum(axis=1) + + # Enforce a minimum number of overlapping bands per pixel-match pair. + ok = n_valid >= min_bands + if not np.any(ok): + continue + + # Compute dot product and norms only over valid overlap. + Xv = np.where(valid, pixel_spectra, 0.0) + yv = np.where(yfinite, y, 0.0) + + dot = (Xv * yv[None, :]).sum(axis=1) + nx = np.sqrt((Xv * Xv).sum(axis=1)) + ny = np.sqrt((yv * yv).sum()) + + # Add a small epsilon to keep the denominator numerically safe. + denom = nx * ny + 1e-12 + cosang = np.clip(dot / denom, -1.0, 1.0) + + # Convert cosine similarity to spectral angle in radians. + ang = np.arccos(cosang).astype(np.float32) + angles[ok, j] = ang[ok] + + # Select the reference with minimum angle for each pixel. + best = np.argmin(angles, axis=1).reshape(height, width) + # Pixels that never reached the min_bands threshold remain as no match. + best[np.isinf(angles).all(axis=1).reshape(height, width)] = -1 + return best + +def compute_sam_wo_nan(img: np.ndarray, usgs_df: pd.DataFrame) -> np.array: + """Compute SAM for fully valid arrays without missing values.""" + bands, height, width = img.shape + inspect(data=img.shape, message="SAM no-NaN path: input image shape", level='debug') + # Flatten to [pixels, bands] to vectorize pairwise comparisons. + pixel_spectra = img.reshape(bands, -1).T + + # Normalize each pixel vector so comparison focuses on spectral shape. + pixel_norms = np.linalg.norm(pixel_spectra, axis=1, keepdims=True) + inspect(data=pixel_norms.shape, message="SAM no-NaN path: pixel norms shape", level='debug') + pixel_spectra_norm = pixel_spectra / (pixel_norms + 1e-10) + inspect(data=pixel_spectra_norm.shape, message="SAM no-NaN path: normalized pixel spectra shape", level='debug') + + # Normalize each reference signature with the same convention. + inspect(data=usgs_df.shape, message="SAM no-NaN path: USGS dataframe shape", level='debug') + usgs_matrix = usgs_df.T.values + inspect(data=usgs_matrix.shape, message="SAM no-NaN path: USGS matrix shape", level='debug') + usgs_norms = np.linalg.norm(usgs_matrix, axis=1, keepdims=True) + inspect(data=usgs_norms.shape, message="SAM no-NaN path: USGS norms shape", level='debug') + usgs_matrix_norm = usgs_matrix / (usgs_norms + 1e-10) + inspect(data=usgs_matrix_norm.shape, message="SAM no-NaN path: normalized USGS matrix shape", level='debug') + # Compute cosine similarity between each pixel and each reference spectrum. + dot_product = np.dot( + pixel_spectra_norm, usgs_matrix_norm.T + ) + inspect(data=dot_product.shape, message="SAM no-NaN path: dot product shape", level='debug') + # Convert similarities to angles for SAM-based ranking. + angles = np.arccos(np.clip(dot_product, -1.0, 1.0)) + inspect(data=angles.shape, message="SAM no-NaN path: angle matrix shape", level='debug') + # Lower angle means closer spectral match. + best_match_indices = np.argmin(angles, axis=1) + inspect(data=best_match_indices.shape, message="SAM no-NaN path: best match index vector shape", level='debug') + # Restore the original 2D spatial layout. + match_index_map = best_match_indices.reshape(height, width) + inspect(data=match_index_map.shape, message="SAM no-NaN path: output map shape", level='debug') + return match_index_map + + + +# UDF entry point +def apply_datacube(cube: xarray.DataArray, context: dict) -> xarray.DataArray: + """openEO UDF entry point. + + Expects an input cube with band-first layout and returns a 2D class-index map. + """ + inspect(data=cube.dims, message="SAM UDF: input cube dimensions", level='debug') + inspect(data=cube.shape, message="SAM UDF: input cube shape", level='debug') + # External reference library used for spectral matching. + csv_url = "https://artifactory.vgt.vito.be/artifactory/auxdata-public/openEO_terravision/S2_USGS_library498.csv" + usgs_df = pd.read_csv(csv_url, sep=";") + inspect(data=usgs_df.shape, message="SAM UDF: loaded USGS dataframe shape", level='debug') + + # Use the NaN-aware path only when needed. + if np.isnan(cube.values).any(): + result_cube = compute_sam_w_nan(cube.values, usgs_df, 12) + else: + result_cube = compute_sam_wo_nan(cube.values, usgs_df) + + # Return only spatial dimensions, with class index per pixel. + return xarray.DataArray(result_cube, dims=cube.dims[1:], coords={dim: cube.coords[dim] for dim in cube.dims[1:]}) \ No newline at end of file diff --git a/algorithm_catalog/terravision/s2_sam/records/s2_sam.json b/algorithm_catalog/terravision/s2_sam/records/s2_sam.json new file mode 100644 index 00000000..ea511c97 --- /dev/null +++ b/algorithm_catalog/terravision/s2_sam/records/s2_sam.json @@ -0,0 +1,133 @@ +{ + "id": "s2_sam", + "type": "Feature", + "conformsTo": [ + "https://www.opengis.net/spec/ogcapi-records-1/1.0/req/record-core", + "https://apex.esa.int/core/openeo-udp" + ], + "geometry": null, + "properties": { + "created": "2026-07-13T00:00:00Z", + "updated": "2026-08-11T00:00:00Z", + "type": "service", + "title": "Sentinel-2 based Spectral Angle Mapping (S2-SAM)", + "description": "Computes Spectral Angle Mapping (SAM) from a Sentinel-2 L2A BAP composite to identify the closest matching signature from a Sentinel-2-resampled USGS spectral library for geology and mining related targets.", + "keywords": [ + "Geology", + "Mining", + "Mineral Mapping", + "Sentinel-2", + "Spectral Angle Mapping (SAM)", + "Remote Sensing" + ], + "language": { + "code": "en-US", + "name": "English (United States)" + }, + "languages": [ + { + "code": "en-US", + "name": "English (United States)" + } + ], + "contacts": [ + { + "name": "Christoph Becker", + "position": "Researcher", + "organization": "VITO", + "links": [ + { + "href": "https://www.vito.be/", + "title": "VITO Website", + "rel": "about", + "type": "text/html" + } + ], + "contactInstructions": "Contact via VITO", + "roles": [ + "principal investigator" + ] + } + ], + "themes": [ + { + "concepts": [ + { + "id": "GEOLOGY" + }, + { + "id": "MINING" + }, + { + "id": "REMOTE SENSING" + }, + { + "id": "Sentinel-2 MSI" + } + ], + "scheme": "https://gcmd.earthdata.nasa.gov/kms/concepts/concept_scheme/sciencekeywords" + } + ], + "formats": [ + { + "name": "NetCDF" + } + ], + "license": "Proprietary", + "version": "1.0.0", + "publisher": { + "name": "TERRAVISION" + } + }, + "linkTemplates": [], + "links": [ + { + "rel": "application", + "type": "application/vnd.openeo+json;type=process", + "title": "openEO Process Definition", + "href": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/c0f8d549a866737a3c276f9838fa0a019e7975be/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json" + }, + { + "rel": "code", + "type": "text/html", + "title": "S2-SAM source code", + "href": "https://github.com/ESA-APEx/apex_algorithms" + }, + { + "rel": "service", + "type": "application/json", + "title": "CDSE openEO Federation", + "href": "https://openeofed.dataspace.copernicus.eu" + }, + { + "rel": "platform", + "type": "application/json", + "title": "CDSE openEO Federation processing backend", + "href": "../../../../platform_catalog/cdse_openeo_federation.json" + }, + { + "rel": "provider", + "type": "application/json", + "title": "TERRAVISION", + "href": "../../record.json" + }, + { + "rel": "notebook", + "type": "application/x-ipynb+json", + "title": "S2-SAM example notebook", + "href": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/dc4fe8715fc01539896db75f751278fe6189a59d/algorithm_catalog/terravision/s2_sam/notebook/s2_sam.ipynb" + }, + { + "rel": "webapp", + "type": "text/html", + "title": "OpenEO Web Editor", + "href": "https://editor.openeo.org/?wizard=UDP&wizard%7Eprocess=s2_sam&wizard%7EprocessUrl=https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/c0f8d549a866737a3c276f9838fa0a019e7975be/algorithm_catalog/terravision/s2_sam/openeo_udp/s2_sam.json&server=https%3A%2F%2Fopeneofed.dataspace.copernicus.eu" + }, + { + "rel": "thumbnail", + "type": "image/png", + "title": "Thumbnail image", + "href": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/dc4fe8715fc01539896db75f751278fe6189a59d/algorithm_catalog/terravision/s2_sam/records/thumbnail.png" + } + ] +} \ No newline at end of file diff --git a/algorithm_catalog/terravision/s2_sam/records/thumbnail.png b/algorithm_catalog/terravision/s2_sam/records/thumbnail.png new file mode 100644 index 00000000..2b028e83 Binary files /dev/null and b/algorithm_catalog/terravision/s2_sam/records/thumbnail.png differ