Repository: opendatacube/odc-geo
Pull request #148 (2024-05-07 04:44:57)

Geobox anchor accessor
- Adding  `GeoBox.anchor`
- Some fixes in COG related tests

github-actions[bot]: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: 2516d918dc82767c089255cf7129f75d9f959e53c07db8d3368972c4656b3dce -->
🚀 Deployed on https://6639b259692e0894fdf68101--odc-geo-docs.netlify.app

Patch:
diff --git a/docs/api-geobox.rst b/docs/api-geobox.rst
index dc3a646..446376f 100644
--- a/docs/api-geobox.rst
+++ b/docs/api-geobox.rst
@@ -9,6 +9,7 @@ GeoBox
    GeoBox.from_geopolygon
    GeoBox.from_bbox
    GeoBox.affine
+   GeoBox.anchor
    GeoBox.boundary
    GeoBox.buffered
    GeoBox.coordinates
diff --git a/odc/geo/geobox.py b/odc/geo/geobox.py
index 5e9433b..0addf16 100644
--- a/odc/geo/geobox.py
+++ b/odc/geo/geobox.py
@@ -2,6 +2,8 @@
 #
 # Copyright (c) 2015-2020 ODC Contributors
 # SPDX-License-Identifier: Apache-2.0
+from __future__ import annotations
+
 import importlib
 import itertools
 import math
