Some cases of geometries that are near and that cross the antimeridian don't work well in odc.geo.overlap._relative_rois

[alexgleith]

The issue is documented in detail in opendatacube/odc-stac#172

The short reproduction of the cause of the issue is this code:

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

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

which results in:

[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),
 XY(x=192.70853967263247, y=-0.9888770583802398),
 XY(x=191.77107980153232, y=63.979562991820785),
 XY(x=190.82837076691794, y=128.97819103086113),
 XY(x=189.88040104990068, y=194.00717585515667),
 XY(x=188.9271590545104, y=259.0666866327829),
 XY(x=140.6501542032056, y=258.34019569222437),
 XY(x=92.37616415832599, y=257.59639652622354),
 XY(x=44.105259828415, y=256.8352942077672),
 XY(x=262139.83751210378, y=256.05689392704153),
 XY(x=262140.88266387527, y=191.01398425493971),
 XY(x=262141.92203548463, y=126.00156417726794),
 XY(x=262142.95563963818, y=61.01946476773992)]

Noting that the large value coords like 262143.xxx should be something like 0.xxx.

The main issue is transforming from Web Mercator to UTM and back to to Web Mercator, which results in ambiguity, and therefore, coordinates that are one on side or the other of the antimeridian.

[alexgleith]

This is the more explicit code illustrating the issue:

transformer_to = Transformer.from_crs(src.crs, dst.crs, always_xy=True)
t_to_wgs = Transformer.from_crs(src.crs, "EPSG:4326", 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, which results in some coords that are on the wrong
# side of the antimeridian
xx_p, yy_p = transformer_to.transform(xx, yy)

# 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
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

[alexgleith]

Another example of unloadable data at the antimeridian:

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

from odc.stac import load

url = 'https://landsatlook.usgs.gov/stac-server/collections/landsat-c2l2-sr/items/LC09_L2SR_073073_20240128_20240130_02_T1_SR'
item = Item.from_file(url)

src_geobox = GeoBox(item.properties["proj:shape"], Affine(*item.properties["proj:transform"]), crs=item.properties["proj:epsg"])

# This is wrong... it's world spanning.
src_geobox.geographic_extent

shape = (5000, 5000)
transform = Affine(0.0003, 0.0, 180.0, 0.0, -0.0003, -18.0)
crs = 'EPSG:4326'

dst_geobox = GeoBox(shape, transform, crs)

# Can't load data in the target geobox
data = load([item], chunks={}, geobox=dst_geobox, bands=["red"]).compute()
data  # Empty