zonal: fix int32 overflow in _strides() that silently corrupts stats/crosstab on huge rasters#2617
Merged
Merged
Conversation
…rflow (#2612) _strides accumulated cumulative zone-boundary counts into an int32 array inside a numba @ngjit kernel. For rasters above ~2.1 billion pixels the count wraps silently to a negative value, corrupting the slice bounds that _calc_stats and crosstab derive from it, so stats() and crosstab() returned silently wrong results on very large numpy-backed (or large-chunk dask) rasters. Allocate the counter as int64 to match the natural array-index type.
brendancol
commented
May 29, 2026
brendancol
left a comment
Contributor
Author
There was a problem hiding this comment.
PR Review: fix int32 overflow in _strides()
Blockers
None.
Suggestions
None.
Nits
- The test pins the dtype and checks correctness on a tiny input. It can't actually trigger the wrap without a ~17 GB array, so the dtype assertion is doing the real work here. A one-line comment to that effect would save the next reader a moment. Not blocking.
What looks good
- Root cause is right. In a numba nopython kernel, storing a count past 2**31-1 into an int32 array wraps silently instead of raising, so int64 is the correct fix and it matches how the value is later used (as a slice index).
- Minimal change, no surrounding churn.
_stridesis the one chokepoint feedingzone_breaksinto_calc_statsand the crosstab paths, so this coversstats()andcrosstab()on numpy, dask+numpy, and the cupy/dask+cupy paths (which route through the numpy code). Nothing else to touch per-backend.- Downstream uses these as plain Python slice bounds, and the index arrays come from
np.argsort(already int64 on 64-bit), so there's no other int32 narrowing waiting to reintroduce this. - Existing stats/crosstab tests still pass (111), so the wider dtype didn't move results at normal sizes.
Checklist
- Algorithm vs reference: n/a (overflow fix)
- Backends consistent: yes, one shared codepath
- NaN handling: unchanged
- Edge cases: dtype + correctness covered; a real overflow isn't allocatable
- Dask boundaries: unchanged, per-block strides still correct
- Premature materialization: none
- Benchmark: not needed
- README matrix: n/a
- Docstrings: accurate
brendancol
commented
May 29, 2026
brendancol
left a comment
Contributor
Author
There was a problem hiding this comment.
Follow-up: addressed the one nit from the previous review. The test now carries a comment explaining that a real >2**31-element input isn't allocatable (~17 GB), so the int64 dtype assertion stands in for the overflow case. No other findings were open. Regression test still passes.
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 #2612
What changed
_strides()accumulated cumulative zone-boundary counts into annp.int32array inside a numba@ngjitkernel. For rasters above ~2.1 billion pixels (about 46341 x 46341) the count wraps silently to a negative int32, so the slice bounds_calc_statsand the crosstab paths derive from it go negative andstats()/crosstab()return silently wrong results.int64, matching the natural array-index type, so it cannot overflow at realistic raster sizes._stridesreturns anint64, non-negative, monotonic result.Backend coverage
_stridesruns in the numpy and dask+numpy paths (and the cupy/dask+cupy paths, which route through the numpy implementation). The dtype fix applies to all of them.Test plan
test_strides_int64_no_overflow_2612passes