@@ -28,6 +30,7 @@ from .crs import CRS, MaybeCRS, SomeCRS, norm_crs
 from .geom import BoundingBox, Geometry, bbox_intersection, bbox_union
 from .math import (
     clamp,
+    extract_anchor,
     is_affine_st,
     is_almost_int,
     maybe_zero,
@@ -195,10 +198,24 @@ class GeoBoxBase:
         Alignment of pixel boundaries in CRS units.
 
         This is usally ``(0,0)``.
+
+        .. seealso:: :py:meth:`~odc.geo.geobox.GeoBoxBase.anchor`
         """
+        # TODO: deprecate alignment method
         rx, _, tx, _, ry, ty, *_ = self._affine
         return xy_(tx % abs(rx), ty % abs(ry))
 
+    @property
+    def anchor(self) -> XY[float] | AnchorEnum:
+        """
+        Anchor point in pixel coordinates.
+
+        This is only valid for linear sources.
+        """
+        if not self.linear:
+            return AnchorEnum.FLOATING
+        return extract_anchor(self._affine)
+
     @property
     def linear(self) -> bool:
         return True
diff --git a/odc/geo/math.py b/odc/geo/math.py
index 92c7da1..74ca86e 100644
--- a/odc/geo/math.py
+++ b/odc/geo/math.py
@@ -7,6 +7,8 @@ Various mathy helpers.
 
 Minimal dependencies in this module.
 """
+from __future__ import annotations
+
 from math import ceil, floor, fmod, isfinite, log2
 from typing import (
     Any,
@@ -24,7 +26,17 @@ import numpy as np
 from affine import Affine
 from numpy.polynomial.polynomial import polygrid2d, polyval2d
 
-from .types import XY, Resolution, SomeResolution, SomeShape, res_, resxy_, shape_, xy_
+from .types import (
+    XY,
+    AnchorEnum,
+    Resolution,
+    SomeResolution,
+    SomeShape,
+    res_,
+    resxy_,
+    shape_,
+    xy_,
+)
 
 AffineX = TypeVar("AffineX", np.ndarray, Affine)
 
@@ -795,3 +807,18 @@ class Poly2d:
         cc = np.vstack([cc, np.asarray([0, 0])]).reshape(2, 2, 2)
 
         return Poly2d(cc, Ain)
+
+
+def extract_anchor(pix2wld: Affine, tol: float = 1e-3) -> AnchorEnum | XY[float]:
+    def _anchor(px: float, tol: float) -> float:
+        _, x = split_float(px)  # x in (-0.5, +0.5)
+        x = (1 + x) if x < 0 else x  # x in [0, 1)
+        x = snap_scale(maybe_zero(x, tol), tol)
+        return 0 if x >= 1 else x
+
+    anchor = tuple(_anchor(px, tol) for px in (~pix2wld) * (0, 0))
+    if anchor == (0, 0):
+        return AnchorEnum.EDGE
+    if anchor == (0.5, 0.5):
+        return AnchorEnum.CENTER
+    return xy_(anchor)


Test patch:
diff --git a/tests/test_cog.py b/tests/test_cog.py
index db740d0..4f2a23e 100644
--- a/tests/test_cog.py
+++ b/tests/test_cog.py
@@ -5,6 +5,10 @@ from typing import Optional, Tuple
 
 import numpy as np
 import pytest
+
+_ = pytest.importorskip("tifffile")
+_ = pytest.importorskip("rasterio")
+
 from dask import array as da
 from rasterio import MemoryFile
 from rasterio import open as rio_open
@@ -543,12 +547,15 @@ def test_cog_with_dask_smoke_test(gbox: GeoBox, tmp_path: Path, dtype):
             id="float with nan, nan nodata",
         ),
         pytest.param(
-            [[1.0, np.nan], [0, 1.0]],
+            [
+                [1.0, np.nan],
+                [0.0, 9.0],
+            ],
             0,
             1,
-            1,
-            1,
-            0,
+            9.0,
+            5.0,
+            4.0,
             50,
             id="float with nan, numeric nodata",
         ),
@@ -559,7 +566,7 @@ def test_stats_from_layer(array, nodata, minimum, maximum, mean, stddev, valid_p
     stats = _stats_from_layer(x, nodata).compute()[0]
 
     assert stats["minimum"] == minimum
-    assert stats["maximum"] == minimum
+    assert stats["maximum"] == maximum
     assert stats["mean"] == mean
     assert stats["stddev"] == stddev
     assert stats["valid_percent"] == valid_percent
diff --git a/tests/test_gcp.py b/tests/test_gcp.py
index 8dc398f..3d263cd 100644
--- a/tests/test_gcp.py
+++ b/tests/test_gcp.py
@@ -6,7 +6,7 @@ import pytest
 
 from odc.geo import geom
 from odc.geo.gcp import GCPGeoBox, GCPMapping
-from odc.geo.geobox import GeoBox
+from odc.geo.geobox import AnchorEnum, GeoBox
 from odc.geo.xr import xr_zeros
 
 rasterio = pytest.importorskip("rasterio")
@@ -105,6 +105,8 @@ def test_gcp_geobox_basics(au_gcp_geobox: GCPGeoBox):
     assert p1 == pytest.approx((-44.50301336231415, 109.39806656168265))
     assert p2 == pytest.approx((-9.47177497427409, 157.04711254391185))
 
+    assert gbox.anchor == AnchorEnum.FLOATING
+
 
 def test_gcp_geobox_xr(au_gcp_geobox: GCPGeoBox):
     gbox = au_gcp_geobox
diff --git a/tests/test_geobox.py b/tests/test_geobox.py
index 2a6b212..2a6b45e 100644
--- a/tests/test_geobox.py
+++ b/tests/test_geobox.py
@@ -7,8 +7,8 @@ import numpy as np
 import pytest
 from affine import Affine
 
-from odc.geo._interop import have
 from odc.geo import CRS, geom, ixy_, resxy_, resyx_, wh_, xy_
+from odc.geo._interop import have
 from odc.geo.geobox import (
     AnchorEnum,
     GeoBox,
@@ -666,7 +666,7 @@ def test_crop(gbox, shape):
     ],
 )
 def test_explore_geobox(geobox):
-    from folium import Map, GeoJson
+    from folium import GeoJson, Map
 
     # Test explore on dataset input and verify that output is a folium map
     # that contains a GeoJson layer
@@ -704,3 +704,34 @@ def test_shape(shape, allowed):
     else:
         with pytest.raises(ValueError):
             GeoBox.from_bbox((-3, -4, 3, 4), epsg4326, shape=shape)
+
+
+@pytest.mark.parametrize(
+    "gbox, expected_anchor",
+    [
+        (
+            GeoBox.from_bbox((3, 4, 5, 7), epsg4326, resolution=0.25, anchor="edge"),
+            AnchorEnum.EDGE,
+        ),
+        (
+            GeoBox.from_bbox((-3, -4, 5, 7), epsg4326, resolution=0.25),
+            AnchorEnum.EDGE,
+        ),
+        (
+            GeoBox.from_bbox(
+                (-3, -4, 5, 7), epsg4326, resolution=0.25, anchor="center"
+            ),
+            AnchorEnum.CENTER,
+        ),
+        (
+            GeoBox.from_bbox((35, 40, 50, 70), epsg3857, resolution=10, tight=True),
+            xy_(0.5, 0.0),
+        ),
+        (
+            GeoBox.from_bbox((-31, 40, 50, 72), epsg3857, resolution=10, tight=True),
+            xy_(0.1, 0.2),
+        ),
+    ],
+)
+def test_geobox_anchor(gbox, expected_anchor):
+    assert gbox.anchor == expected_anchor