Track CRS info in BoundingBox and changes to GeoBox constructors

docs/geobox.rst

@@ -5,22 +5,28 @@ GeoBox
 .. autosummary::
    :toctree: _api/
 
+   GeoBox.from_geopolygon
+   GeoBox.from_bbox
    GeoBox.affine
-   GeoBox.alignment
    GeoBox.boundary
    GeoBox.buffered
    GeoBox.coordinates
-   GeoBox.coords
    GeoBox.crs
    GeoBox.dimensions
    GeoBox.dims
-   GeoBox.extent
    GeoBox.flipx
    GeoBox.flipy
-   GeoBox.from_geopolygon
+   GeoBox.extent
+   GeoBox.footprint
+   GeoBox.boundingbox
    GeoBox.geographic_extent
+   GeoBox.overlap_roi
    GeoBox.height
    GeoBox.is_empty
+   GeoBox.left
+   GeoBox.right
+   GeoBox.top
+   GeoBox.bottom
    GeoBox.pad
    GeoBox.pad_wh
    GeoBox.resolution

odc/geo/geobox.py

@@ -5,13 +5,14 @@
 import itertools
 import math
 from collections import OrderedDict, namedtuple
+from enum import Enum
 from textwrap import dedent
 from typing import Dict, Iterable, List, Literal, Optional, Tuple, Union
 
 import numpy
 from affine import Affine
 
-from .crs import CRS, MaybeCRS, norm_crs
+from .crs import CRS, MaybeCRS, SomeCRS, norm_crs
 from .geom import (
     BoundingBox,
     Geometry,
@@ -23,13 +24,14 @@
     point,
     polygon_from_transform,
 )
-from .math import clamp, is_affine_st, is_almost_int
+from .math import clamp, is_affine_st, is_almost_int, snap_grid
 from .roi import Tiles as RoiTiles
 from .roi import align_up, polygon_path, roi_normalise, roi_shape
 from .types import (
     XY,
     Index2d,
     MaybeInt,
+    NormalizedROI,
     Resolution,
     Shape2d,
     SomeIndex2d,
@@ -46,27 +48,21 @@
     Literal["native"], Literal["pixel"], Literal["geo"], Literal["auto"]
 ]
 
-# pylint: disable=invalid-name,too-many-public-methods,too-many-lines
-Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
 
+class AnchorEnum(Enum):
+    """
+    Defines which way to snap geobox pixel grid.
+    """
 
-def _align_pix(left: float, right: float, res: float, off: float) -> Tuple[float, int]:
-    if res < 0:
-        res = -res
-        val = math.ceil((right - off) / res) * res + off
-        width = max(1, int(math.ceil((val - left - 0.1 * res) / res)))
-    else:
-        val = math.floor((left - off) / res) * res + off
-        width = max(1, int(math.ceil((right - val - 0.1 * res) / res)))
-    return val, width
+    EDGE = 0
+    CENTER = 1
+    FLOATING = 2
 
 
-def _align_pix_tight(left: float, right: float, res: float) -> Tuple[float, int]:
-    if res < 0:
-        res = -res
-        return right, max(1, int(math.ceil((right - left - 0.1 * res) / res)))
+GeoboxAnchor = Union[AnchorEnum, XY[float]]
 
-    return left, max(1, int(math.ceil((right - left - 0.1 * res) / res)))
+# pylint: disable=invalid-name,too-many-public-methods,too-many-lines
+Coordinate = namedtuple("Coordinate", ("values", "units", "resolution"))
 
 
 class GeoBox:
@@ -97,39 +93,57 @@ def from_bbox(
         tight: bool = False,
         shape: Optional[SomeShape] = None,
         resolution: Optional[SomeResolution] = None,
-        align: Optional[XY[float]] = None,
+        anchor: GeoboxAnchor = AnchorEnum.EDGE,
+        tol: float = 0.01,
     ) -> "GeoBox":
         """
         Construct :py:class:`~odc.geo.geobox.GeoBox` from a bounding box.
 
         :param bbox: Bounding box in CRS units, lonlat is assumed when ``crs`` is not supplied
         :param crs: CRS of the bounding box (defaults to EPSG:4326)
-        :param tight: Supplying ``tight=True`` turns off pixel snapping.
         :param shape: Span that many pixels.
         :param resolution: Use specified resolution
-        :param align:
-        :return: :py:class:`~odc.geo.geobox.GeoBox` that covers supplied bounding box.
+        :param tight: Supplying ``tight=True`` turns off pixel snapping.
+        :param anchor:
+            By default snaps grid such that pixel edges fall on X/Y axis. Ignored when tight mode is
+            used.
+        :param tol:
+            Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output
+            geobox is allowed to be smaller than supplied bounding box by that amount.
+
+        :return:
+           :py:class:`~odc.geo.geobox.GeoBox` that covers supplied bounding box.
         """
