geotiff: validate coord regularity in _coords_to_transform (#1720)

`_coords_to_transform` used `x[1] - x[0]` and `y[1] - y[0]` as the
affine pixel sizes without checking that the rest of the spacing
matched. GeoTIFF only supports an affine transform, so non-uniform
coords cannot be expressed faithfully; the writer silently used the
first-step spacing and produced a wrong transform.

Validate `np.diff` against the median step with a 1e-6 relative
tolerance and raise `ValueError` naming the offending dim and the
worst residual. Median (not mean) keeps a single bad sample from
shifting the reference step.

Also fix `test_pixel_is_point_round_trip`, whose hand-typed coords
had ~3.6e-3 relative deviation and were exactly the silent-wrong-
transform case the new check now catches.

Closes #1720

Author: brendancol

xrspatial/geotiff/__init__.py

@@ -299,6 +299,30 @@ def _coords_to_transform(da: xr.DataArray) -> GeoTransform | None:
     if len(x) < 2 or len(y) < 2:
         return None
 
+    # GeoTIFF only supports an affine transform; non-uniform spacing
+    # cannot be expressed faithfully. Validate up-front instead of
+    # silently writing a transform that only matches the first step.
+    def _is_regular(coord, name):
+        diffs = np.diff(coord)
+        # Use median (not mean) so a single bad sample doesn't shift
+        # the reference step. The 1e-6 relative tolerance is forgiving
+        # for float artifacts in otherwise-uniform coords.
+        step = float(np.median(diffs))
+        if step == 0:
+            raise ValueError(
+                f"{name} coords are constant; cannot infer pixel size"
+            )
+        rel = float(np.max(np.abs(diffs - step)) / abs(step))
+        if rel > 1e-6:
+            raise ValueError(
+                f"{name} coords are not uniformly spaced "
+                f"(max relative deviation {rel:.3e} exceeds 1e-6); "
+                f"GeoTIFF requires an affine transform."
+            )
+
+    _is_regular(x, "x")
+    _is_regular(y, "y")
+
     pixel_width = float(x[1] - x[0])
     pixel_height = float(y[1] - y[0])
 

xrspatial/geotiff/tests/test_coord_regularity_1720.py

@@ -0,0 +1,106 @@
+"""Regression test for issue #1720.
+
+``_coords_to_transform`` previously read ``x[1] - x[0]`` and
+``y[1] - y[0]`` as the pixel sizes without checking that the rest of
+the spacing matched. GeoTIFF supports only an affine transform, so
+non-uniform coords cannot be expressed faithfully; the writer would
+silently use the first-step spacing and produce a wrong transform.
+
+The fix validates ``np.diff(x)`` and ``np.diff(y)`` against a relative
+tolerance of ``1e-6`` of the median step and raises ``ValueError`` if
+spacing is non-uniform.
+"""
+from __future__ import annotations
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from xrspatial.geotiff import _coords_to_transform, to_geotiff
+
+
+def _make_da(x_coords, y_coords):
+    arr = np.zeros((len(y_coords), len(x_coords)), dtype=np.uint8)
+    return xr.DataArray(
+        arr,
+        dims=['y', 'x'],
+        coords={'y': np.asarray(y_coords), 'x': np.asarray(x_coords)},
+    )
+
+
+def test_uniform_coords_ok():
+    """Uniform coords write successfully (no regression)."""
+    da = _make_da(np.linspace(500.0, 700.0, 20), np.linspace(100.0, 200.0, 10))
+    gt = _coords_to_transform(da)
+    assert gt is not None
+    np.testing.assert_allclose(gt.pixel_width, (700.0 - 500.0) / 19)
+    np.testing.assert_allclose(gt.pixel_height, (200.0 - 100.0) / 9)
+
+
+def test_uniform_coords_roundtrip_to_geotiff_1720(tmp_path):
+    """End-to-end write succeeds on uniform coords."""
+    da = _make_da(np.linspace(500.0, 700.0, 20), np.linspace(100.0, 200.0, 10))
+    p = str(tmp_path / 'uniform_1720.tif')
+    to_geotiff(da, p)
+
+
+def test_non_uniform_x_raises_1720():
+    """Non-uniform x coords raise ValueError naming x."""
+    # Mostly-uniform x with one stretched gap to expose the bug.
+    x = np.array([0.0, 1.0, 2.0, 3.0, 5.0])
+    y = np.linspace(0.0, 10.0, 11)
+    da = _make_da(x, y)
+    with pytest.raises(ValueError, match=r"\bx coords are not uniformly spaced"):
+        _coords_to_transform(da)
+
+
+def test_non_uniform_y_raises_1720():
+    """Non-uniform y coords raise ValueError naming y."""
+    x = np.linspace(0.0, 10.0, 11)
+    y = np.array([0.0, 1.0, 2.0, 3.0, 5.0])
+    da = _make_da(x, y)
+    with pytest.raises(ValueError, match=r"\by coords are not uniformly spaced"):
+        _coords_to_transform(da)
+
+
+def test_jitter_within_tolerance_ok_1720():
+    """Float jitter within 1e-6 relative tolerance writes successfully."""
+    x = np.linspace(0.0, 100.0, 11)
+    # Add jitter at ~1e-8 relative scale — well below the 1e-6 threshold.
+    rng = np.random.default_rng(0)
+    x = x + rng.uniform(-1e-7, 1e-7, size=x.shape)
+    y = np.linspace(0.0, 50.0, 6)
+    da = _make_da(x, y)
+    gt = _coords_to_transform(da)
+    assert gt is not None
+
+
+def test_jitter_just_above_tolerance_raises_1720():
+    """Jitter just above the 1e-6 relative tolerance raises ValueError."""
+    # Step size is 10; one diff is 10 + 1e-4 -> relative deviation ~1e-5,
+    # which exceeds the 1e-6 threshold.
+    x = np.array([0.0, 10.0, 20.0001, 30.0001, 40.0001])
+    y = np.linspace(0.0, 50.0, 6)
+    da = _make_da(x, y)
+    with pytest.raises(ValueError, match=r"max relative deviation"):
+        _coords_to_transform(da)
+
+
+def test_two_sample_coords_ok_1720():
+    """Two-sample coords (one diff) trivially pass the regularity check."""
+    x = np.array([0.0, 10.0])
+    y = np.array([0.0, 5.0])
+    da = _make_da(x, y)
+    gt = _coords_to_transform(da)
+    assert gt is not None
+    np.testing.assert_allclose(gt.pixel_width, 10.0)
+    np.testing.assert_allclose(gt.pixel_height, 5.0)
+
+
+def test_constant_coords_raises_1720():
+    """Constant coords (step == 0) raise ValueError."""
+    x = np.array([0.0, 0.0, 0.0, 0.0])
+    y = np.linspace(0.0, 10.0, 4)
+    da = _make_da(x, y)
+    with pytest.raises(ValueError, match=r"x coords are constant"):
+        _coords_to_transform(da)

xrspatial/geotiff/tests/test_edge_cases.py

@@ -512,8 +512,11 @@ def test_crs_projected(self, tmp_path):
 
     def test_pixel_is_point_round_trip(self, tmp_path):
         """PixelIsPoint raster_type preserves origin correctly."""
-        y = np.array([41.0, 40.999722, 40.999444, 40.999167])
-        x = np.array([-75.0, -74.999722, -74.999444, -74.999167])
+        # Use truly uniform coords; hand-typed fractional decimals had
+        # ~3.6e-3 relative deviation between steps and silently produced
+        # a wrong transform before issue #1720 enforced regularity.
+        y = np.linspace(41.0, 40.999167, 4)
+        x = np.linspace(-75.0, -74.999167, 4)
         da = xr.DataArray(
             np.arange(16, dtype=np.float32).reshape(4, 4),
             dims=['y', 'x'], coords={'y': y, 'x': x},