Inability to load pixels near/at antimeridian

[alexgleith]

We've been working to try to resolve issue with loading data at the antimeridian.

Currently, there seems to be an issue in both Datacube Core and ODC STAC when loading an antimeridian-touching geobox.

Below is some code that reproduces the issue.

from affine import Affine
from pystac import Item
from odc.geo.geobox import GeoBox
from odc.stac import load

affine = Affine(
    152.87405657034833,
    0.0,
    -20037508.342789244,
    0.0,
    -152.87405657034833,
    -1995923.6825825237,
)
geobox = GeoBox((256, 256), affine, "EPSG:3857")

item = Item.from_file("https://earth-search.aws.element84.com/v1/collections/sentinel-2-c1-l2a/items/S2A_T01KAA_20240825T221936_L2A")
data = load([item], geobox=geobox, measurements=["red", "green", "blue"], chunks={})

data.squeeze().to_array().plot.imshow()

The geobox looks fine, and is directly adjacent to the antimeridian.

The rendered data (white) should go all the way to the left edge of the frame in this image below.

[SpacemanPaul]

Something similar happens for geoboxes abutting the antimeridian from the other side (i.e. with the anti-meridian along the right edge).

[alexgleith]

Ok, some kind of minimal reproduction of the cause here:

from odc.geo.types import xy_
import numpy as np
from pyproj import Transformer
from odc.geo.geobox import GeoBox
from affine import Affine

from odc.geo.overlap import roi_boundary, unstack_xy, stack_xy, gbox_boundary, roi_from_points, native_pix_transform

pts_per_side = 5
padding = 1
align = True

dst_affine = Affine(
    152.87405657034833,
    0.0,
    -20037508.342789244,
    0.0,
    -152.87405657034833,
    -1995923.6825825237,
)
dst = GeoBox((256, 256), dst_affine, "EPSG:3857")

src_affine = Affine(10.0, 0.0, 99960.0, 0.0, -10.0, 8100040.0)
src = GeoBox((10980, 10980), src_affine, "EPSG:32701")

# Do this the internal odc.geo.overlay way
# NOTE: These functions hide the world/pixel conversion.
tr = native_pix_transform(src, dst)

xy = tr.back(unstack_xy(gbox_boundary(dst, pts_per_side)))
roi_src = roi_from_points(stack_xy(xy), src.shape, padding, align=align)

xy_pix_src = unstack_xy(roi_boundary(roi_src, pts_per_side))

xx, yy = np.asarray([pt.xy for pt in xy_pix_src]).T
xx,yy

# This goes via world transform to a pixel space
xys = tr([xy_(x, y) for x, y in zip(xx, yy)])
xys

[XY(x=262143.98348895763, y=-3.9324825294115726), XY(x=48.16012841704651, y=-3.171199267373595), XY(x=96.33989687502617, y=-2.4268339181071497), XY(x=144.52272405748954, y=-1.699391535272298), ...

that first coord is wrong, clearly. It should be near 0

And deeper simpler reproduction:

transformer_to = Transformer.from_crs(src.crs, dst.crs, always_xy=True)

# xy_pix_src is in pixel space of source
xx, yy = np.asarray([pt.xy for pt in xy_pix_src]).T

# Convert to world coordinates
xx, yy = src.pix2wld(xx, yy)

# Transform to destination crs
xx_p, yy_p = transformer_to.transform(xx, yy, errcheck=True)

# print(f"This should be negative {xx_p[0]}")

# Uncommenting this fixes the whole thing, but it's not the right way to do it
# for i, n in enumerate(xx_p):
#     if n > 0:
#         xx_p[i] = n * -1

# Convert back to pixel space of destination
# NOTE: this is the step that is resulting in incorrect coordinates
xx, yy = dst.wld2pix(*(xx_p, yy_p))
xys_new = [xy_(x, y) for x, y in zip(xx, yy)]

# Coords are borked
xys_new

[XY(x=262143.98348895763, y=-3.9324825294115726), XY(x=48.16012841704651, y=-3.171199267373595), XY(x=96.33989687502617, y=-2.4268339181071497) ...

Note that fixing the broken coordinates coming out of the transformer fixes everything:

# for i, n in enumerate(xx_p):
#     if n > 0:
#         xx_p[i] = n * -1

[alexgleith]

At the core, the issue is caused by code in odc-geo but I'll leave this here for now.

[alexgleith]

This code is where the breakage is: https://github.com/opendatacube/odc-geo/blob/94ad126571cdddcf4884b0f9ea4024dde15d0208/odc/geo/overlap.py#L138

And I can't think of any way to fix the issue!