diff --git a/algorithm_catalog/argans/record.json b/algorithm_catalog/argans/record.json new file mode 100644 index 000000000..cd1b846f7 --- /dev/null +++ b/algorithm_catalog/argans/record.json @@ -0,0 +1,68 @@ +{ + "id": "argans", + "type": "Feature", + "conformsTo": [ + "https://www.opengis.net/spec/ogcapi-records-1/1.0/req/record-core" + ], + "properties": { + "created": "2026-04-01T00:00:00Z", + "updated": "2026-04-01T00:00:00Z", + "type": "algorithm_provider", + "title": "Argans Ltd", + "description": "Argans specializes in satellite-based Earth observation, remote sensing, and GIS for monitoring marine, atmospheric, and land environments.", + "keywords": [], + "language": { + "code": "en-US", + "name": "English (United States)" + }, + "languages": [ + { + "code": "en-US", + "name": "English (United States)" + } + ], + "contacts": [ + { + "name": "Argans Ltd", + "emails": [ + { + "value": "enquiries@argans.co.uk" + } + ], + "links": [ + { + "href": "https://argans.co.uk", + "type": "text/html", + "title": "Argans Ltd", + "rel": "website" + } + ] + } + ], + "themes": [], + "acl": { + "admin": ["@argans.co.uk"] + } + }, + "linkTemplates": [], + "links": [ + { + "rel": "website", + "type": "text/html", + "title": "Argans Ltd", + "href": "https://argans.co.uk" + }, + { + "rel": "logo-light", + "type": "image/png", + "title": "Logo", + "href": "https://argans.co.uk/img/logo.png" + }, + { + "rel": "logo-dark", + "type": "image/png", + "title": "Logo", + "href": "https://argans.co.uk/img/logos/argans_white_new.png" + } + ] +} diff --git a/algorithm_catalog/argans/waterlines/benchmark_scenarios/waterlines.json b/algorithm_catalog/argans/waterlines/benchmark_scenarios/waterlines.json new file mode 100644 index 000000000..1e58e59cd --- /dev/null +++ b/algorithm_catalog/argans/waterlines/benchmark_scenarios/waterlines.json @@ -0,0 +1,30 @@ +[ + { + "id": "waterlines", + "type": "openeo", + "backend": "openeofed.dataspace.copernicus.eu", + "process_graph": { + "waterlines1": { + "process_id": "waterlines_v1", + "namespace": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/43d0b17c4ae5bed9db5f9cc5df38efcdcdfb7ecf/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json", + "arguments": { + "temporal_extent": ["2024-06-01", "2024-06-30"], + "spatial_extent": { + "west": -95.13, + "south": 29.078, + "east": -95.12, + "north": 29.082, + "crs": "EPSG:4326" + }, + "iterations": 2, + "max_cloud_coverage": 5, + "ndwi_threshold": 0.01, + "simplify_tolerance": 10 + }, + "result": true + } + }, + "reference_data": { + } + } +] diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/README.md b/algorithm_catalog/argans/waterlines/openeo_udp/README.md new file mode 100644 index 000000000..1ab61cf2f --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/README.md @@ -0,0 +1,48 @@ +# Waterlines openEO UDP +## Purpose +Extract coastline waterlines from Sentinel-2 imagery using NDWI-based water detection, morphological refinement, and UDF-based conversion from water polygons to coast waterlines. + +## Methodology +### Water/Land Classification +Water masks are generated using the **Normalized Difference Water Index (NDWI)**, where pixels are classified as water when **NDWI > threshold**. The default threshold is **0.01**, but it can be adjusted using the `ndwi_threshold` parameter. + +NDWI is computed as the normalized difference between the Sentinel-2 **Green** band (B03) and **Near-Infrared** band (B08), defined as the difference between these bands divided by their sum. + +*This MVP supports only one method (**S2_NDWI**). Originally, multiple methods were selectable via a parameter, but this required openEO `if_()` logic, which converts the result into a `ProcessBuilder` instead of a `DataCube`. +This breaks the `raster_to_vector()` step needed for waterline extraction.* + +### Morphological Processing +For each timestamp, the water/land mask is refined using morphological operations to remove small isolated objects, fill small holes, smooth boundaries and reduce artifacts such as narrow bridges and estuaries. This improves the quality and stability of the resulting waterlines. + +### Waterline Extraction +The cleaned masks are vectorized using the built-in openEO function `raster_to_vector()`. The resulting water polygons are then transformed into waterlines via a UDF, producing time-resolved geometries for each timestep. + +The output is a vector cube of coastline waterlines with the following properties: +- **time**: Acquisition timestamp (Sentinel-2 datetime) +- **type**: Feature type (`waterline_segment`) +- **sea_direction_8**: Sea direction (N, NE, E, SE, S, SW, W, NW) +- **sea_azimuth_deg**: Sea direction in degrees (azimuth, clockwise from north) +- **geometry**: Waterline geometry (LineString or MultiLineString) in EPSG:3857 + +The **sea_azimuth_deg** property is particularly useful for downstream processing, as it can be used to shift the waterline and derive a shoreline (*a waterline normalized for beach slope and tidal conditions*). + +## Authors / Contact +- **Milena Napiorkowska** (openEO UDP) Argans Ltd + mnapiorkowska@argans.co.uk + +- **Martin Jones** (Project Manager) Argans Ltd + mjones@argans.co.uk + +- **Holly Baxter** (Methodology) Argans Ltd + hbaxtar@argans.co.uk + +- **Cameron Mackenzie** (Methodology, openEO UDP) Argans Ltd + cmackenzie@argans.co.uk + +## Acknowledgments +This work was developed as part of an ESA-funded **Fast Track** project. + +## Known Limitations +- Results are most reliable for scenes with low cloud coverage +- NoData areas may introduce artifacts, particularly along boundaries between valid and invalid pixels +- NDWI might be less reliable in turbid waters \ No newline at end of file diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/generate.py b/algorithm_catalog/argans/waterlines/openeo_udp/generate.py new file mode 100644 index 000000000..e5c796c90 --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/generate.py @@ -0,0 +1,179 @@ +import json +from pathlib import Path + +import openeo +from openeo import UDF +from openeo.api.process import Parameter +from openeo.rest.connection import Connection +from openeo.rest.datacube import DataCube +from openeo.rest.udp import build_process_dict + +from s2_index import ( + s2_index_mask, + DEFAULT_S2_COLLECTION, + DEFAULT_MAX_CLOUD_COVER, + WATERLAND_THRESHOLDS, +) + + +def apply_morphology(cube: DataCube, iterations: int) -> DataCube: + udf = UDF.from_file( + Path(__file__).parent / "udf_morph_operations.py", + context={"from_parameter": "context"}, + ) + return cube.apply_dimension( + process=udf, + dimension="t", + context={"iterations": iterations}, + ) + + +def create_waterlines(cube: DataCube, simplify_tolerance: float = 10) -> DataCube: + """ + Extract waterlines from a water/land mask using a UDF. + + The input must remain a DataCube because the workflow relies on + `raster_to_vector()` before applying the vector-based waterline UDF. + """ + cube = cube.raster_to_vector() + + udf = UDF.from_file( + Path(__file__).parent / "udf_waterlines_from_water_land_mask.py", + context={"from_parameter": "context"}, + ) + return cube.apply_dimension( + process=udf, + dimension="geometry", + context={"simplify_tolerance": simplify_tolerance}, + ) + + +def build_water_land_mask_cube( + con: Connection, + bbox, + time_range, + max_cloud_coverage, + iterations, + ndwi_threshold, +) -> DataCube: + """ + Build a water/land mask using Sentinel-2 NDWI only. + + MVP rationale: + Multiple selectable methods were intentionally removed from this UDP. + Selecting between whole DataCubes through nested openEO `if_()` expressions + turns the selected result into a ProcessBuilder rather than a DataCube. + That breaks the next step of the workflow, because `raster_to_vector()` is + used as a DataCube method when preparing input for the vector-based + waterline UDF. + + S2_NDWI was chosen for the MVP because it is a standard water-detection + index, fits the existing index-mask pipeline, and allows validation of the + complete end-to-end workflow without introducing graph-selection complexity. + + Future extensions can add the other methods either as separate UDPs or by + refactoring the waterline UDF to work directly on raster input. + """ + _, cube = s2_index_mask( + con=con, + collection_id=DEFAULT_S2_COLLECTION, + bbox=bbox, + time_range=time_range, + index_name="S2_NDWI", + threshold=ndwi_threshold, + mode="gt", + max_cloud_coverage=max_cloud_coverage, + ) + return apply_morphology(cube, iterations) + + +def generate() -> dict: + """ + Create the MVP UDP for extracting waterlines from Sentinel-2 imagery. + + Workflow: + 1. Load Sentinel-2 data + 2. Create water/land mask using S2_NDWI + 3. Apply morphology + 4. Vectorize the mask + 5. Extract waterlines + + Why only S2_NDWI? + The original multi-method design used a runtime UDP parameter to choose + between several masking methods. In practice, selecting between whole cubes + with openEO graph logic (`if_`) produced a ProcessBuilder instead of a + DataCube. Because the downstream workflow needs `raster_to_vector()`, that + design blocked the current implementation. + + Restricting the MVP to a single method keeps the graph in DataCube form and + allows the existing vector-based waterline UDF to work unchanged. + """ + + conn = openeo.connect(url="openeo.dataspace.copernicus.eu") + + spatial_extent = Parameter.bounding_box( + name="spatial_extent", + description="Bounding box of the area of interest. Defined as west, south, east, north in EPSG:4326.", + ) + + temporal_extent = Parameter.temporal_interval( + name="temporal_extent", + default=["2025-01-01", "2025-12-31"], + description="Date range over which to extract waterlines.", + ) + + max_cloud_coverage = Parameter.number( + name="max_cloud_coverage", + default=DEFAULT_MAX_CLOUD_COVER, + description="Maximum allowed cloud coverage.", + ) + + iterations = Parameter.integer( + name="iterations", + default=2, + description="Number of iterations for morphological operations.", + ) + + ndwi_threshold = Parameter.number( + name="ndwi_threshold", + default=WATERLAND_THRESHOLDS["S2_NDWI"].defaults["threshold"], + description=WATERLAND_THRESHOLDS["S2_NDWI"].description, + ) + + simplify_tolerance = Parameter.number( + name="simplify_tolerance", + default=10, + description="Tolerance used to simplify vectorized water polygons before extracting waterlines.", + ) + + water_land_mask = build_water_land_mask_cube( + con=conn, + bbox=spatial_extent, + time_range=temporal_extent, + max_cloud_coverage=max_cloud_coverage, + iterations=iterations, + ndwi_threshold=ndwi_threshold, + ) + + waterlines_cube = create_waterlines(water_land_mask, simplify_tolerance=simplify_tolerance) + + return build_process_dict( + process_graph=waterlines_cube, + process_id="waterlines_v1", + summary="Waterlines extracted from Sentinel-2 using NDWI.", + description=(Path(__file__).parent / "README.md").read_text(), + parameters=[ + spatial_extent, + temporal_extent, + max_cloud_coverage, + iterations, + ndwi_threshold, + simplify_tolerance, + ], + categories=["sentinel-2", "coastline", "waterlines"], + ) + + +if __name__ == "__main__": + with open(Path(__file__).parent / "waterlines.json", "w") as f: + json.dump(generate(), f, indent=2) diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/s2_index.py b/algorithm_catalog/argans/waterlines/openeo_udp/s2_index.py new file mode 100644 index 000000000..0f2c6d531 --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/s2_index.py @@ -0,0 +1,295 @@ +"""Land/water mask extraction methods using openEO.""" + +from __future__ import annotations +from functools import reduce +from operator import or_ +from enum import Enum +from typing import Mapping, Optional +from dataclasses import dataclass +from openeo import collection_property + +from openeo.rest.connection import Connection +from openeo.rest.datacube import DataCube +from openeo.processes import gt, lt + + +# region defaults + +DEFAULT_S2_COLLECTION = "SENTINEL2_L2A" +DEFAULT_TARGET_EPSG: int = 3857 +DEFAULT_MAX_CLOUD_COVER = 10 + +# endregion + + +# region types + + +class ThresholdMode(str, Enum): + """Comparison direction used when applying a threshold to an index.""" + + GT = "gt" + LT = "lt" + + +@dataclass(frozen=True) +class ThresholdSpec: + """Default threshold(s) and help text for UI/CLI usage.""" + + defaults: Mapping[str, Optional[float]] # keys are method-arg names, e.g. {"threshold": 0.1} + mode: ThresholdMode + description: str + + +# endregion + + +# region registries + + +WATERLAND_THRESHOLDS: dict[str, ThresholdSpec] = { + "S2_NDWI": ThresholdSpec( + defaults={"threshold": 0.01}, + mode=ThresholdMode.GT, + description="NDWI threshold (water if NDWI > threshold).", + ), + "S2_MNDWI": ThresholdSpec( + defaults={"threshold": 0.1}, + mode=ThresholdMode.GT, + description="MNDWI threshold (water if MNDWI > threshold).", + ), + "S2_SCL": ThresholdSpec( + defaults={}, + mode=ThresholdMode.GT, # unused (no threshold), kept for consistency + description="No threshold (SCL class 6 = water).", + ), + "S2_NDVI": ThresholdSpec( + defaults={"threshold": 0.03}, + mode=ThresholdMode.LT, + description="NDVI threshold (water if NDVI < threshold).", + ), + "S2_BNDVI": ThresholdSpec( + defaults={"threshold": 0.03}, + mode=ThresholdMode.LT, + description="BNDVI threshold (water if BNDVI < threshold).", + ), + "S2_GNDVI": ThresholdSpec( + defaults={"threshold": 0.03}, + mode=ThresholdMode.LT, + description="GNDVI threshold (water if GNDVI < threshold).", + ), +} + +# endregion + + +_NORMDIFF_S2: dict[str, tuple[str, str]] = { + # index_name: (band_pos, band_neg) used in (pos - neg) / (pos + neg) + "ndwi": ("B03", "B08"), + "mndwi": ("B03", "B11"), + "ndvi": ("B08", "B04"), + "bndvi": ("B08", "B02"), + "gndvi": ("B08", "B03"), +} + + +def load_collection( + con: Connection, + collection_id: str, + bbox: dict, + time_range: list[str] | None, + bands: list[str] | None = None, + max_cloud_cover: float | None = None, + target_epsg: int | None = None, + resolution: float | tuple[float, float] | None = None, + method: str = "near", + grid_ids: list[str] | None = None, +) -> DataCube: + """Generic openEO collection loader. + + max_cloud_cover is only applied to Sentinel-2 collections. + + Args: + target_epsg: Reproject output to this EPSG code. + resolution: Optional output resolution (single value or (x, y)). + method: Resampling method for reprojection. + """ + + load_kwargs = { + "collection_id": collection_id, + "spatial_extent": bbox, + "temporal_extent": time_range, + "bands": bands, + } + + properties = [] + if grid_ids: + grid_prop = collection_property("grid:code") + grid_filter = reduce(or_, [(grid_prop == gid) for gid in grid_ids]) + properties.append(grid_filter) + + if properties: + load_kwargs["properties"] = properties + + # Apply cloud cover filter only for Sentinel-2 + if "SENTINEL2" in collection_id.upper() and max_cloud_cover is not None: + load_kwargs["max_cloud_cover"] = max_cloud_cover + + cube = con.load_collection(**load_kwargs) + + # Optional reprojection + if target_epsg is not None: + reproj_kwargs = {"projection": target_epsg, "method": method} + if resolution is not None: + reproj_kwargs["resolution"] = resolution + + cube = cube.process("resample_spatial", data=cube, **reproj_kwargs) + + return cube + + +def s2_clear_mask_from_scl(cube: DataCube) -> DataCube: + """Create a boolean clear-pixel mask from the Sentinel-2 SCL band.""" + scl = cube.band("SCL") + # SCL codes: 3 - cloud shadows, 8 - cloud medium prob, 9 - cloud high prob + return (scl == 3 - 1000) | (scl == 8 - 1000) | (scl == 9 - 1000) + + +def load_s2( + con: Connection, + collection_id: str, + bbox: dict, + time_range: list[str], + bands: list[str], + max_cloud_coverage: float | None = DEFAULT_MAX_CLOUD_COVER, + target_epsg: int | None = DEFAULT_TARGET_EPSG, + grid_ids: list[str] | None = None, +) -> tuple[DataCube, DataCube]: + """Load Sentinel-2 with SCL-based cloud masking before optional temporal reduction.""" + bands_with_scl = list(dict.fromkeys(bands + ["SCL"])) # keep order, unique + + # Load S2 bands + SCL + cube = load_collection( + con, + collection_id, + bbox, + time_range, + bands_with_scl, + max_cloud_cover=max_cloud_coverage, + target_epsg=target_epsg, + grid_ids=grid_ids, + ) + + # Cloud masking + clear = s2_clear_mask_from_scl(cube) + cube = cube.mask(clear) + + # Remove SCL band + cube = cube.filter_bands(bands) + + return cube, clear + + +# region private (processing helpers) + + +def _bin(cube: DataCube) -> DataCube: + """Convert a boolean condition cube to a 0/1 cube using an openEO `if` process.""" + return cube.apply(lambda x: x.process("if", arguments={"value": x, "accept": 1, "reject": 0})) + + +# endregion + + +# region public +def s2_scl( + con: Connection, + collection_id: str, + bbox: dict, + time_range: list[str], + max_cloud_coverage: float, +) -> tuple[DataCube, DataCube]: + """Load Sentinel-2 SCL and return a boolean mask selecting SCL class 6. + + Args: + con: openEO connection object. + collection_id: Sentinel-2 collection identifier to load. + bbox: Spatial extent (bounding box) to load. + time_range: Temporal extent as `[start, end]`. + + Returns: + Sentinel-2 cube and + Boolean cube where pixels equal to SCL class 6 are True. + """ + s2_cube, _ = load_s2(con, collection_id, bbox, time_range, bands=["SCL"], max_cloud_coverage=max_cloud_coverage) + scl_water_mask = s2_cube.band("SCL") == 6 + return s2_cube, scl_water_mask + + +def s2_index_mask( + con: Connection, + collection_id: str, + bbox: dict, + time_range: list[str], + index_name: str, + threshold: float | None, + mode: str = "gt", + max_cloud_coverage: float | None = None, + grid_ids: list[str] | None = None, +) -> tuple[DataCube, DataCube]: + """Compute a supported Sentinel-2 norm-diff index and return a binary mask using a threshold. + + Args: + con: openEO connection object. + collection_id: Sentinel-2 collection identifier to load. + bbox: Spatial extent (bounding box) to load. + time_range: Temporal extent as `[start, end]`. + index_name: Name of the index to compute. Supported values: `ndwi`, `mndwi`, + `ndvi`, `bndvi`, `gndvi`. + threshold: Threshold applied to the index. + mode: Threshold mode: + - `"gt"`: pixels where index > threshold become 1 + - `"lt"`: pixels where index < threshold become 1 + max_cloud_coverage: Max allowed cloud coverage. + grid_ids: Filter collection output to these grid IDs. + + Returns: + Sentinel-2 cube and + 0/1 cube representing the index threshold mask. The land pixels + are labelled with 0 and the water pixels are labelled with 1. + + Raises: + ValueError: If `index_name` is not supported or `mode` is not `gt`/`lt`. + """ + key = index_name.lower().split("_")[1] + if key not in _NORMDIFF_S2: + raise ValueError(f"Unsupported index_name={index_name!r}. Supported: {sorted(_NORMDIFF_S2)}") + + band_pos, band_neg = _NORMDIFF_S2[key] + s2_cube, clear = load_s2( + con, + collection_id, + bbox, + time_range, + [band_pos, band_neg], + max_cloud_coverage, + grid_ids=grid_ids, + ) + + pos = s2_cube.band(band_pos) + neg = s2_cube.band(band_neg) + + idx = (pos - neg) / (pos + neg) + idx = idx.mask(clear) + + if mode == "gt": + mask = idx.apply(lambda x: x > threshold) + elif mode == "lt": + mask = idx.apply(lambda x: x < threshold) + else: + raise ValueError(f"Unsupported mode: {mode}") + + return s2_cube, _bin(mask) + + +# endregion diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/udf_morph_operations.py b/algorithm_catalog/argans/waterlines/openeo_udp/udf_morph_operations.py new file mode 100644 index 000000000..ad73574f4 --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/udf_morph_operations.py @@ -0,0 +1,82 @@ +# /// script +# dependencies = [ +# "scikit-image", +# "scipy", +# ] +# /// + +from openeo.udf import XarrayDataCube, inspect +import numpy as np +import xarray as xr +from scipy.ndimage import binary_fill_holes, binary_opening + +DEFAULT_WATER_VALUE = 1 + + +def _build_coastal_water_mask( + arr: np.ndarray, + water_value: int = DEFAULT_WATER_VALUE, + nodata: float | None = 999, + iterations: int = 1, +) -> np.ndarray: + """ + Build coastal-water-only mask from land/water mask. + Inland water is filled. + + Returns: + 0/1 mask where 1 is coastal water and 0 is land + """ + water = arr == water_value + if nodata is not None: + water = water & (arr != nodata) + + # Remove small estuaries + water = binary_opening(water, iterations=iterations) + + # Remove bridges + land = ~water + land = binary_opening(land, iterations=iterations) + + land_filled = binary_fill_holes(land) + water_filled = ~land_filled + return water_filled.astype(np.uint8) + + +def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube: + """Apply morphological algorithms on DataCube""" + + cube_array: xr.DataArray = cube.get_array() + inspect(data=[cube_array.shape], message="Input UDF cube_array shape") + + cube_array = cube_array.astype(np.uint8) + + cube_array_3d = cube_array.squeeze(dim="bands") + + modified = xr.apply_ufunc( + _build_coastal_water_mask, + cube_array_3d, + input_core_dims=[["y", "x"]], + output_core_dims=[["y", "x"]], + vectorize=True, + dask="parallelized", + output_dtypes=[np.uint8], + kwargs={ + "water_value": DEFAULT_WATER_VALUE, + "nodata": 999, + "iterations": context["iterations"], + }, + ) + + modified_da = xr.DataArray( + modified, + coords={ + "t": cube_array.coords["t"], + "y": cube_array.coords["y"], + "x": cube_array.coords["x"], + }, + dims=["t", "y", "x"], + ) + modified_da = modified_da.expand_dims(dim={"bands": cube_array.coords["bands"]}) + modified_da = modified_da.transpose("t", "bands", "y", "x") + + return XarrayDataCube(modified_da) diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/udf_waterlines_from_water_land_mask.py b/algorithm_catalog/argans/waterlines/openeo_udp/udf_waterlines_from_water_land_mask.py new file mode 100644 index 000000000..fd7241dfb --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/udf_waterlines_from_water_land_mask.py @@ -0,0 +1,377 @@ +from typing import Iterable, Union, Any +import numpy as np +import geopandas as gpd +from shapely.geometry import ( + box, + LineString, + MultiLineString, + Polygon, + Point, + GeometryCollection, +) +from openeo.udf import inspect +from openeo.udf.feature_collection import FeatureCollection +from openeo.udf.udf_data import UdfData + +GeometryLike = Union[LineString, MultiLineString, GeometryCollection] + +DEFAULT_OUT_LAYER = "waterline" +DEFAULT_TIME_DIM = "time" +DEFAULT_VAR_NAME = "var" +DEFAULT_SEA_DIRECTION_8_COLUMN = "sea_direction_8" +DEFAULT_SEA_AZIMUTH_DEG_COLUMN = "sea_azimuth_deg" +DEFAULT_MIN_DANGLING_LENGTH = 10000 +DEFAULT_MIN_HOLE_AREA = 1000000 + + +def _iter_lines(geom: GeometryLike) -> Iterable[LineString]: + """Recursively yield all LineString objects contained in a geometry.""" + if geom.is_empty: + return + + if isinstance(geom, LineString): + yield geom + elif isinstance(geom, (MultiLineString, GeometryCollection)): + for subgeom in geom.geoms: + yield from _iter_lines(subgeom) + + +def split_into_segments(geom: GeometryLike) -> list[LineString]: + """ + Split a geometry into 2-point LineStrings representing individual segments + between consecutive vertices. + + Args: + geom: Input geometry. + + Returns: List of non-zero-length segments. + """ + segments: list[LineString] = [] + + for line in _iter_lines(geom): + coords = list(line.coords) + for start, end in zip(coords[:-1], coords[1:]): + if start != end: # Avoid zero-length segments + segments.append(LineString([start, end])) + + return segments + + +def _remove_small_interiors(geom: Polygon, min_hole_area: float = DEFAULT_MIN_HOLE_AREA) -> Polygon: + """Remove small interior rings from polygon.""" + if geom.is_empty: + return geom + + kept_holes = [] + for ring in geom.interiors: + if Polygon(ring).area >= min_hole_area: + kept_holes.append(ring) + + return Polygon(geom.exterior, holes=kept_holes) + + +def _remove_extent_intersections(waterline: LineString, bounds, buffer: float = 600) -> list[LineString]: + """Return 2-point segments that do NOT intersect the raster extent boundary.""" + extent_edge = box(*bounds).boundary + edges = split_into_segments(waterline) + extent_edge_buffered = extent_edge.buffer(buffer) + return [e for e in edges if not e.within(extent_edge_buffered)] + + +def _remove_short_dangling_segments( + segments: list[LineString], + min_dangling_length: float = 0.0, +) -> list[LineString]: + """Remove short isolated segments (both endpoints occur only once).""" + if not segments or min_dangling_length <= 0: + return segments + + endpoint_counts: dict[Any, int] = {} + for seg in segments: + a, b = seg.coords[0], seg.coords[-1] + endpoint_counts[a] = endpoint_counts.get(a, 0) + 1 + endpoint_counts[b] = endpoint_counts.get(b, 0) + 1 + + kept: list[LineString] = [] + for seg in segments: + if seg.length >= min_dangling_length: + kept.append(seg) + continue + + a, b = seg.coords[0], seg.coords[-1] + if endpoint_counts.get(a, 0) > 1 or endpoint_counts.get(b, 0) > 1: + kept.append(seg) + + return kept + + +def _clean_waterline_segments( + waterline: LineString, + bounds, + min_dangling_length: float = DEFAULT_MIN_DANGLING_LENGTH, +) -> list[LineString]: + """ + Clean waterline and return as a *list of 2-point segments* (one per edge). + """ + segments = _remove_extent_intersections(waterline, bounds) + if not segments: + return [] + + segments = _remove_short_dangling_segments( + segments, + min_dangling_length=min_dangling_length, + ) + if not segments: + return [] + + return segments + + +def _get_sea_direction_for_segment(water_poly: Polygon, seg: LineString) -> tuple[str, float | None]: + """ + Determine where the sea (water polygon side) lies relative to a segment. + + Returns: + sea_dir: General sea direction (N, S, NS etc) and detailed sea dir in degrees. + """ + if seg.is_empty or seg.length == 0: + return "unknown", None + + # Get first and last coordinates of the segment + a = np.asarray(seg.coords[0], dtype=float) + b = np.asarray(seg.coords[-1], dtype=float) + + # Compute direction vector and length + v = b - a + norm = np.linalg.norm(v) + if norm == 0: + return "unknown", None + + # Unit tangent and left normal vector + t = v / norm + n_left = np.array([-t[1], t[0]], dtype=float) + + # Segment midpoint + mid = (a + b) / 2.0 + + # How far to step away from the segment (1 perc of segment len.) + eps = max(0.5, min(5.0, float(seg.length) * 0.01)) + + # Probe points + left_pt = mid + eps * n_left + right_pt = mid - eps * n_left + + left_in = water_poly.contains(Point(left_pt)) + right_in = water_poly.contains(Point(right_pt)) + + if left_in and not right_in: + sea_vec = n_left + elif right_in and not left_in: + sea_vec = -n_left + else: + return "unknown", None + + # Map-based 8-way direction from sea_vec (x=east, y=north) + x, y = float(sea_vec[0]), float(sea_vec[1]) + + # Compute angle + angle = np.degrees(np.arctan2(y, x)) + angle = (angle + 360.0) % 360.0 + + # 8-sector compass, centered on E=0°, NE=45°, N=90°, ... + dirs = ["E", "NE", "N", "NW", "W", "SW", "S", "SE"] + idx = int(((angle + 22.5) % 360) // 45) + + return dirs[idx], float(angle) + + +def _segments_for_water_mask( + gdf_water_one_timestamp, + bounds, + simplify_tolerance: float | None = None, +) -> tuple[list[LineString], Polygon] | None: + """Converts water land mask for single timestamp to cleaned waterline segments.""" + + # Remove small interiors + gdf_water_one_timestamp["geometry"] = gdf_water_one_timestamp["geometry"].apply(_remove_small_interiors) + + # Remove small polygons + gdf_water_one_timestamp = gdf_water_one_timestamp[gdf_water_one_timestamp["geometry"].area > DEFAULT_MIN_HOLE_AREA] + + # Merge all polygons + water_poly = gdf_water_one_timestamp.union_all() + + if simplify_tolerance is not None: + water_poly = water_poly.simplify(simplify_tolerance, preserve_topology=True) + + boundary = water_poly.boundary + + # boundary can be MultiLineString/LineString -> convert to a single multilinestring-ish + if isinstance(boundary, LineString): + cleaned_segments = _clean_waterline_segments(boundary, bounds=bounds) + elif isinstance(boundary, MultiLineString): + cleaned_segments = [] + for part in boundary.geoms: + cleaned_segments.extend(_clean_waterline_segments(part, bounds=bounds)) + else: + return None + + return cleaned_segments, water_poly + + +def waterline_from_vectorized_water_raster( + gdf: gpd.GeoDataFrame, + simplify_tolerance: float | None = None, +) -> gpd.GeoDataFrame: + """ + Generate waterline segments for each time step from a vectorized land/water mask. + The input gdf is the output of openEO raster_to_vector() process. + Args: + gdf: GeoDataFrame containing polygon geometries and one non-geometry column + per timestamp. For each timestamp column, values indicate whether a + polygon belongs to water or land at that time: + - null: polygon not present for that timestamp + - 0: land polygon + - non-zero: water polygon + simplify_tolerance: Optional tolerance for geometry simplification. If + provided, merged water geometries are simplified before extracting + boundary segments. + + Returns: + A GeoDataFrame with columns: + - time: Time step associated with each geometry. + - type: Feature classification. + - sea_direction_8: Direction toward the sea expressed as one of + eight cardinal/inter-cardinal directions (N, NE, E, SE, S, SW, W, NW). + - sea_azimuth_deg: Direction toward the sea in degrees (azimuth, + typically measured clockwise from north). + - geometry: Waterline geometry (LineString or MultiLineString). + Raises: + ValueError if no waterlines segments extracted. + """ + + records: list[dict[str, Any]] = [] + + # Identify time columns + time_stamps = gdf.loc[:, gdf.columns != "geometry"].columns.to_list() + inspect(data=[time_stamps], message="Input time stamps.") + + if not time_stamps: + inspect(data=["No time columns found"], message="Empty input") + return gpd.GeoDataFrame(geometry=[], crs=gdf.crs) + + bounds = gdf.total_bounds + + for time_stamp in time_stamps: + inspect(data=[time_stamp], message="Processing timestamp") + + # Drop NaNs + one_time_stamp_gdf = gdf[[time_stamp, "geometry"]].dropna(subset=[time_stamp]) + + if one_time_stamp_gdf.empty: + inspect(data=[time_stamp], message="All values NaN, skipping") + continue + + # Keep only water polygons + water_gdf = one_time_stamp_gdf[one_time_stamp_gdf[time_stamp] != 0] + + if water_gdf.empty: + inspect(data=[time_stamp], message="No water polygons, skipping") + continue + + # Remove invalid geometries early + water_gdf = water_gdf[~water_gdf.geometry.is_empty & water_gdf.geometry.notna()] + + if water_gdf.empty: + inspect(data=[time_stamp], message="No valid geometries after cleaning") + continue + + # Process water mask → segments + res = _segments_for_water_mask( + water_gdf.copy(), + bounds=bounds, + simplify_tolerance=simplify_tolerance, + ) + + if res is None: + inspect(data=[time_stamp], message="No boundary extracted (None), skipping") + continue + + segments, water_poly = res + + if not segments: + inspect(data=[time_stamp], message="No segments found, skipping") + continue + + inspect(data=[len(segments)], message="Segments found") + + for seg in segments: + if seg.is_empty or seg.length == 0: + continue + + sea_direction = _get_sea_direction_for_segment(water_poly, seg) + + records.append( + { + "time": time_stamp, + "type": "waterline_segment", + DEFAULT_SEA_DIRECTION_8_COLUMN: sea_direction[0], + DEFAULT_SEA_AZIMUTH_DEG_COLUMN: sea_direction[1], + "geometry": seg, + } + ) + + # 🔑 CRITICAL CHANGE: never crash on empty result + if len(records) == 0: + inspect( + data=["No waterline segments found for any timestamp"], + message="Returning empty GeoDataFrame", + ) + + return gpd.GeoDataFrame( + columns=[ + "time", + "type", + DEFAULT_SEA_DIRECTION_8_COLUMN, + DEFAULT_SEA_AZIMUTH_DEG_COLUMN, + "geometry", + ], + geometry="geometry", + crs=gdf.crs, + ) + + inspect(data=[len(records)], message="Total segments created") + + result_gdf = gpd.GeoDataFrame(records, geometry="geometry", crs=gdf.crs) + result_gdf = result_gdf[ + ~result_gdf.geometry.isna() & ~result_gdf.geometry.is_empty + ].reset_index(drop=True) + + return result_gdf + +def apply_udf_data(udf_data: UdfData) -> UdfData: + + feature_collection = udf_data.get_feature_collection_list()[0] + gdf = feature_collection.data + inspect(data=[gdf], message="Input gdf data inspection") + + user_context = udf_data.user_context or {} + simplify_tolerance = user_context.get("simplify_tolerance", 10) + + inspect( + data=[simplify_tolerance], + message="Simplify tolerance resolved" + ) + + gdf = waterline_from_vectorized_water_raster( + gdf=gdf, + simplify_tolerance=simplify_tolerance, + ) + + inspect(data=[gdf], message="Output gdf data inspection") + + udf_data.set_feature_collection_list([ + FeatureCollection(id="_", data=gdf) + ]) + + return udf_data diff --git a/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json b/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json new file mode 100644 index 000000000..0e075d3dd --- /dev/null +++ b/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json @@ -0,0 +1,468 @@ +{ + "process_graph": { + "loadcollection1": { + "process_id": "load_collection", + "arguments": { + "bands": [ + "B03", + "B08", + "SCL" + ], + "id": "SENTINEL2_L2A", + "properties": { + "eo:cloud_cover": { + "process_graph": { + "lte1": { + "process_id": "lte", + "arguments": { + "x": { + "from_parameter": "value" + }, + "y": { + "from_parameter": "max_cloud_coverage" + } + }, + "result": true + } + } + } + }, + "spatial_extent": { + "from_parameter": "spatial_extent" + }, + "temporal_extent": { + "from_parameter": "temporal_extent" + } + } + }, + "resamplespatial1": { + "process_id": "resample_spatial", + "arguments": { + "data": { + "from_node": "loadcollection1" + }, + "method": "near", + "projection": 3857 + } + }, + "reducedimension1": { + "process_id": "reduce_dimension", + "arguments": { + "data": { + "from_node": "resamplespatial1" + }, + "dimension": "bands", + "reducer": { + "process_graph": { + "arrayelement1": { + "process_id": "array_element", + "arguments": { + "data": { + "from_parameter": "data" + }, + "index": 2 + } + }, + "eq1": { + "process_id": "eq", + "arguments": { + "x": { + "from_node": "arrayelement1" + }, + "y": -997 + } + }, + "eq2": { + "process_id": "eq", + "arguments": { + "x": { + "from_node": "arrayelement1" + }, + "y": -992 + } + }, + "or1": { + "process_id": "or", + "arguments": { + "x": { + "from_node": "eq1" + }, + "y": { + "from_node": "eq2" + } + } + }, + "eq3": { + "process_id": "eq", + "arguments": { + "x": { + "from_node": "arrayelement1" + }, + "y": -991 + } + }, + "or2": { + "process_id": "or", + "arguments": { + "x": { + "from_node": "or1" + }, + "y": { + "from_node": "eq3" + } + }, + "result": true + } + } + } + } + }, + "mask1": { + "process_id": "mask", + "arguments": { + "data": { + "from_node": "resamplespatial1" + }, + "mask": { + "from_node": "reducedimension1" + } + } + }, + "filterbands1": { + "process_id": "filter_bands", + "arguments": { + "bands": [ + "B03", + "B08" + ], + "data": { + "from_node": "mask1" + } + } + }, + "reducedimension2": { + "process_id": "reduce_dimension", + "arguments": { + "data": { + "from_node": "filterbands1" + }, + "dimension": "bands", + "reducer": { + "process_graph": { + "arrayelement2": { + "process_id": "array_element", + "arguments": { + "data": { + "from_parameter": "data" + }, + "index": 0 + } + }, + "arrayelement3": { + "process_id": "array_element", + "arguments": { + "data": { + "from_parameter": "data" + }, + "index": 1 + } + }, + "subtract1": { + "process_id": "subtract", + "arguments": { + "x": { + "from_node": "arrayelement2" + }, + "y": { + "from_node": "arrayelement3" + } + } + }, + "add1": { + "process_id": "add", + "arguments": { + "x": { + "from_node": "arrayelement2" + }, + "y": { + "from_node": "arrayelement3" + } + } + }, + "divide1": { + "process_id": "divide", + "arguments": { + "x": { + "from_node": "subtract1" + }, + "y": { + "from_node": "add1" + } + }, + "result": true + } + } + } + } + }, + "mask2": { + "process_id": "mask", + "arguments": { + "data": { + "from_node": "reducedimension2" + }, + "mask": { + "from_node": "reducedimension1" + } + } + }, + "apply1": { + "process_id": "apply", + "arguments": { + "data": { + "from_node": "mask2" + }, + "process": { + "process_graph": { + "gt1": { + "process_id": "gt", + "arguments": { + "x": { + "from_parameter": "x" + }, + "y": { + "from_parameter": "ndwi_threshold" + } + }, + "result": true + } + } + } + } + }, + "apply2": { + "process_id": "apply", + "arguments": { + "data": { + "from_node": "apply1" + }, + "process": { + "process_graph": { + "if1": { + "process_id": "if", + "arguments": { + "accept": 1, + "reject": 0, + "value": { + "from_parameter": "x" + } + }, + "result": true + } + } + } + } + }, + "applydimension1": { + "process_id": "apply_dimension", + "arguments": { + "context": { + "iterations": { + "from_parameter": "iterations" + } + }, + "data": { + "from_node": "apply2" + }, + "dimension": "t", + "process": { + "process_graph": { + "runudf1": { + "process_id": "run_udf", + "arguments": { + "context": { + "from_parameter": "context" + }, + "data": { + "from_parameter": "data" + }, + "runtime": "Python", + "udf": "# /// script\n# dependencies = [\n# \"scikit-image\",\n# \"scipy\",\n# ]\n# ///\n\nfrom openeo.udf import XarrayDataCube, inspect\nimport numpy as np\nimport xarray as xr\nfrom scipy.ndimage import binary_fill_holes, binary_opening\n\nDEFAULT_WATER_VALUE = 1\n\n\ndef _build_coastal_water_mask(\n arr: np.ndarray,\n water_value: int = DEFAULT_WATER_VALUE,\n nodata: float | None = 999,\n iterations: int = 1,\n) -> np.ndarray:\n \"\"\"\n Build coastal-water-only mask from land/water mask.\n Inland water is filled.\n\n Returns:\n 0/1 mask where 1 is coastal water and 0 is land\n \"\"\"\n water = arr == water_value\n if nodata is not None:\n water = water & (arr != nodata)\n\n # Remove small estuaries\n water = binary_opening(water, iterations=iterations)\n\n # Remove bridges\n land = ~water\n land = binary_opening(land, iterations=iterations)\n\n land_filled = binary_fill_holes(land)\n water_filled = ~land_filled\n return water_filled.astype(np.uint8)\n\n\ndef apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:\n \"\"\"Apply morphological algorithms on DataCube\"\"\"\n\n cube_array: xr.DataArray = cube.get_array()\n inspect(data=[cube_array.shape], message=\"Input UDF cube_array shape\")\n\n cube_array = cube_array.astype(np.uint8)\n\n cube_array_3d = cube_array.squeeze(dim=\"bands\")\n\n modified = xr.apply_ufunc(\n _build_coastal_water_mask,\n cube_array_3d,\n input_core_dims=[[\"y\", \"x\"]],\n output_core_dims=[[\"y\", \"x\"]],\n vectorize=True,\n dask=\"parallelized\",\n output_dtypes=[np.uint8],\n kwargs={\n \"water_value\": DEFAULT_WATER_VALUE,\n \"nodata\": 999,\n \"iterations\": context[\"iterations\"],\n },\n )\n\n modified_da = xr.DataArray(\n modified,\n coords={\n \"t\": cube_array.coords[\"t\"],\n \"y\": cube_array.coords[\"y\"],\n \"x\": cube_array.coords[\"x\"],\n },\n dims=[\"t\", \"y\", \"x\"],\n )\n modified_da = modified_da.expand_dims(dim={\"bands\": cube_array.coords[\"bands\"]})\n modified_da = modified_da.transpose(\"t\", \"bands\", \"y\", \"x\")\n\n return XarrayDataCube(modified_da)\n" + }, + "result": true + } + } + } + } + }, + "rastertovector1": { + "process_id": "raster_to_vector", + "arguments": { + "data": { + "from_node": "applydimension1" + } + } + }, + "applydimension2": { + "process_id": "apply_dimension", + "arguments": { + "context": { + "simplify_tolerance": { + "from_parameter": "simplify_tolerance" + } + }, + "data": { + "from_node": "rastertovector1" + }, + "dimension": "geometry", + "process": { + "process_graph": { + "runudf2": { + "process_id": "run_udf", + "arguments": { + "context": { + "from_parameter": "context" + }, + "data": { + "from_parameter": "data" + }, + "runtime": "Python", + "udf": "from typing import Iterable, Union, Any\nimport numpy as np\nimport geopandas as gpd\nfrom shapely.geometry import (\n box,\n LineString,\n MultiLineString,\n Polygon,\n Point,\n GeometryCollection,\n)\nfrom openeo.udf import inspect\nfrom openeo.udf.feature_collection import FeatureCollection\nfrom openeo.udf.udf_data import UdfData\n\nGeometryLike = Union[LineString, MultiLineString, GeometryCollection]\n\nDEFAULT_OUT_LAYER = \"waterline\"\nDEFAULT_TIME_DIM = \"time\"\nDEFAULT_VAR_NAME = \"var\"\nDEFAULT_SEA_DIRECTION_8_COLUMN = \"sea_direction_8\"\nDEFAULT_SEA_AZIMUTH_DEG_COLUMN = \"sea_azimuth_deg\"\nDEFAULT_MIN_DANGLING_LENGTH = 10000\nDEFAULT_MIN_HOLE_AREA = 1000000\n\n\ndef _iter_lines(geom: GeometryLike) -> Iterable[LineString]:\n \"\"\"Recursively yield all LineString objects contained in a geometry.\"\"\"\n if geom.is_empty:\n return\n\n if isinstance(geom, LineString):\n yield geom\n elif isinstance(geom, (MultiLineString, GeometryCollection)):\n for subgeom in geom.geoms:\n yield from _iter_lines(subgeom)\n\n\ndef split_into_segments(geom: GeometryLike) -> list[LineString]:\n \"\"\"\n Split a geometry into 2-point LineStrings representing individual segments\n between consecutive vertices.\n\n Args:\n geom: Input geometry.\n\n Returns: List of non-zero-length segments.\n \"\"\"\n segments: list[LineString] = []\n\n for line in _iter_lines(geom):\n coords = list(line.coords)\n for start, end in zip(coords[:-1], coords[1:]):\n if start != end: # Avoid zero-length segments\n segments.append(LineString([start, end]))\n\n return segments\n\n\ndef _remove_small_interiors(geom: Polygon, min_hole_area: float = DEFAULT_MIN_HOLE_AREA) -> Polygon:\n \"\"\"Remove small interior rings from polygon.\"\"\"\n if geom.is_empty:\n return geom\n\n kept_holes = []\n for ring in geom.interiors:\n if Polygon(ring).area >= min_hole_area:\n kept_holes.append(ring)\n\n return Polygon(geom.exterior, holes=kept_holes)\n\n\ndef _remove_extent_intersections(waterline: LineString, bounds, buffer: float = 600) -> list[LineString]:\n \"\"\"Return 2-point segments that do NOT intersect the raster extent boundary.\"\"\"\n extent_edge = box(*bounds).boundary\n edges = split_into_segments(waterline)\n extent_edge_buffered = extent_edge.buffer(buffer)\n return [e for e in edges if not e.within(extent_edge_buffered)]\n\n\ndef _remove_short_dangling_segments(\n segments: list[LineString],\n min_dangling_length: float = 0.0,\n) -> list[LineString]:\n \"\"\"Remove short isolated segments (both endpoints occur only once).\"\"\"\n if not segments or min_dangling_length <= 0:\n return segments\n\n endpoint_counts: dict[Any, int] = {}\n for seg in segments:\n a, b = seg.coords[0], seg.coords[-1]\n endpoint_counts[a] = endpoint_counts.get(a, 0) + 1\n endpoint_counts[b] = endpoint_counts.get(b, 0) + 1\n\n kept: list[LineString] = []\n for seg in segments:\n if seg.length >= min_dangling_length:\n kept.append(seg)\n continue\n\n a, b = seg.coords[0], seg.coords[-1]\n if endpoint_counts.get(a, 0) > 1 or endpoint_counts.get(b, 0) > 1:\n kept.append(seg)\n\n return kept\n\n\ndef _clean_waterline_segments(\n waterline: LineString,\n bounds,\n min_dangling_length: float = DEFAULT_MIN_DANGLING_LENGTH,\n) -> list[LineString]:\n \"\"\"\n Clean waterline and return as a *list of 2-point segments* (one per edge).\n \"\"\"\n segments = _remove_extent_intersections(waterline, bounds)\n if not segments:\n return []\n\n segments = _remove_short_dangling_segments(\n segments,\n min_dangling_length=min_dangling_length,\n )\n if not segments:\n return []\n\n return segments\n\n\ndef _get_sea_direction_for_segment(water_poly: Polygon, seg: LineString) -> tuple[str, float | None]:\n \"\"\"\n Determine where the sea (water polygon side) lies relative to a segment.\n\n Returns:\n sea_dir: General sea direction (N, S, NS etc) and detailed sea dir in degrees.\n \"\"\"\n if seg.is_empty or seg.length == 0:\n return \"unknown\", None\n\n # Get first and last coordinates of the segment\n a = np.asarray(seg.coords[0], dtype=float)\n b = np.asarray(seg.coords[-1], dtype=float)\n\n # Compute direction vector and length\n v = b - a\n norm = np.linalg.norm(v)\n if norm == 0:\n return \"unknown\", None\n\n # Unit tangent and left normal vector\n t = v / norm\n n_left = np.array([-t[1], t[0]], dtype=float)\n\n # Segment midpoint\n mid = (a + b) / 2.0\n\n # How far to step away from the segment (1 perc of segment len.)\n eps = max(0.5, min(5.0, float(seg.length) * 0.01))\n\n # Probe points\n left_pt = mid + eps * n_left\n right_pt = mid - eps * n_left\n\n left_in = water_poly.contains(Point(left_pt))\n right_in = water_poly.contains(Point(right_pt))\n\n if left_in and not right_in:\n sea_vec = n_left\n elif right_in and not left_in:\n sea_vec = -n_left\n else:\n return \"unknown\", None\n\n # Map-based 8-way direction from sea_vec (x=east, y=north)\n x, y = float(sea_vec[0]), float(sea_vec[1])\n\n # Compute angle\n angle = np.degrees(np.arctan2(y, x))\n angle = (angle + 360.0) % 360.0\n\n # 8-sector compass, centered on E=0\u00b0, NE=45\u00b0, N=90\u00b0, ...\n dirs = [\"E\", \"NE\", \"N\", \"NW\", \"W\", \"SW\", \"S\", \"SE\"]\n idx = int(((angle + 22.5) % 360) // 45)\n\n return dirs[idx], float(angle)\n\n\ndef _segments_for_water_mask(\n gdf_water_one_timestamp,\n bounds,\n simplify_tolerance: float | None = None,\n) -> tuple[list[LineString], Polygon] | None:\n \"\"\"Converts water land mask for single timestamp to cleaned waterline segments.\"\"\"\n\n # Remove small interiors\n gdf_water_one_timestamp[\"geometry\"] = gdf_water_one_timestamp[\"geometry\"].apply(_remove_small_interiors)\n\n # Remove small polygons\n gdf_water_one_timestamp = gdf_water_one_timestamp[gdf_water_one_timestamp[\"geometry\"].area > DEFAULT_MIN_HOLE_AREA]\n\n # Merge all polygons\n water_poly = gdf_water_one_timestamp.union_all()\n\n if simplify_tolerance is not None:\n water_poly = water_poly.simplify(simplify_tolerance, preserve_topology=True)\n\n boundary = water_poly.boundary\n\n # boundary can be MultiLineString/LineString -> convert to a single multilinestring-ish\n if isinstance(boundary, LineString):\n cleaned_segments = _clean_waterline_segments(boundary, bounds=bounds)\n elif isinstance(boundary, MultiLineString):\n cleaned_segments = []\n for part in boundary.geoms:\n cleaned_segments.extend(_clean_waterline_segments(part, bounds=bounds))\n else:\n return None\n\n return cleaned_segments, water_poly\n\n\ndef waterline_from_vectorized_water_raster(\n gdf: gpd.GeoDataFrame,\n simplify_tolerance: float | None = None,\n) -> gpd.GeoDataFrame:\n \"\"\"\n Generate waterline segments for each time step from a vectorized land/water mask.\n The input gdf is the output of openEO raster_to_vector() process.\n Args:\n gdf: GeoDataFrame containing polygon geometries and one non-geometry column\n per timestamp. For each timestamp column, values indicate whether a\n polygon belongs to water or land at that time:\n - null: polygon not present for that timestamp\n - 0: land polygon\n - non-zero: water polygon\n simplify_tolerance: Optional tolerance for geometry simplification. If\n provided, merged water geometries are simplified before extracting\n boundary segments.\n\n Returns:\n A GeoDataFrame with columns:\n - time: Time step associated with each geometry.\n - type: Feature classification.\n - sea_direction_8: Direction toward the sea expressed as one of\n eight cardinal/inter-cardinal directions (N, NE, E, SE, S, SW, W, NW).\n - sea_azimuth_deg: Direction toward the sea in degrees (azimuth,\n typically measured clockwise from north).\n - geometry: Waterline geometry (LineString or MultiLineString).\n Raises:\n ValueError if no waterlines segments extracted.\n \"\"\"\n\n records: list[dict[str, Any]] = []\n\n # Identify time columns\n time_stamps = gdf.loc[:, gdf.columns != \"geometry\"].columns.to_list()\n inspect(data=[time_stamps], message=\"Input time stamps.\")\n\n if not time_stamps:\n inspect(data=[\"No time columns found\"], message=\"Empty input\")\n return gpd.GeoDataFrame(geometry=[], crs=gdf.crs)\n\n bounds = gdf.total_bounds\n\n for time_stamp in time_stamps:\n inspect(data=[time_stamp], message=\"Processing timestamp\")\n\n # Drop NaNs\n one_time_stamp_gdf = gdf[[time_stamp, \"geometry\"]].dropna(subset=[time_stamp])\n\n if one_time_stamp_gdf.empty:\n inspect(data=[time_stamp], message=\"All values NaN, skipping\")\n continue\n\n # Keep only water polygons\n water_gdf = one_time_stamp_gdf[one_time_stamp_gdf[time_stamp] != 0]\n\n if water_gdf.empty:\n inspect(data=[time_stamp], message=\"No water polygons, skipping\")\n continue\n\n # Remove invalid geometries early\n water_gdf = water_gdf[~water_gdf.geometry.is_empty & water_gdf.geometry.notna()]\n\n if water_gdf.empty:\n inspect(data=[time_stamp], message=\"No valid geometries after cleaning\")\n continue\n\n # Process water mask \u2192 segments\n res = _segments_for_water_mask(\n water_gdf.copy(),\n bounds=bounds,\n simplify_tolerance=simplify_tolerance,\n )\n\n if res is None:\n inspect(data=[time_stamp], message=\"No boundary extracted (None), skipping\")\n continue\n\n segments, water_poly = res\n\n if not segments:\n inspect(data=[time_stamp], message=\"No segments found, skipping\")\n continue\n\n inspect(data=[len(segments)], message=\"Segments found\")\n\n for seg in segments:\n if seg.is_empty or seg.length == 0:\n continue\n\n sea_direction = _get_sea_direction_for_segment(water_poly, seg)\n\n records.append(\n {\n \"time\": time_stamp,\n \"type\": \"waterline_segment\",\n DEFAULT_SEA_DIRECTION_8_COLUMN: sea_direction[0],\n DEFAULT_SEA_AZIMUTH_DEG_COLUMN: sea_direction[1],\n \"geometry\": seg,\n }\n )\n\n # \ud83d\udd11 CRITICAL CHANGE: never crash on empty result\n if len(records) == 0:\n inspect(\n data=[\"No waterline segments found for any timestamp\"],\n message=\"Returning empty GeoDataFrame\",\n )\n\n return gpd.GeoDataFrame(\n columns=[\n \"time\",\n \"type\",\n DEFAULT_SEA_DIRECTION_8_COLUMN,\n DEFAULT_SEA_AZIMUTH_DEG_COLUMN,\n \"geometry\",\n ],\n geometry=\"geometry\",\n crs=gdf.crs,\n )\n\n inspect(data=[len(records)], message=\"Total segments created\")\n\n result_gdf = gpd.GeoDataFrame(records, geometry=\"geometry\", crs=gdf.crs)\n result_gdf = result_gdf[\n ~result_gdf.geometry.isna() & ~result_gdf.geometry.is_empty\n ].reset_index(drop=True)\n\n return result_gdf\n\ndef apply_udf_data(udf_data: UdfData) -> UdfData:\n\n feature_collection = udf_data.get_feature_collection_list()[0]\n gdf = feature_collection.data\n inspect(data=[gdf], message=\"Input gdf data inspection\")\n\n user_context = udf_data.user_context or {}\n simplify_tolerance = user_context.get(\"simplify_tolerance\", 10)\n\n inspect(\n data=[simplify_tolerance],\n message=\"Simplify tolerance resolved\"\n )\n\n gdf = waterline_from_vectorized_water_raster(\n gdf=gdf,\n simplify_tolerance=simplify_tolerance,\n )\n\n inspect(data=[gdf], message=\"Output gdf data inspection\")\n\n udf_data.set_feature_collection_list([\n FeatureCollection(id=\"_\", data=gdf)\n ])\n\n return udf_data\n" + }, + "result": true + } + } + } + }, + "result": true + } + }, + "id": "waterlines_v1", + "summary": "Waterlines extracted from Sentinel-2 using NDWI.", + "description": "# Waterlines openEO UDP\n## Purpose\nExtract coastline waterlines from Sentinel-2 imagery using NDWI-based water detection, morphological refinement, and UDF-based conversion from water polygons to coast waterlines.\n\n## Methodology\n### Water/Land Classification\nWater masks are generated using the **Normalized Difference Water Index (NDWI)**, where pixels are classified as water when **NDWI > threshold**. The default threshold is **0.01**, but it can be adjusted using the `ndwi_threshold` parameter.\n\nNDWI is computed as the normalized difference between the Sentinel-2 **Green** band (B03) and **Near-Infrared** band (B08), defined as the difference between these bands divided by their sum.\n\n*This MVP supports only one method (**S2_NDWI**). Originally, multiple methods were selectable via a parameter, but this required openEO `if_()` logic, which converts the result into a `ProcessBuilder` instead of a `DataCube`.\nThis breaks the `raster_to_vector()` step needed for waterline extraction.*\n\n### Morphological Processing\nFor each timestamp, the water/land mask is refined using morphological operations to remove small isolated objects, fill small holes, smooth boundaries and reduce artifacts such as narrow bridges and estuaries. This improves the quality and stability of the resulting waterlines.\n\n### Waterline Extraction\nThe cleaned masks are vectorized using the built-in openEO function `raster_to_vector()`. The resulting water polygons are then transformed into waterlines via a UDF, producing time-resolved geometries for each timestep.\n\nThe output is a vector cube of coastline waterlines with the following properties:\n- **time**: Acquisition timestamp (Sentinel-2 datetime) \n- **type**: Feature type (`waterline_segment`) \n- **sea_direction_8**: Sea direction (N, NE, E, SE, S, SW, W, NW) \n- **sea_azimuth_deg**: Sea direction in degrees (azimuth, clockwise from north) \n- **geometry**: Waterline geometry (LineString or MultiLineString) in EPSG:3857 \n\nThe **sea_azimuth_deg** property is particularly useful for downstream processing, as it can be used to shift the waterline and derive a shoreline (*a waterline normalized for beach slope and tidal conditions*).\n\n## Authors / Contact\n- **Milena Napiorkowska** (openEO UDP) Argans Ltd \n mnapiorkowska@argans.co.uk \n\n- **Martin Jones** (Project Manager) Argans Ltd \n mjones@argans.co.uk \n\n- **Holly Baxter** (Methodology) Argans Ltd \n hbaxtar@argans.co.uk \n\n- **Cameron Mackenzie** (Methodology, openEO UDP) Argans Ltd\n cmackenzie@argans.co.uk \n\n## Acknowledgments\nThis work was developed as part of an ESA-funded **Fast Track** project.\n\n## Known Limitations\n- Results are most reliable for scenes with low cloud coverage \n- NoData areas may introduce artifacts, particularly along boundaries between valid and invalid pixels \n- NDWI might be less reliable in turbid waters", + "categories": [ + "sentinel-2", + "coastline", + "waterlines" + ], + "parameters": [ + { + "name": "spatial_extent", + "description": "Bounding box of the area of interest. Defined as west, south, east, north in EPSG:4326.", + "schema": { + "type": "object", + "subtype": "bounding-box", + "required": [ + "west", + "south", + "east", + "north" + ], + "properties": { + "west": { + "type": "number", + "description": "West (lower left corner, coordinate axis 1)." + }, + "south": { + "type": "number", + "description": "South (lower left corner, coordinate axis 2)." + }, + "east": { + "type": "number", + "description": "East (upper right corner, coordinate axis 1)." + }, + "north": { + "type": "number", + "description": "North (upper right corner, coordinate axis 2)." + }, + "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": [ + { + "type": "integer", + "subtype": "epsg-code", + "title": "EPSG Code", + "minimum": 1000 + }, + { + "type": "string", + "subtype": "wkt2-definition", + "title": "WKT2 definition" + } + ], + "default": 4326 + } + } + } + }, + { + "name": "temporal_extent", + "description": "Date range over which to extract waterlines.", + "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" + } + ] + } + }, + "default": [ + "2025-01-01", + "2025-12-31" + ], + "optional": true + }, + { + "name": "max_cloud_coverage", + "description": "Maximum allowed cloud coverage.", + "schema": { + "type": "number" + }, + "default": 10, + "optional": true + }, + { + "name": "iterations", + "description": "Number of iterations for morphological operations.", + "schema": { + "type": "integer" + }, + "default": 2, + "optional": true + }, + { + "name": "ndwi_threshold", + "description": "NDWI threshold (water if NDWI > threshold).", + "schema": { + "type": "number" + }, + "default": 0.01, + "optional": true + }, + { + "name": "simplify_tolerance", + "description": "Tolerance used to simplify vectorized water polygons before extracting waterlines.", + "schema": { + "type": "number" + }, + "default": 10, + "optional": true + } + ] +} \ No newline at end of file diff --git a/algorithm_catalog/argans/waterlines/records/thumbnail.png b/algorithm_catalog/argans/waterlines/records/thumbnail.png new file mode 100644 index 000000000..4b69c2ee9 Binary files /dev/null and b/algorithm_catalog/argans/waterlines/records/thumbnail.png differ diff --git a/algorithm_catalog/argans/waterlines/records/waterlines.json b/algorithm_catalog/argans/waterlines/records/waterlines.json new file mode 100644 index 000000000..7205ee9fd --- /dev/null +++ b/algorithm_catalog/argans/waterlines/records/waterlines.json @@ -0,0 +1,123 @@ +{ + "id": "waterlines", + "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-04-01T00:00:00Z", + "updated": "2026-04-16T00:00:00Z", + "type": "service", + "title": "Waterlines from Sentinel-2", + "description": "Extracts coastal waterlines from Sentinel-2 time series by generating NDWI-based land-water masks, vectorising them, and converting them into waterlines.", + "keywords": [ + "Normalized Difference Water Index (NDWI)", + "Natural Hazards", + "Sentinel-2", + "Coastal" + ], + "language": { + "code": "en-US", + "name": "English (United States)" + }, + "languages": [ + { + "code": "en-US", + "name": "English (United States)" + } + ], + "contacts": [ + { + "name": "Argans Ltd", + "organization": "Argans Ltd", + "links": [ + { + "href": "https://argans.co.uk", + "title": "Argans Ltd", + "rel": "about", + "type": "text/html" + } + ], + "contactInstructions": "Contact via Argans Ltd", + "roles": [ + "processor" + ] + } + ], + "themes": [ + { + "concepts": [ + { + "id": "COASTAL PROCESSES" + }, + { + "id": "REMOTE SENSING" + }, + { + "id": "Sentinel-2 MSI" + } + ], + "scheme": "https://gcmd.earthdata.nasa.gov/kms/concepts/concept_scheme/sciencekeywords" + } + ], + "formats": [ + { + "name": "GeoJSON" + }, + { + "name": "Parquet" + }, + { + "name": "GPKG" + } + ], + "license": "CC-BY-4.0" + }, + "linkTemplates": [], + "links": [ + { + "rel": "application", + "type": "application/vnd.openeo+json;type=process", + "title": "openEO Process Definition", + "href": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/refs/heads/argans_waterlines/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json" + }, + { + "rel": "code", + "type": "text/html", + "title": "openeo-udp repository", + "href": "" + }, + { + "rel": "webapp", + "type": "text/html", + "title": "OpenEO Web Editor", + "href": "https://editor.openeo.org/?wizard=UDP&wizard%7Eprocess=waterlines_v1&wizard%7EprocessUrl=https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/refs/heads/argans_waterlines/algorithm_catalog/argans/waterlines/openeo_udp/waterlines.json&server=openeofed.dataspace.copernicus.eu" + }, + { + "rel": "service", + "type": "application/json", + "title": "CDSE openEO Federation", + "href": "https://openeofed.dataspace.copernicus.eu" + }, + { + "rel": "platform", + "type": "application/json", + "title": "CDSE openEO federation", + "href": "../../../../platform_catalog/cdse_openeo_federation.json" + }, + { + "rel": "provider", + "type": "application/json", + "title": "Argans Ltd", + "href": "../../record.json" + }, + { + "rel": "thumbnail", + "type": "image/png", + "title": "Thumbnail image", + "href": "https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/refs/heads/argans_waterlines/algorithm_catalog/argans/waterlines/records/thumbnail.png" + } + ] +} diff --git a/schemas/record.json b/schemas/record.json index c077d594d..00feae192 100644 --- a/schemas/record.json +++ b/schemas/record.json @@ -135,7 +135,8 @@ "Land Surface Temperature", "Climate", "Soils", - "Radar" + "Radar", + "Coastal" ] }, "minItems": 1,