Replace alignment with anchor

`align=` was too confusing and hard to document, concept of anchor I
feel is more reasonable. You can anchor pixel edges to 0,0 or pixel
center to 0,0, or you can leave it floating. Or you can pick any
location within a pixel and say that 0,0 should coincide with that.

center == (0.5, 0.5)
edge == (0, 0)

Also replacing _align_pix methods, that were hard to decipher, with
snap_grid method in odc.geo.math

odc/geo/geobox.py

@@ -5,6 +5,7 @@
 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
 
@@ -23,7 +24,7 @@
  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 (
@@ -47,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:
@@ -98,22 +93,35 @@ def from_bbox(
  tight: bool = False,
  shape: Optional[SomeShape] = None,
  resolution: Optional[SomeResolution] = None,
- align: Optional[XY[float]] = None,
+ anchor: GeoboxAnchor = AnchorEnum.EDGE,
  ) -> "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 pixel edges are align with axis. Ignored when tight mode is
+ used.
+
+ :return:
+ :py:class:`~odc.geo.geobox.GeoBox` that covers supplied bounding box.
  """
 
- if align is None and tight is False:
- align = xy_(0, 0)
+ _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, crs=(crs or "epsg:4326"))
@@ -122,12 +130,12 @@ def from_bbox(
 
  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)
+ offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None)
  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)
+ offy, ny = snap_grid(bbox.bottom, bbox.top, ry, _snap.y)
 
  affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
  return GeoBox((ny, nx), crs=bbox.crs, affine=affine)
@@ -140,11 +148,11 @@ 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)
+ offy, _ = snap_grid(bbox.bottom, bbox.top, ry, _snap.y)
 
  affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
  return GeoBox((ny, nx), crs=bbox.crs, affine=affine)
@@ -154,7 +162,7 @@ def from_geopolygon(
  geopolygon: Geometry,
  resolution: SomeResolution,
  crs: MaybeCRS = None,
- align: Optional[XY[float]] = None,
+ anchor: GeoboxAnchor = AnchorEnum.EDGE,
  ) -> "GeoBox":
  """
  Construct :py:class:`~odc.geo.geobox.GeoBox` from a polygon.
@@ -163,27 +171,15 @@ 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"
-
  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
  )
 
  def buffered(self, xbuff: float, ybuff: Optional[float] = None) -> "GeoBox":

odc/geo/math.py

@@ -7,7 +7,7 @@
 
 Minimal dependencies in this module.
 """
-from math import floor, fmod
+from math import ceil, floor, fmod
 from typing import Literal, Optional, Tuple, Union
 
 import numpy as np
@@ -123,6 +123,53 @@ def is_almost_int(x: float, tol: float) -> bool:
  return x < tol
 
 
+def _snap_edge_pos(x0: float, x1: float, res: float) -> Tuple[float, int]:
+ assert res > 0
+ assert x1 >= x0
+ _x0, _x1 = floor(x0 / res), ceil(x1 / res)
+ nx = max(1, _x1 - _x0)
+ return _x0 * res, nx
+
+
+def _snap_edge(x0: float, x1: float, res: float) -> Tuple[float, int]:
+ assert x1 >= x0
+ if res > 0:
+ return _snap_edge_pos(x0, x1, res)
+ _tx, nx = _snap_edge_pos(x0, x1, -res)
+ tx = _tx + nx * (-res)
+ return tx, nx
+
+
+def snap_grid(
+ x0: float, x1: float, res: float, off_pix: Optional[float] = 0
+) -> Tuple[float, int]:
+ """
+ Compute grid snapping for single axis.
+
+ :param x0: In point ``x0 <= x1``
+ :param x1: Out point ``x0 <= x1``
+ :param res: Pixel size and direction (can be negative)
+ :param off_pix:
+ Pixel fraction to align to ``x=0``.
+ 0 - edge aligned
+ 0.5 - center aligned
+ None - don't snap
+
+ :return: ``tx, nx`` that defines 1-d grid, such that ``x0`` and ``x1`` are within edge pixels.
+ """
+ assert (off_pix is None) or (0 <= off_pix < 1)
+ if off_pix is None:
+ if res > 0:
+ nx = ceil((x1 - x0) / res)
+ return x0, max(1, nx)
+ nx = ceil((x1 - x0) / (-res))
+ return x1, max(nx, 1)
+
+ off = off_pix * abs(res)
+ _tx, nx = _snap_edge(x0 - off, x1 - off, res)
+ return _tx + off, nx
+
+
 def data_resolution_and_offset(
  data, fallback_resolution: Optional[float] = None
 ) -> Tuple[float, float]:

tests/test_geobox.py

@@ -7,8 +7,9 @@
 import pytest
 from affine import Affine
 
-from odc.geo import CRS, geom, ixy_, resyx_, wh_
+from odc.geo import CRS, geom, ixy_, resyx_, wh_, xy_
 from odc.geo.geobox import (
+ AnchorEnum,
  GeoBox,
  bounding_box_in_pixel_domain,
  gbox_boundary,
@@ -308,6 +309,18 @@ def test_from_bbox():
  == "epsg:3857"
  )
 
+ assert GeoBox.from_bbox(
+ bbox, shape=shape, anchor=AnchorEnum.FLOATING
+ ) == GeoBox.from_bbox(bbox, shape=shape, tight=True)
+
+ assert GeoBox.from_bbox(
+ bbox, shape=shape, anchor=AnchorEnum.EDGE
+ ) == GeoBox.from_bbox(bbox, shape=shape, anchor=xy_(0, 0))
+
+ assert GeoBox.from_bbox(
+ bbox, shape=shape, anchor=AnchorEnum.CENTER
+ ) == GeoBox.from_bbox(bbox, shape=shape, anchor=xy_(0.5, 0.5))
+
  # one of resolution= or shape= must be supplied
  with pytest.raises(ValueError):
  _ = GeoBox.from_bbox(bbox)

tests/test_geom.py

@@ -13,7 +13,7 @@
 
 from odc.geo import CRS, BoundingBox, CRSMismatchError, geom, wh_
 from odc.geo.crs import norm_crs, norm_crs_or_error
-from odc.geo.geobox import GeoBox, _align_pix, _round_to_res
+from odc.geo.geobox import GeoBox, _round_to_res
 from odc.geo.geom import (
  bbox_intersection,
  bbox_union,
@@ -827,25 +827,6 @@ def test_geom_clone():
  assert b.geom is not geom.Geometry(b).geom
 
 
-@pytest.mark.parametrize(
- "left, right, res, off, expect",
- [
- (20, 30, 10, 0, (20, 1)),
- (20, 30.5, 10, 0, (20, 1)),
- (20, 31.5, 10, 0, (20, 2)),
- (20, 30, 10, 3, (13, 2)),
- (20, 30, 10, -3, (17, 2)),
- (20, 30, -10, 0, (30, 1)),
- (19.5, 30, -10, 0, (30, 1)),
- (18.5, 30, -10, 0, (30, 2)),
- (20, 30, -10, 3, (33, 2)),
- (20, 30, -10, -3, (37, 2)),
- ],
-)
-def test_align_pix(left, right, res, off, expect):
- assert _align_pix(left, right, res, off) == expect
-
-
 def test_lonlat_bounds():
  # example from landsat scene: spans lon=180
  poly = geom.box(618300, -1876800, 849000, -1642500, "EPSG:32660")

tests/test_math.py

@@ -16,6 +16,7 @@
  maybe_int,
  maybe_zero,
  snap_affine,
+ snap_grid,
  snap_scale,
  split_translation,
 )
@@ -216,3 +217,20 @@ def test_snap_affine():
  ttol=0.2,
  stol=1e-3,
  ) == mkA(scale=(1 / 2, 1 / 3), translation=(10, 20))
+
+
+@pytest.mark.parametrize(
+ "left, right, res, off, expect",
+ [
+ (20, 30, 10, 0, (20, 1)),
+ (20, 30.5, 10, 0, (20, 2)),
+ (20, 31.5, 10, 0, (20, 2)),
+ (20, 30, 10, 3, (13, 2)),
+ (20, 30, -10, 0, (30, 1)),
+ (19.5, 30, -10, 0, (30, 2)),
+ (18.5, 30, -10, 0, (30, 2)),
+ (20, 30, -10, 3, (33, 2)),
+ ],
+)
+def test_snap_grid(left, right, res, off, expect):
+ assert snap_grid(left, right, res, off / abs(res)) == expect