Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 25 additions & 16 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -46,56 +46,65 @@ 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;

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 \
Expand Down
4 changes: 1 addition & 3 deletions ldn/cli_geomad.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
56 changes: 0 additions & 56 deletions ldn/geomad.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__)
Expand Down Expand Up @@ -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):

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@caitlinadams I could split this change to another PR. It is an mostly unrelated refactor. Just moving this shared class to a different file

"""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
4 changes: 3 additions & 1 deletion ldn/lulc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand Down
79 changes: 78 additions & 1 deletion ldn/raster.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
27 changes: 27 additions & 0 deletions ldn/tests/test_raster_antimeridian_bbox.py
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading