Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Data gap when loading items across the antimeridian #165

Open
jessjaco opened this issue Aug 15, 2024 · 1 comment
Open

Data gap when loading items across the antimeridian #165

jessjaco opened this issue Aug 15, 2024 · 1 comment

Comments

@jessjaco
Copy link

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
Copy link
Contributor

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants