Fix geodesic coord resolution to prefer real lat/lon over numeric y/x#2903
Merged
Conversation
…#2896) _find_coord accepted the numeric dimension coord before scanning for a coord named latitude/longitude, so a curvilinear raster with numeric y/x index coords plus real 2-D lat/lon coords used the pixel indices as lat/lon. The range check did not catch it because small indices fall inside the geographic ranges. Prefer an explicitly named lat/lon coord over the dimension coord. This covers slope, aspect, and surface_distance, which share _extract_latlon_coords. Add curvilinear coverage (dims y/x + 2-D lat/lon + numeric y/x) across the four backends.
brendancol
commented
Jun 3, 2026
Contributor
Author
brendancol
left a comment
There was a problem hiding this comment.
PR Review: Fix geodesic coord resolution to prefer real lat/lon over numeric y/x
Blockers (must fix before merge)
None.
Suggestions (should fix, not blocking)
None.
Nits (optional improvements)
- xrspatial/utils.py:967-972: The explicit-name scan picks the first matching coord in
agg.coordsiteration order. If a DataArray carries both alatand alatitudecoord (orlon/longitude) with different dimensionality, the winner depends on coord insertion order, and a 1-D winner on one axis paired with a 2-D winner on the other would trip the "both 1-D or both 2-D" error in_extract_latlon_coords. This is an unusual input and the resulting error is clear, so it does not need handling now, but a one-line docstring note that duplicate explicit names resolve by coord order would help a future reader.
What looks good
- Root cause fixed at the right layer: reordering
_find_coordprecedence fixes slope, aspect, and surface_distance at once, since they all route through_extract_latlon_coords. - The
y/xfallback is preserved, so rasters that legitimately use y/x as geographic coords still resolve (confirmed by the still-passingtest_projected_coords_raises). - Tests assert against the 1-D lat/lon reference rather than just "runs without error", and the
> 0.05interior check would catch a regression back to the pixel-index behavior (~0.007). - Backend parity covered across numpy, cupy, dask+numpy, dask+cupy via the general_checks markers.
Checklist
- Algorithm matches reference: n/a (coord resolution, not a numeric kernel)
- All implemented backends produce consistent results: yes, parity tests added
- NaN handling correct: unchanged by this PR
- Edge cases covered by tests: curvilinear path covered; duplicate-explicit-name case noted as a nit
- Dask chunk boundaries handled correctly: unchanged; coord resolution is host-side
- No premature materialization: coord
.valuesread is host-side and pre-existing - Benchmark exists or not needed: not needed (no perf-sensitive path changed)
- README feature matrix updated: not applicable (bug fix, no new API)
- Docstrings present and accurate: yes,
_find_coorddocstring updated
brendancol
commented
Jun 3, 2026
Contributor
Author
brendancol
left a comment
There was a problem hiding this comment.
PR Review (follow-up pass)
Re-reviewed after commit 40a0065.
Blockers
None.
Suggestions
None.
Nits
None outstanding. The earlier nit on _find_coord (utils.py) is resolved: the docstring now states that when several explicit lat/lon names are present, the first one in coord order wins.
What looks good
- The fix and tests are unchanged from the first pass and still green: 103 passed across geodesic slope, geodesic aspect, and surface_distance suites locally (cupy / dask+cupy classes skip without a GPU).
- The diff stays scoped to coord resolution plus its tests. No unrelated edits.
Nothing further to address. Recommending merge once CI is green.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #2896
Problem
_find_coord()inxrspatial/utils.pyaccepted the numeric dimension coordinate before it scanned for a coordinate named latitude/longitude. For a curvilinear raster with dims ('y', 'x'), numeric y/x pixel-index coords, and real 2-D lat/lon coords, the geodesic path used the pixel indices as lat/lon. The range check did not catch it because small indices (0..N) fall inside the accepted geographic ranges. The result was a wrong slope (~0.007 vs ~0.067 on the reproduction grid) with no warning.Change
_find_coord()now prefers an explicitly named lat/lon coord (lat/latitude/lon/longitude) over the dimension coord. The dimension coord and the broadery/xname set remain as fallbacks, so rasters with y/x as the geographic coords still work._extract_latlon_coords/_find_coord. The logic is coord resolution only, so it is backend-agnostic.Backend coverage
numpy, cupy, dask+numpy, dask+cupy. Coord resolution runs on the host before dispatch; the new curvilinear tests assert parity across all four (cupy / dask+cupy skip without a GPU).
Test plan
TestGeodesicSlopeCurvilinear*classes: dims y/x + 2-D lat/lon + numeric y/x, asserting correct slope and cross-backend paritytest_projected_coords_raises, which relies on the y/x fallback)