+        # pylint: disable=too-many-locals
 
-        if align is None and tight is False:
-            align = xy_(0, 0)
-        if crs is None:
-            crs = "EPSG:4326"
+        _snap: Optional[XY[float]] = None
+
+        if tight:
+            anchor = AnchorEnum.FLOATING
+
+        if isinstance(anchor, XY):
+            _snap = anchor
+        if anchor == AnchorEnum.EDGE:
+            _snap = xy_(0, 0)
+        elif anchor == AnchorEnum.CENTER:
+            _snap = xy_(0.5, 0.5)
 
         if not isinstance(bbox, BoundingBox):
-            bbox = BoundingBox(*bbox)
+            bbox = BoundingBox(*bbox, crs=(crs or "epsg:4326"))
+        elif bbox.crs is None:
+            bbox = BoundingBox(*bbox.bbox, crs=(crs or "epsg:4326"))
 
         if resolution is not None:
             rx, ry = res_(resolution).xy
-            if align is None:
-                offx, nx = _align_pix_tight(bbox.left, bbox.right, rx)
-                offy, ny = _align_pix_tight(bbox.bottom, bbox.top, ry)
+            if _snap is None:
+                offx, nx = snap_grid(bbox.left, bbox.right, rx, None, tol=tol)
+                offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None, tol=tol)
             else:
-                offx, nx = _align_pix(bbox.left, bbox.right, rx, align.x)
-                offy, ny = _align_pix(bbox.bottom, bbox.top, ry, align.y)
+                offx, nx = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol)
+                offy, ny = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol)
 
             affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
-            return GeoBox((ny, nx), crs=crs, affine=affine)
+            return GeoBox((ny, nx), crs=bbox.crs, affine=affine)
 
         if shape is None:
             raise ValueError("Must supply shape or resolution")
@@ -139,21 +153,22 @@ def from_bbox(
         rx = bbox.span_x / nx
         ry = -bbox.span_y / ny
 
-        if align is None:
+        if _snap is None:
             offx, offy = bbox.left, bbox.top
         else:
-            offx, _ = _align_pix(bbox.left, bbox.right, rx, align.x)
-            offy, _ = _align_pix(bbox.bottom, bbox.top, ry, align.y)
+            offx, _ = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol)
+            offy, _ = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol)
 
         affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
-        return GeoBox((ny, nx), crs=crs, affine=affine)
+        return GeoBox((ny, nx), crs=bbox.crs, affine=affine)
 
     @staticmethod
     def from_geopolygon(
         geopolygon: Geometry,
         resolution: SomeResolution,
         crs: MaybeCRS = None,
-        align: Optional[XY[float]] = None,
+        anchor: GeoboxAnchor = AnchorEnum.EDGE,
+        tol: float = 0.01,
     ) -> "GeoBox":
         """
         Construct :py:class:`~odc.geo.geobox.GeoBox` from a polygon.
@@ -162,27 +177,23 @@ def from_geopolygon(
            Either a single number or a :py:class:`~odc.geo.types.Resolution` object.
         :param crs:
            CRS to use, if different from the geopolygon
-        :param align:
-           Align geobox such that point 'align' lies on the pixel boundary.
-        """
-        resolution = res_(resolution)
-        if align is None:
-            align = xy_(0.0, 0.0)
 
-        assert (
-            0.0 <= align.x <= abs(resolution.x)
-        ), "X align must be in [0, abs(x_resolution)] range"
-        assert (
-            0.0 <= align.y <= abs(resolution.y)
-        ), "Y align must be in [0, abs(y_resolution)] range"
+        :param anchor:
+            By default snaps grid such that pixel edges fall on X/Y axis.
 
