odc-geo
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`
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.