Skip to content

geotiff: preserve transform on 1xN / Nx1 georeferenced writes (#1945)#1953

Merged
brendancol merged 2 commits into
mainfrom
issue-1945
May 15, 2026
Merged

geotiff: preserve transform on 1xN / Nx1 georeferenced writes (#1945)#1953
brendancol merged 2 commits into
mainfrom
issue-1945

Conversation

@brendancol
Copy link
Copy Markdown
Contributor

Summary

Closes #1945.

  • coords_to_transform now recovers the per-axis pixel size from the non-degenerate axis on 1xN / Nx1 inputs and assumes square pixels for the degenerate one, instead of bailing out and letting the writer produce a non-georeferenced TIFF.
  • For 1x1 inputs (no axis carries a step), attrs['transform'] is honoured; if none is supplied, the writer fails closed with a ValueError rather than silently dropping spatial metadata.
  • A new require_transform_for_georeferenced helper enforces the fail-closed rule from each writer entry point: to_geotiff eager, to_geotiff dask, to_geotiff with a .vrt path (VRT-tiled), and write_geotiff_gpu.

Why this is a behaviour change

Before:

da = xr.DataArray(
    np.arange(8, dtype='float32').reshape(1, 8),
    dims=('y', 'x'),
    coords={'x': np.linspace(-120.0, -113.0, 8), 'y': np.array([45.0])},
    attrs={'crs': 4326},
)
to_geotiff(da, '/tmp/strip.tif')
open_geotiff('/tmp/strip.tif').coords['x'].values[:3]  # array([0, 1, 2]) int64

After: the coords round-trip back as floats, and attrs['transform'] is populated.

A 1x1 DataArray without an explicit attrs['transform'] now raises ValueError instead of silently producing a non-georeferenced file.

Test plan

  • New tests in xrspatial/geotiff/tests/test_degenerate_georef_1945.py: 1xN, Nx1, and 1x1 round-trips across numpy, dask+numpy, GPU (cupy), dask+cupy, and the VRT-tiled writer (15 tests).
  • Full geotiff test suite: 3109 passed locally. The 8 failures on this branch already fail on main (unrelated read_to_array monkeypatch and signature-parity tests).

Copilot AI review requested due to automatic review settings May 15, 2026 15:29
@github-actions github-actions Bot added the performance PR touches performance-sensitive code label May 15, 2026
Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot encountered an error and was unable to review this pull request. You can try again by re-requesting a review.

@brendancol
Copy link
Copy Markdown
Contributor Author

PR Review: geotiff: preserve transform on 1xN / Nx1 georeferenced writes (#1945)

Blockers (must fix before merge)

  • None

Suggestions (should fix, not blocking)

  • xrspatial/geotiff/_coords.py:264require_transform_for_georeferenced error message says "at least one axis has length 1, so the pixel size cannot be recovered from coord differences", but the new coords_to_transform does recover from 1xN / Nx1 just fine; this helper only fires for the 1x1 case (both axes len < 2). Wording could read "both axes have length 1" to match the actual trigger condition and avoid misleading users who hit it.
  • xrspatial/geotiff/tests/test_degenerate_georef_1945.py:75 — only the default raster_type='area' is exercised. A raster_type='point' 1xN / Nx1 round-trip would lock in the origin offset path (no half-pixel shift) alongside the area path the tests already cover.
  • xrspatial/geotiff/_coords.py:367 — borrow comments say "x increases left-to-right by convention" and "y decreases top-to-bottom by convention" but the borrow uses abs(...), which discards the sign carried by the non-degenerate axis. For a Nx1 raster stored bottom-up (positive pixel_height), the synthesised pixel_width is still positive, which happens to match the convention but is the wrong answer if the caller is mirroring an unusual source raster. A test pinning the sign behaviour on a y-ascending Nx1 input would document the choice.

Nits (optional improvements)

  • xrspatial/geotiff/_coords.py:267require_transform_for_georeferenced calls da.dims[-2] for any non-3D input. A 1D DataArray would raise an IndexError before the spatial-coord check. The writers already reject 1D earlier, but a defensive if da.ndim < 2: return would make the helper self-contained.
  • xrspatial/geotiff/_coords.py:352 — the inline comment uses an em dash; project style elsewhere in the file uses hyphens or parentheses for asides.

What looks good

  • 1x1 fail-closed wired into all four writer entry points (eager, dask via to_geotiff, VRT, GPU) through a single helper, so the policy stays in sync.
  • The flip from or to and on len(x) < 2 keeps the 1x1 short-circuit while letting the 1xN / Nx1 fallback run.
  • Round-trip tests cover numpy, dask+numpy, GPU, dask+cupy, and VRT paths plus the 1x1 raise.
  • attrs['transform'] precedence is preserved so the 1x1 case has a clean opt-in.

Checklist

  • Algorithm matches reference/paper
  • All implemented backends produce consistent results
  • NaN handling is correct
  • Edge cases are covered by tests
  • Dask chunk boundaries handled correctly
  • No premature materialization or unnecessary copies
  • Benchmark exists or is not needed
  • README feature matrix updated (if applicable)
  • Docstrings present and accurate

@brendancol
Copy link
Copy Markdown
Contributor Author

PR Review: geotiff: preserve transform on 1xN / Nx1 georeferenced writes (#1945) (re-review after rebase)

Re-reviewed at head 81d9704 after a clean rebase onto origin/main (which now carries the Windows fix 9ce0e60). Local run of test_degenerate_georef_1945.py is 19 passed in 0.58s on the CPU paths; GPU classes are skipped here. Force-pushed with lease.

Blockers

  • None

Suggestions

  • None

Nits

  • The error message in require_transform_for_georeferenced calls out the 1x1 case in detail. For an unusual future shape (say a 0-length axis sneaking through) the wording would still read 1x1. Not worth changing for this PR, but a phrasing like "both spatial axes are degenerate (len < 2)" would future-proof it.
  • _strip_1xN/_strip_Nx1 accept a raster_type argument but mutate the attrs dict in place. Building the attrs dict per call (as the code does) is fine, just noting that callers reusing the literal would share state if that pattern ever changed.

What looks good

  • The split is right: coords_to_transform now returns a usable transform whenever at least one axis carries a step, and require_transform_for_georeferenced is the single fail-closed gate that catches the 1x1 case before a non-georeferenced TIFF gets written. The guard is wired into all three writer entry points (to_geotiff, _write_vrt_tiled, write_geotiff_gpu) with matching comments, so the behaviour cannot drift between CPU, dask, GPU, and VRT paths.
  • The borrow convention (pixel_width = abs(pixel_height), pixel_height = -abs(pixel_width)) is pinned by TestCoordsToTransformBorrowSignPinning against both an ascending-y Nx1 and a descending-x 1xN. That nails down the sign choice so a refactor cannot quietly flip it.
  • Test matrix covers numpy, dask+numpy, cupy, dask+cupy, and VRT, plus a raster_type='point' parametrisation that pins the no-half-pixel-shift path on degenerate writes. The round-trip assertion checks both float coord dtype and that attrs['transform'] came back, which is exactly the regression geotiff: 1xN / Nx1 georeferenced writes silently strip transform #1945 reported.
  • Reader-side behaviour is untouched; this is a writer-only change with a clear error path for the genuinely ambiguous shape.

Checklist

  • Rebase clean on origin/main (no conflicts)
  • test_degenerate_georef_1945.py passes locally on CPU (19/19, GPU skipped)
  • Fail-closed guard reaches CPU, dask, GPU, and VRT writers
  • Sign convention on borrowed axis covered by a dedicated test
  • raster_type='point' covered for both 1xN and Nx1
  • CI green on Windows 3.12 after the rebase (depends on 9ce0e60; verify on GH)

coords_to_transform returned None whenever either spatial coord had
length < 2, and the writers fell through and produced a non-georeferenced
TIFF. A single-row or single-column DataArray with real x/y coords
round-tripped back with integer pixel coords and no transform attribute,
silently dropping spatial metadata.

Recover the per-axis pixel size from the non-degenerate axis when only
one axis is degenerate (1xN / Nx1) and assume square pixels for the
degenerate one. For the 1x1 case where no axis carries a step, prefer
attrs['transform']; if no explicit transform is supplied, fail closed
with a clear ValueError rather than writing a non-georeferenced file.

Wire the same fail-closed guard into the eager, dask, GPU, and VRT-tiled
writer paths via a new require_transform_for_georeferenced helper.

Add round-trip tests covering 1xN, Nx1, and 1x1 across numpy, dask+numpy,
GPU (cupy), dask+cupy, and the VRT-tiled writer.
@brendancol
Copy link
Copy Markdown
Contributor Author

Rebased onto current main at 2b72bba to pick up the Python 3.14 CI matrix limit (b1579f8).

@brendancol brendancol merged commit 9a8687a into main May 15, 2026
5 checks passed
brendancol added a commit that referenced this pull request May 15, 2026
…el parity with coords_to_transform) (#1969)

PR #1953 added require_transform_for_georeferenced which raises when
both spatial dims are in da.coords and the resolved transform is None.
PR #1954 made coords_to_transform return None for integer x/y coords
as the no-georef sentinel (#1949). Their interaction broke writer
calls against int-coord arrays. Mirror the int/uint exemption from
coords_to_transform inside the guard so the writer accepts int coords
silently. Also fixes the error wording, which previously claimed
"both axes are degenerate (1x1)" even when the array was not 1x1.

Adds a regression test covering 2D and 3D int-coord round-trips.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

performance PR touches performance-sensitive code

Projects

None yet

Development

Successfully merging this pull request may close these issues.

geotiff: 1xN / Nx1 georeferenced writes silently strip transform

2 participants