+         :param tol:
+            Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output
+            geobox is allowed to be smaller than supplied bounding box by that amount.
+
+        """
+        resolution = res_(resolution)
         if crs is None:
             crs = geopolygon.crs
         else:
             geopolygon = geopolygon.to_crs(crs)
 
         return GeoBox.from_bbox(
-            geopolygon.boundingbox, crs, resolution=resolution, align=align
+            geopolygon.boundingbox, crs, resolution=resolution, anchor=anchor, tol=tol
         )
 
     def buffered(self, xbuff: float, ybuff: Optional[float] = None) -> "GeoBox":
@@ -246,6 +257,22 @@ def __bool__(self) -> bool:
     def __hash__(self):
         return hash((*self._shape, self._crs, self._affine))
 
+    def overlap_roi(self, other: "GeoBox", tol: float = 1e-8) -> NormalizedROI:
+        """
+        Compute overlap as ROI.
+
+        Figure out slice into this geobox that shares pixels with the ``other`` geobox with
+        consistent pixel grid.
+
+        :raises:
+            :py:class:`ValueError` when two geoboxes are not pixel-aligned.
+        """
+        nx, ny = self._shape.xy
+        x0, y0, x1, y1 = map(int, bounding_box_in_pixel_domain(other, self, tol))
+        x0, y0 = max(0, x0), max(0, y0)
+        x1, y1 = min(x1, nx), min(y1, ny)
+        return numpy.s_[y0:y1, x0:x1]
+
     @property
     def transform(self) -> Affine:
         """Linear mapping from pixel space to CRS."""
@@ -342,20 +369,41 @@ def extent(self) -> Geometry:
         self._extent = _extent
         return _extent
 
+    @property
+    def boundingbox(self) -> BoundingBox:
+        """GeoBox bounding box in the native CRS."""
+        return BoundingBox.from_transform(self._shape, self._affine, crs=self._crs)
+
     def _reproject_resolution(self, npoints: int = 100):
         bbox = self.extent.boundingbox
         span = max(bbox.span_x, bbox.span_y)
         return span / npoints
 
+    def footprint(
+        self, crs: SomeCRS, buffer: float = 0, npoints: int = 100
+    ) -> Geometry:
+        """
+        Compute footprint in foreign CRS.
+
+        :param crs: CRS of the destination
+        :param buffer: amount to buffer in source pixels before transforming
+        :param npoints: number of points per-side to use, higher number
+                        is slower but more accurate
+        """
+        assert self.crs is not None
+        ext = self.extent
+        if buffer > 0:
+            buffer = buffer * max(*self.resolution.xy)
+            ext = ext.buffer(buffer)
+
+        return ext.to_crs(crs, resolution=self._reproject_resolution(npoints))
+
     @property
     def geographic_extent(self) -> Geometry:
         """GeoBox extent in EPSG:4326."""
         if self._crs is None or self._crs.geographic:
             return self.extent
-
-        return self.extent.to_crs(
-            CRS("EPSG:4326"), resolution=self._reproject_resolution()
-        )
+        return self.footprint("epsg:4326")
 
     coords = coordinates
     dims = dimensions
@@ -800,7 +848,7 @@ def bounding_box_in_pixel_domain(
 
     tx, ty = round(tx), round(ty)
     ny, nx = geobox.shape
-    return BoundingBox(tx, ty, tx + nx, ty + ny)
+    return BoundingBox(tx, ty, tx + nx, ty + ny, None)
 
 
 def geobox_union_conservative(geoboxes: List[GeoBox]) -> GeoBox:
@@ -840,11 +888,19 @@ def geobox_intersection_conservative(geoboxes: List[GeoBox]) -> GeoBox:
     # standardise empty geobox representation
     if bbox.left > bbox.right:
         bbox = BoundingBox(
-            left=bbox.left, bottom=bbox.bottom, right=bbox.left, top=bbox.top
+            left=bbox.left,
+            bottom=bbox.bottom,
+            right=bbox.left,
+            top=bbox.top,
+            crs=bbox.crs,
         )
     if bbox.bottom > bbox.top:
         bbox = BoundingBox(
-            left=bbox.left, bottom=bbox.bottom, right=bbox.right, top=bbox.bottom
+            left=bbox.left,
+            bottom=bbox.bottom,
+            right=bbox.right,
+            top=bbox.bottom,
+            crs=bbox.crs,
         )
 
     affine = reference.affine * Affine.translation(*bbox[:2])