opendatacube / odc-stac

Load STAC items into xarray Datasets.
Apache License 2.0
136 stars 19 forks source link

Data gap when loading items across the antimeridian #165

Open jessjaco opened 4 weeks ago

jessjaco commented 4 weeks ago

I am trying to load data for a GeoBox that crosses the antimeridian

from affine import Affine
from odc.geo.cog import write_cog
from odc.geo.geobox import GeoBox
import odc.stac
from pyproj import CRS
from pystac import Item, ItemCollection
import stackstac

item_urls = [
    "https://earth-search.aws.element84.com/v1/collections/cop-dem-glo-30/items/Copernicus_DSM_COG_10_S18_00_W180_00_DEM",
    "https://earth-search.aws.element84.com/v1/collections/cop-dem-glo-30/items/Copernicus_DSM_COG_10_S18_00_E179_00_DEM",
    "https://earth-search.aws.element84.com/v1/collections/cop-dem-glo-30/items/Copernicus_DSM_COG_10_S17_00_W180_00_DEM",
    "https://earth-search.aws.element84.com/v1/collections/cop-dem-glo-30/items/Copernicus_DSM_COG_10_S17_00_E179_00_DEM",
]

items = ItemCollection([Item.from_file(url) for url in item_urls])

geobox = GeoBox(
    (2491, 3200), Affine(30.0, 0.0, 3336000.0, 0.0, -30.0, -1887990.0), CRS("EPSG:3832")
)

odc_ds = odc.stac.load(items, geobox=geobox, chunks=dict(x=2048, y=2048))
write_cog(odc_ds.data, "odc_load.tif")

Unfortunately this produces a gap at the antimeridian: image

I can confirm this is not an issue with the assets because stackstac can load them correctly.

ss_da = (
    stackstac.stack(
        items, epsg=3832, chunksize=dict(x=2048, y=2048), bounds=geobox.boundingbox.bbox
    )
    .groupby("time")
    .first()
    .squeeze()
)

write_cog(ss_da, "ss_load.tif")

image

I suspect it may be because the bounding box of two of the stac items crosses -180 (see e.g. https://earth-search.aws.element84.com/v1/collections/cop-dem-glo-30/items/Copernicus_DSM_COG_10_S18_00_W180_00_DEM) but that's just a guess.
image

odc-stac==0.3.10

alexgleith commented 2 weeks ago

Looks like an issue with the GeoBox.

from affine import Affine
from odc.geo.geobox import GeoBox
from pyproj import CRS

geobox = GeoBox(
    (2491, 3200), Affine(30.0, 0.0, 3336000.0, 0.0, -30.0, -1887990.0), CRS("EPSG:3832")
)

list(geobox.extent.to_crs("epsg:4326").boundingbox)
[-179.169819449018,
 -17.469444408354423,
 179.96779787822723,
 -16.824027928472105]

this is world spanning!