perf: stream the remaining eager Dataset methods (#969) - #971
Merged
Conversation
…d of always reading whole to_feature_collection read the whole raster into memory by default (tile=False). Make tile default to None and auto-select: read the whole array when it is at most 256 MiB, else read it tile by tile, so a huge or /vsicurl source is never materialised whole (#969). The tiled path already restores row-major order, so the result is byte-identical; an explicit tile=True/False still overrides the choice. Refs #969
…d copy sieve() seeded its in-memory target with dst_band.WriteArray(src_band.ReadAsArray()), materialising the whole source band as a NumPy array on top of the MEM buffer (peak memory ~2x the band). Replace the manual MEM Create + NumPy round trip with a gdal.Translate to a single-band MEM dataset: GDAL copies the band block by block in the C layer, carrying the geotransform, CRS, dtype, and no-data across, so no NumPy array is ever allocated. The gdal.SieveFilter call is unchanged; output is byte-identical. Refs #969
…sk whole _crop_aligned() read the full source (all bands) AND the full mask band into memory before stamping the mask's no-data layout onto the source (peak ~ source + mask). For a raster mask with no gap-filling, apply the mask tile by tile via a new _crop_aligned_tiled helper: each 256x256 window reads its own source and mask block, applies the no-data value, and writes straight into the destination, so the full arrays are never held together. The mask-apply is purely per-pixel, so the result is byte-identical to the whole-array path. The numpy-array mask (already in memory) and the fill_gaps path (interpolates across the whole array) keep the eager path. Refs #969
apply() read the whole band and applied func to every domain value at once. Add an opt-in elementwise=True that applies func one 256x256 tile at a time via a new _apply_elementwise_tiled helper, so a very large or /vsicurl source is never materialised whole. For a genuine per-pixel func the tiled result is byte-identical to the whole-array pass; the shared domain-apply logic is factored into _apply_func_to_domain so both paths behave identically. func that depend on the whole array (normalisation, ranking, any global reduction) must keep the default elementwise=False, whose behaviour is unchanged. Refs #969
plot_histogram() flattened the whole band and footprint() read the whole band to build the coverage mask, materialising a very large raster just to reduce it to a plot or a polygon. Add an opt-in max_samples cap to both: when set and the band has more cells than that, GDAL reads a nearest-neighbour decimated grid (~max_samples cells) in the C layer instead of the full band. A shared _read_decimated helper does the read; footprint scales the mask geotransform to the coarser grid via _scaled_geotransform so the polygon coordinates stay geographically correct. The result is then approximate (a subsample / coarser trace), documented as such. max_samples=None (default) reads every pixel, so both methods stay exact and byte-identical to before. Refs #969
change_no_data_value read each band whole (a full NumPy array per band) to swap the old no-data value for the new one. Keep GDAL's block-based CreateCopy (so every band's colour table, description, scale/offset, RAT and metadata survive exactly) but stream the swap one 256x256 tile at a time via a new _swap_no_data_tiled helper, so a full band is never held as a NumPy array. Add a path= option: when given, the clone is a disk-backed GeoTIFF and the tiled swap makes the whole operation out-of-core; FlushCache persists it before reopen. Output is byte-identical to before. Because NetCDF.read_array takes variable as its first positional, the old positional read_array(band) call incidentally raised a "pinned" ValueError on a variable-view NetCDF; the windowed read (band= keyword) removes that accident, so change_no_data_value now works on a variable view like fill() and the other band ops already do. The test that asserted the accidental guard is updated to assert it works. Refs #969
MAfarrag
changed the title
…ks on a NetCDF view apply()'s default path called read_array(band) positionally; because NetCDF.read_array puts variable first, that mis-bound band->variable and raised a spurious "pinned" error on a variable view, while apply(elementwise=True) (which passes band=) worked. Both paths now read with band= as a keyword, so apply behaves consistently on a variable view like the other band ops. The test that asserted the accidental guard now asserts both paths work and agree. Refs #969
… fails mid-swap If the tiled no-data swap raised mid-stream (e.g. a dtype-mismatch NoDataValueError) on the disk-backed path, a half-written GeoTIFF was left at path with the write handle still open. Wrap the swap in try/except that, on the disk path, closes the wrapper handle, drops the local dst reference (both are needed or GDAL keeps the file locked on Windows), and deletes the partial file plus its .aux.xml sidecar before re-raising. Refs #969
Add regression tests for the tiled/CreateCopy redesign: a per-band old_value list on a multi-band raster, the dtype guard firing on a band with zero matching cells (proving the swap is attempted per tile), RAT preservation through CreateCopy, inplace=True + path= updating the source and persisting to disk, and a mid-swap failure removing the partial disk file. Refs #969
…ts memory claim Sum the per-band itemsizes instead of assuming band 0's dtype, so a mixed-dtype stack is estimated correctly for the tile=None auto-select. Reword the comments/docstring: the tiled path avoids the full-band ndarray, but the resulting DataFrame still scales with the surviving non-no-data cells, so 'never materialised whole' overstated the bound. Refs #969
…ll raster The existing max_samples footprint tests used axis-aligned square-cell rasters, so the geotransform rotation terms (gt[2]/gt[4]) were always 0 and a regression dropping their scaling in _scaled_geotransform would pass. Add a fully-covered raster with non-zero rotation and non-square cells, asserting exact and decimated footprints share total_bounds. Refs #969
_crop_aligned_tiled called _apply_mask_nodata per tile, which re-ran the per-band _check_no_data_value dtype coercion on every window. Precompute the coerced list once before the loop and pass it in; _apply_mask_nodata still validates it itself when called without one (the non-tiled path), so behaviour is unchanged. Refs #969
When apply(elementwise=True) streamed a fully-no-data tile, domain_values was size 0; a func that raises on array input (a legitimate per-pixel map like lambda v: 1 if v > 5 else 0) fell back to np.vectorize(func), which raises "cannot call 'vectorize' on size 0 inputs unless otypes is set". The whole-array path never hit this (the band domain is non-empty), so the tiled path diverged. Return early on an empty domain -- out_array is already the no-data fill -- restoring byte-identity and also hardening the whole-array path for an all-no-data band. Refs #969
…ation max_samples is a public keyword on footprint and plot_histogram; max_samples=0 reached the decimation factor and raised ZeroDivisionError, and a negative value produced a complex factor and a 'complex has no __round__' TypeError. Reject max_samples < 1 at the read boundary with a clear ValueError. Refs #969
…t never masks the real error The disk-path failure cleanup unlinked the partial file with a bare unlink; if a GDAL/OS build still held the GeoTIFF locked, that unlink would raise PermissionError and replace the original NoDataValueError. Wrap each unlink in try/except OSError so a residual lock only lets the file linger and the original exception always propagates. Refs #969
…ta copy Tighten the loose decimation checks: the decimated footprint now must track the covered block's extent (not merely stay inside the raster) and its area within rel=0.3; the histogram test asserts the decimated samples are real band values (nearest-neighbour), not just a small count. Document that sieve's gdal.Translate seed also carries the band color table / RAT / scale-offset onto the result. Refs #969
Split the `assert X and Y` checks added in this PR into separate asserts so a failure points at the exact condition: the decimated-footprint extent bounds, the both-footprints existence check, the multi-band non-matching-cell preservation, and the variable-view apply both-paths check. Refs #969
…narCloud S3776) The added path=/cleanup logic pushed change_no_data_value's cognitive complexity to 22. Extract three helpers -- _normalize_no_data_arg (per-band arg validation), _swap_all_bands (the per-band swap loop), and _discard_partial_output (the disk-path failure cleanup) -- so the method reads as a linear clone -> set -> swap -> flush with a single try/except. Behaviour is unchanged. Refs #969
test_peak_memory_is_bounded_by_the_tile asserted a hard-coded Python-heap budget (dense_bytes // 4) that does not hold across GDAL builds: on the macOS cibuildwheel GDAL the streamed pass traced ~3.9 MB (vs the 0.5 MB budget) even though it tiles correctly, because tracemalloc only sees the Python heap, not GDAL's C buffers, and the wheel build's read/write path allocates differently. Replace the absolute budget with a relative one: measure a whole-array pass (same read -> transform -> write, materialised at once) in the same process and assert the streamed peak stays below it. The baseline shares the streamed path's machinery, so any per-build overhead cancels and only the full-array data streaming avoids remains -- locally the streamed peak is ~2% of the whole-array peak.
|
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
1 participant
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.
Description
Optimizes the remaining eager, full-array
Datasetmethods from the #969 audit so a very large or/vsicurl-backed raster is no longer materialised whole. Each method keeps its exact output (byte-identical)unless a new opt-in flag is passed; the streaming technique is matched to the shape of each method.
to_feature_collectiontiledefaults toNone(auto): read whole below 256 MiB, else tile by tilesievegdal.Translate(block-based C copy) instead ofWriteArray(ReadAsArray())crop(mask=<raster>)_crop_aligned(no-gap-fill case)applyelementwise=Truestreamsfunctile by tilefunc; default unchangedplot_histogrammax_samples=reads a decimated gridfootprintmax_samples=builds the coverage mask from a decimated read (geotransform scaled)change_no_data_valueCreateCopyfor metadata fidelity); addpath=for a disk-backed, out-of-core resultDesign notes:
apply,plot_histogram, andfootprintonly change behaviour when their newopt-in (
elementwise=True/max_samples=) is passed; otherwise the output is identical to before. The otherfour are byte-identical unconditionally.
change_no_data_valuekeeps GDAL'sCreateCopydeliberately — it clones every band's colour table,description, scale/offset, RAT, and metadata exactly, which an
empty_like+ manual metadata copy would riskdropping. The memory win comes from streaming the no-data swap one tile at a time; the new
path=makes theclone a disk-backed GeoTIFF so the whole operation is genuinely out-of-core.
change_no_data_valuecalledread_array(band)positionally, and since
NetCDF.read_arraytakesvariablefirst, that incidentally raised a "pinned"ValueErroron a variable view. The windowed read (band=keyword) removes that accident, sochange_no_data_valuenow works on a variable view likefill()and the other band ops already do. Thetest that asserted the accidental guard is updated to assert it works.
Issues
Type of change
Check relevant points.
Notes: default output is byte-identical for every method; the new parameters (
tile=None,elementwise=,max_samples=,path=) are additive and backward-compatible. The one behaviour change is thatchange_no_data_valuenow succeeds on a NetCDF variable view instead of raising an (accidental) error.How Has This Been Tested?
sieve: existing 13test_sieve.pycases (nodata both ways, geotransform/CRS, band selection, mask) pinthe byte-identical output through the new
gdal.Translateseeding.crop(mask=raster):TestCropAlignedTiling(new) proves the tiled apply matches a direct NumPy applyacross 256-px tile seams (1- and 2-band) and matches the eager numpy-mask path; the existing
TestCropcasesstill pass.
apply:test_apply_elementwise.py(new) — streamed vs whole-array equal across tile boundaries forabs/affine/sqrt, no-data preserved, band selection, in-place.plot_histogram/footprint: new tests assert the exact default is unchanged and thatmax_samples=decimates (fewer samples / coarser polygon within bounds).
change_no_data_value:TestChangeNoDataValueStreaming(new) — multi-tile swap matches a direct swap,colour table + scale/offset + description survive, and
path=persists a disk-backed result; the NetCDFvariable-view test now asserts it works.
tests/dataset/{analysis,spatial,unit,create}+ the touched NetCDF tests: 1288passed.
mypyclean (136 source files); doctests green foranalysis.py,spatial.py,bands.py,vectorize.py.Reproduce:
pixi run -e dev pytest tests/dataset/analysis tests/dataset/spatial tests/dataset/unit tests/dataset/create -m "not plot" -qChecklist: