diff --git a/Makefile b/Makefile index 351744e2..a6b06c2e 100644 --- a/Makefile +++ b/Makefile @@ -46,41 +46,50 @@ print-tasks-test-dep-staging: ldn print-tasks \ --years="2000" \ --region="pacific" \ - --geomad-version 0-3-0-test \ + --geomad-version 0-3-1-test \ --dataset geomad \ --no-overwrite \ --bucket dep-public-staging; geomad-test-ausp: ldn geomad run \ - --tile-id 031_038 \ + --tile-id 064_020 \ --region pacific \ --year 2000 \ - --version 0-3-0-test \ + --version 0-3-1-test \ --decimated \ --bucket data.ldn.auspatious.com \ --overwrite; + +# To test Antimeridian, use Pacific tiles: +# 065_020 (just before the antimeridian) +# 066_020 (crosses the antimeridian) +# 067_020 (just after the antimeridian) geomad-test-dep-staging: - ldn geomad run \ - --tile-id 031_038 \ - --region pacific \ - --year 2000 \ - --version 0-3-0-test \ - --collection-url-root="https://stac.staging.digitalearthpacific.io/collections" \ - --decimated \ - --bucket dep-public-staging \ - --overwrite; + for col in 064 065 066 067; do \ + for row in 020 021 022; do \ + ldn geomad run \ + --tile-id $${col}_$${row} \ + --region pacific \ + --year 2025 \ + --version 0-3-1-test \ + --collection-url-root="https://stac.staging.digitalearthpacific.io/collections" \ + --decimated \ + --bucket dep-public-staging \ + --no-overwrite; \ + done; \ + done; index-geomad-test-ausp: ldn index-to-stac-geoparquet \ --dataset geomad \ - --geomad-version 0-3-0-test \ + --geomad-version 0-3-1-test \ --no-single-region \ --bucket data.ldn.auspatious.com; index-geomad-test-dep-staging: ldn index-to-stac-geoparquet \ --dataset geomad \ - --geomad-version 0-3-0-test \ + --geomad-version 0-3-1-test \ --single-region \ --product-owner dep \ --bucket dep-public-staging; @@ -88,14 +97,14 @@ index-geomad-test-dep-staging: collection-geomad-test-ausp: ldn collection create-collection \ --dataset geomad \ - --geomad-version 0-3-0-test \ + --geomad-version 0-3-1-test \ --no-single-region \ --bucket data.ldn.auspatious.com \ --no-has-stac-api; collection-geomad-test-dep-staging: ldn collection create-collection \ --dataset geomad \ - --geomad-version 0-3-0-test \ + --geomad-version 0-3-1-test \ --url-root="https://stac.staging.digitalearthpacific.io" \ --single-region \ --product-owner dep \ diff --git a/ldn/cli_geomad.py b/ldn/cli_geomad.py index d134080f..d291428c 100644 --- a/ldn/cli_geomad.py +++ b/ldn/cli_geomad.py @@ -20,10 +20,8 @@ GeoMADProcessor, InsufficientScenesError, ) -from ldn.geomad import ( - AwsStacTask as Task, -) from ldn.grids import get_gridspec +from ldn.raster import AwsStacTask as Task from ldn.raster import build_pipeline_components from ldn.utils import ( GEOMAD_DATASET_ID, diff --git a/ldn/geomad.py b/ldn/geomad.py index 2f38b32b..36ea9f30 100644 --- a/ldn/geomad.py +++ b/ldn/geomad.py @@ -4,17 +4,10 @@ import numpy as np from datacube_compute import geomedian_with_mads -from dep_tools.loaders import StacLoader from dep_tools.processors import Processor -from dep_tools.searchers import Searcher -from dep_tools.stac_utils import StacCreator -from dep_tools.task import AreaTask -from dep_tools.writers import AwsDsCogWriter, AwsStacWriter from odc.algo import mask_cleanup -from odc.geo import GeoBox from xarray import DataArray, Dataset -from ldn.raster import PrefixedS3ItemPath from ldn.utils import LdnError logger = logging.getLogger(__name__) @@ -372,52 +365,3 @@ def process(self, ds: Dataset) -> Dataset: geomad = geomad.compute() return _set_stac_properties(data, geomad) - - -# This is a generic function used be geomad creation and lulc classification tasks. -# TODO: Move it to raster.py -class AwsStacTask(AreaTask): - """Area task with search + STAC creation/writing for AWS workflows.""" - - def __init__( - self, - itempath: PrefixedS3ItemPath, - id: tuple[int, int], - area: GeoBox, - searcher: Searcher, - loader: StacLoader, - processor: Processor, - post_processor: Processor | None = None, - logger: logging.Logger = logger, - **kwargs, - ): - writer = kwargs.pop("writer", AwsDsCogWriter(itempath)) - stac_creator = kwargs.pop("stac_creator", StacCreator(itempath)) - stac_writer = kwargs.pop("stac_writer", AwsStacWriter(itempath)) - - super().__init__(id, area, loader, processor, writer, logger) - self.id = id - self.searcher = searcher - self.post_processor = post_processor - self.stac_creator = stac_creator - self.stac_writer = stac_writer - - def run(self): - items = self.searcher.search(self.area) - logger.info(f"Found {len(items)} items for this tile/year") - input_data = self.loader.load(items, self.area) - logger.info(f"Loaded {len(input_data.time.values)} items for this tile/year") - - processor_kwargs = dict(area=self.area) if self.processor.send_area_to_processor else dict() - output_data = self.processor.process(input_data, **processor_kwargs) - - if self.post_processor is not None: - output_data = self.post_processor.process(output_data) - - paths = self.writer.write(output_data, self.id) - - if self.stac_creator is not None and self.stac_writer is not None: - stac_item = self.stac_creator.process(output_data, self.id) - self.stac_writer.write(stac_item, self.id) - - return paths diff --git a/ldn/lulc.py b/ldn/lulc.py index 258440d5..829d933d 100644 --- a/ldn/lulc.py +++ b/ldn/lulc.py @@ -23,7 +23,6 @@ from typing_extensions import Annotated from ldn.aws import configure_s3_access_profile, s3_client -from ldn.geomad import AwsStacTask as Task from ldn.grids import get_gridspec from ldn.raster import ( GEOMAD_BANDS, @@ -32,6 +31,9 @@ load_dem_terrain, scale_offset_landsat, ) +from ldn.raster import ( + AwsStacTask as Task, +) from ldn.utils import ( GEOMAD_DATASET_ID, GEOMAD_VERSION, diff --git a/ldn/raster.py b/ldn/raster.py index d3be0a27..980b2e22 100644 --- a/ldn/raster.py +++ b/ldn/raster.py @@ -3,11 +3,16 @@ import numpy as np import xarray as xr +from antimeridian import fix_shape from dep_tools.aws import BaseClient, object_exists +from dep_tools.loaders import StacLoader from dep_tools.namers import S3ItemPath +from dep_tools.processors import Processor +from dep_tools.searchers import Searcher from dep_tools.stac_utils import StacCreator +from dep_tools.task import AreaTask from dep_tools.utils import bbox_across_180, join_path_or_url, search_across_180 -from dep_tools.writers import AwsDsCogWriter +from dep_tools.writers import AwsDsCogWriter, AwsStacWriter from geopandas import GeoDataFrame from odc.geo.geobox import GeoBox from odc.stac import load as stac_load @@ -346,3 +351,75 @@ def build_pipeline_components( ) return itempath, stac_creator, writer + + +def _antimeridian_safe_bbox(area: GeoBox, fallback_bbox: list[float]) -> list[float]: + """Return a STAC-correct bbox for a tile, fixing antimeridian-crossing cases. + + rio_stac derives the STAC bbox from the raster's own reprojected corners, + which for a tile crossing the antimeridian collapses to a near-global + (-180, 180) box instead of the STAC-recommended [east_edge, ymin, + west_edge, ymax] form (bbox[0] > bbox[2]). `bbox_across_180` already knows + how to split such a tile into two non-crossing halves; here we just glue + the outer edges of those halves back into one antimeridian-flipped bbox. + """ + halves = bbox_across_180(area) + if not isinstance(halves, tuple): + return fallback_bbox + east_bbox, west_bbox = halves + return [ + east_bbox[0], + min(east_bbox[1], west_bbox[1]), + west_bbox[2], + max(east_bbox[3], west_bbox[3]), + ] + + +# Shared by both GeoMAD creation and LULC classification tasks. +class AwsStacTask(AreaTask): + """Area task with search + STAC creation/writing for AWS workflows.""" + + def __init__( + self, + itempath: PrefixedS3ItemPath, + id: tuple[int, int], + area: GeoBox, + searcher: Searcher, + loader: StacLoader, + processor: Processor, + post_processor: Processor | None = None, + logger: logging.Logger = logger, + **kwargs, + ): + writer = kwargs.pop("writer", AwsDsCogWriter(itempath)) + stac_creator = kwargs.pop("stac_creator", StacCreator(itempath)) + stac_writer = kwargs.pop("stac_writer", AwsStacWriter(itempath)) + + super().__init__(id, area, loader, processor, writer, logger) + self.id = id + self.searcher = searcher + self.post_processor = post_processor + self.stac_creator = stac_creator + self.stac_writer = stac_writer + + def run(self): + items = self.searcher.search(self.area) + logger.info(f"Found {len(items)} items for this tile/year") + input_data = self.loader.load(items, self.area) + logger.info(f"Loaded {len(input_data.time.values)} items for this tile/year") + + processor_kwargs = dict(area=self.area) if self.processor.send_area_to_processor else dict() + output_data = self.processor.process(input_data, **processor_kwargs) + + if self.post_processor is not None: + output_data = self.post_processor.process(output_data) + + paths = self.writer.write(output_data, self.id) + + if self.stac_creator is not None and self.stac_writer is not None: + stac_item = self.stac_creator.process(output_data, self.id) + stac_item.bbox = _antimeridian_safe_bbox(self.area, stac_item.bbox) + stac_item.geometry = fix_shape(stac_item.geometry) + self.stac_writer.write(stac_item, self.id) + + return paths diff --git a/ldn/tests/test_raster_antimeridian_bbox.py b/ldn/tests/test_raster_antimeridian_bbox.py new file mode 100644 index 00000000..a4eb270e --- /dev/null +++ b/ldn/tests/test_raster_antimeridian_bbox.py @@ -0,0 +1,27 @@ +from odc.geo.geobox import GeoBox + +from ldn.raster import _antimeridian_safe_bbox + +# EPSG:3832 (WGS 84 / PDC Mercator) x=3_330_000 is roughly where longitude +# wraps from +180 to -180 at this latitude. +_AM_CROSSING_BBOX_3832 = (3_285_000, -2_015_000, 3_375_000, -1_925_000) +_NON_CROSSING_BBOX_3832 = (0, -2_015_000, 90_000, -1_925_000) + + +def test_antimeridian_safe_bbox_flips_crossing_tile() -> None: + """A tile straddling the antimeridian gets a flipped bbox (bbox[0] > bbox[2]).""" + gb = GeoBox.from_bbox(_AM_CROSSING_BBOX_3832, crs="EPSG:3832", resolution=100) + + fixed = _antimeridian_safe_bbox(gb, fallback_bbox=[-180, -18, 180, -17]) + + assert fixed[0] > fixed[2], "expected STAC antimeridian convention: east edge > west edge" + assert fixed[0] == 179.50965708332626 + assert fixed[2] == -179.68185916096616 + + +def test_antimeridian_safe_bbox_passes_through_non_crossing_tile() -> None: + """A tile that doesn't cross the antimeridian keeps rio_stac's original bbox.""" + gb = GeoBox.from_bbox(_NON_CROSSING_BBOX_3832, crs="EPSG:3832", resolution=100) + fallback = [149.9, -18, 150.1, -17] + + assert _antimeridian_safe_bbox(gb, fallback_bbox=fallback) == fallback diff --git a/notebooks/antimeridian_stac_test.ipynb b/notebooks/antimeridian_stac_test.ipynb new file mode 100644 index 00000000..1955c481 --- /dev/null +++ b/notebooks/antimeridian_stac_test.ipynb @@ -0,0 +1,306 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7fb27b941602401d91542211134fc71a", + "metadata": {}, + "source": [ + "# Antimeridian bbox test: pystac → stac-geoparquet → rustac search\n", + "\n", + "Creates a `pystac.Item` with an antimeridian-crossing bbox, writes it to a local\n", + "stac-geoparquet file, then uses `rustac` to run point-intersection searches.\n", + "\n", + "**Note on the bbox:** a literal `[170, -10, 170, 10]` has `west == east` (a\n", + "zero-width box), which can't be the intended antimeridian-crossing shape. This\n", + "notebook uses the reversed-longitude convention from RFC 7946 §5.2 / the STAC\n", + "API spec instead: **`[170, -10, -170, 10]`** (west=170°, east=-170°), which\n", + "covers the ~20°-wide sliver straddling ±180° between 170°E and 170°W —\n", + "consistent with the Fiji-straddling example used throughout this conversation.\n", + "\n", + "Expected results for the three test points:\n", + "- `[175, -5]` → **True** (170° to 180° side of the crossing)\n", + "- `[-175, 5]` → **True** (-180° to -170° side of the crossing)\n", + "- `[0, 0]` → **False** (nowhere near the antimeridian)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acae54e37e7d407bbb7b55eff062a284", + "metadata": {}, + "outputs": [], + "source": [ + "# Install dependencies (skip if already installed)\n", + "%pip install -q pystac rustac" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "9a63283cbaf04dbcab1f6479b197f3a8", + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime, timezone\n", + "\n", + "import pystac\n", + "import rustac" + ] + }, + { + "cell_type": "markdown", + "id": "8dd0d8092fe74a7c96281538738b07e2", + "metadata": {}, + "source": [ + "## 1. Build the pystac Item\n", + "\n", + "`bbox` uses the reversed convention (`west=170 > east=-170`). `geometry` is a\n", + "`MultiPolygon` split at the seam — one piece from 170° to 180°, one piece from\n", + "-180° to -170° — which is the standard way to represent an antimeridian-crossing\n", + "footprint without a self-intersecting ring. pystac does not compute or validate\n", + "this for you (see earlier discussion): both `bbox` and `geometry` are supplied\n", + "explicitly here." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "72eea5119410473aa328ad9291626812", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Item id: antimeridian-test-item\n", + "Item bbox: [170, -10, -170, 10]\n" + ] + } + ], + "source": [ + "bbox = [170, -10, -170, 10] # west, south, east, north (reversed: west > east)\n", + "\n", + "geometry = {\n", + " \"type\": \"MultiPolygon\",\n", + " \"coordinates\": [\n", + " # West piece: 170°E to 180°\n", + " [[[170, -10], [180, -10], [180, 10], [170, 10], [170, -10]]],\n", + " # East piece: -180° to -170° (i.e. 170°W)\n", + " [[[-180, -10], [-170, -10], [-170, 10], [-180, 10], [-180, -10]]],\n", + " ],\n", + "}\n", + "\n", + "item = pystac.Item(\n", + " id=\"antimeridian-test-item\",\n", + " geometry=geometry,\n", + " bbox=bbox,\n", + " datetime=datetime.now(timezone.utc),\n", + " properties={},\n", + ")\n", + "\n", + "item_dict = item.to_dict()\n", + "print(\"Item id: \", item_dict[\"id\"])\n", + "print(\"Item bbox:\", item_dict[\"bbox\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "12831183", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'type': 'Feature',\n", + " 'stac_version': '1.1.0',\n", + " 'stac_extensions': [],\n", + " 'id': 'antimeridian-test-item',\n", + " 'geometry': {'type': 'MultiPolygon',\n", + " 'coordinates': [[[[170, -10], [180, -10], [180, 10], [170, 10], [170, -10]]],\n", + " [[[-180, -10], [-170, -10], [-170, 10], [-180, 10], [-180, -10]]]]},\n", + " 'bbox': [170, -10, -170, 10],\n", + " 'properties': {'datetime': '2026-09-21T04:25:53.629291Z'},\n", + " 'links': [],\n", + " 'assets': {}}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "item_dict" + ] + }, + { + "cell_type": "markdown", + "id": "8edb47106e1a46a883d545849b8ab81b", + "metadata": {}, + "source": [ + "## 2. Write to a local stac-geoparquet file\n", + "\n", + "`rustac.write` takes a list of STAC item dicts and writes them to a\n", + "stac-geoparquet file. `rustac`'s async functions can be awaited directly at\n", + "the top level of a notebook cell." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "10185d26023b46108eb7d9f57d49d2b3", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Wrote 1 item to items.parquet\n" + ] + } + ], + "source": [ + "parquet_path = \"items.parquet\"\n", + "\n", + "await rustac.write(parquet_path, [item_dict])\n", + "print(f\"Wrote 1 item to {parquet_path}\")" + ] + }, + { + "cell_type": "markdown", + "id": "8763a12b2bbd4a93a75aff182afb95dc", + "metadata": {}, + "source": [ + "### Sanity check: read it back\n", + "\n", + "Confirms the reversed bbox and split `MultiPolygon` geometry survive the\n", + "geoparquet round-trip unchanged." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "7623eae2785240b9bd12b16a66d81610", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bbox from geoparquet: (170.0, -10.0, -170.0, 10.0)\n", + "geometry from geoparquet: {'type': 'MultiPolygon', 'coordinates': [[[[170.0, -10.0], [180.0, -10.0], [180.0, 10.0], [170.0, 10.0], [170.0, -10.0]]], [[[-180.0, -10.0], [-170.0, -10.0], [-170.0, 10.0], [-180.0, 10.0], [-180.0, -10.0]]]]}\n" + ] + } + ], + "source": [ + "read_back = await rustac.read(parquet_path)\n", + "feature = read_back[\"features\"][0]\n", + "\n", + "print(\"bbox from geoparquet: \", feature[\"bbox\"])\n", + "print(\"geometry from geoparquet:\", feature[\"geometry\"])" + ] + }, + { + "cell_type": "markdown", + "id": "7cdc8c89c7104fffa095e18ddfef8986", + "metadata": {}, + "source": [ + "## 3. Search with rustac for each test point\n", + "\n", + "`rustac.search` can query a local stac-geoparquet file directly (no STAC API\n", + "server needed) using `intersects` with a GeoJSON `Point`. Under the hood this\n", + "uses DuckDB's spatial extension, which is downloaded automatically on first\n", + "use — that first call needs internet access." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "b118ea5561624da68c537baed56e602f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "✅ [175, -5] (east hemisphere piece, 170→180) matched=1 expected=1\n", + "✅ [-175, 5] (west hemisphere piece, -180→-170) matched=1 expected=1\n", + "✅ [165, 5] (not intersecting) matched=0 expected=0\n", + "✅ [180, 0] (interesection on AM) matched=1 expected=1\n" + ] + } + ], + "source": [ + "test_points = {\n", + " \"[175, -5] (east hemisphere piece, 170→180)\": ([175, -5], 1),\n", + " \"[-175, 5] (west hemisphere piece, -180→-170)\": ([-175, 5], 1),\n", + " \"[165, 5] (not intersecting)\": ([165, 5], 0),\n", + " \"[180, 0] (interesection on AM)\": ([180, 0], 1),\n", + "}\n", + "\n", + "results = {}\n", + "\n", + "for label, (coords, expected) in test_points.items():\n", + " result = await rustac.search(\n", + " parquet_path,\n", + " intersects={\"type\": \"Point\", \"coordinates\": coords},\n", + " )\n", + " matched = len(result)\n", + " results[label] = matched\n", + " status = \"✅\" if matched == expected else \"❌\"\n", + " print(f\"{status} {label:38s} matched={matched} expected={expected}\")" + ] + }, + { + "cell_type": "markdown", + "id": "938c804e27f84196a10c8828c723f798", + "metadata": {}, + "source": [ + "## 4. Assert the expected pattern: True, True, False" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "504fb2a444614c0babb325280ed9130a", + "metadata": {}, + "outputs": [], + "source": [ + "expected_pattern = [True, True, False]\n", + "actual_pattern = list(results.values())\n", + "\n", + "print(\"Expected:\", expected_pattern)\n", + "print(\"Actual: \", actual_pattern)\n", + "\n", + "assert actual_pattern == expected_pattern, (\n", + " \"rustac's search over the antimeridian-crossing item did not match the \"\n", + " \"expected True/True/False pattern — the reversed bbox/geometry may not be \"\n", + " \"handled correctly.\"\n", + ")\n", + "print(\"\\nAll assertions passed — rustac correctly resolves the antimeridian-crossing item.\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ldn (3.14.6)", + "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 +}