odc-geo icon indicating copy to clipboard operation
odc-geo copied to clipboard

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

Open alexgleith opened this issue 5 months ago • 1 comments

The issue is documented in detail in https://github.com/opendatacube/odc-stac/issues/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 avatar Sep 06 '24 04:09 alexgleith