diff --git a/odc/geo/crs.py b/odc/geo/crs.py index 8553b69c..5401af0f 100644 --- a/odc/geo/crs.py +++ b/odc/geo/crs.py @@ -76,12 +76,12 @@ def _make_crs( return (crs, crs_str, epsg) -def _make_crs_transform_key(from_crs, to_crs, always_xy): - return (id(from_crs), id(to_crs), always_xy) +def _make_crs_transform_key(from_crs, to_crs, always_xy, force_over=False): + return (id(from_crs), id(to_crs), always_xy, force_over) @cachetools.cached({}, key=_make_crs_transform_key) -def _make_crs_transform(from_crs: _CRS, to_crs: _CRS, always_xy: bool) -> Transformer: +def _make_crs_transform(from_crs: _CRS, to_crs: _CRS, always_xy: bool, force_over: bool = False) -> Transformer: return Transformer.from_crs(from_crs, to_crs, always_xy=always_xy, force_over=force_over) @@ -299,7 +299,7 @@ def crs_str(self) -> str: return self._str def transformer_to_crs( - self, other: "CRS", always_xy: bool = True + self, other: "CRS", always_xy: bool = True, force_over: bool = False ) -> Callable[[Any, Any], Tuple[Any, Any]]: """ Build coordinate transformer to other projection. @@ -318,7 +318,7 @@ def transformer_to_crs( """ # pylint: disable=protected-access - tr = _make_crs_transform(self._crs, other._crs, always_xy=always_xy) + tr = _make_crs_transform(self._crs, other._crs, always_xy=always_xy, force_over=force_over) def result(x, y, **kw): rx, ry = tr.transform(x, y, **kw) # pylint: disable=unpacking-non-sequence diff --git a/odc/geo/overlap.py b/odc/geo/overlap.py index 759b9116..036bb77d 100644 --- a/odc/geo/overlap.py +++ b/odc/geo/overlap.py @@ -103,7 +103,7 @@ def __init__( self._src = src self._dst = dst self._back = back - self._tr = src.crs.transformer_to_crs(dst.crs) + self._tr = src.crs.transformer_to_crs(dst.crs, force_over=True) self._clamps: Optional[Tuple[Tuple[float, float], Tuple[float, float]]] = None if src.crs.geographic: self._clamps = ((-180, 180), (-90, 90)) diff --git a/tests/test_geom.py b/tests/test_geom.py index 7ccea762..565e2ea2 100644 --- a/tests/test_geom.py +++ b/tests/test_geom.py @@ -522,6 +522,53 @@ def test_gen_test_image_xy(): assert isinstance(A, Affine) +def test_geobox_overlap(): + from odc.geo.types import xy_ + 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") + + 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)]) + + # Results should be within a resonable range in pixel space + # Not sure how to test it better. + for xy in xys: + assert xy.x >= 0 - 25.6 + assert xy.y <= 256 + 25.6 + + def test_fixed_point(): aa = np.asarray([0, 0.5, 1]) uu = to_fixed_point(aa, "uint8")