diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 70f35984..a229830e 100644 --- a/xrspatial/geotiff/__init__.py +++ b/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]) diff --git a/xrspatial/geotiff/tests/test_coord_regularity_1720.py b/xrspatial/geotiff/tests/test_coord_regularity_1720.py new file mode 100644 index 00000000..615fbc64 --- /dev/null +++ b/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) diff --git a/xrspatial/geotiff/tests/test_edge_cases.py b/xrspatial/geotiff/tests/test_edge_cases.py index edbdab56..01ab163e 100644 --- a/xrspatial/geotiff/tests/test_edge_cases.py +++ b/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},