GeoTIFF reader/writer and reproject with GPU + datum transforms#1062
Merged
brendancol merged 48 commits intomasterfrom Mar 23, 2026
Merged
GeoTIFF reader/writer and reproject with GPU + datum transforms#1062brendancol merged 48 commits intomasterfrom
brendancol merged 48 commits intomasterfrom
Conversation
Reads and writes GeoTIFF and Cloud Optimized GeoTIFF files using only numpy, numba, xarray, and the standard library. No GDAL required. What it does: - Parses TIFF/BigTIFF headers and IFDs via struct - Reads GeoTIFF tags: CRS (EPSG), affine transforms, GeoKeys, PixelIsArea vs PixelIsPoint - Deflate (zlib), LZW (Numba JIT), horizontal and floating-point predictor codecs - Strip and tiled layouts, windowed reads - COG writer with overview generation and IFD-first layout - HTTP range-request reader for remote COGs - mmap for local file access (zero-copy) - Nodata sentinel masking to NaN on read - Metadata round-trip: CRS, transform, nodata, raster type, pixels Read performance vs rioxarray/GDAL: - 5-7x faster on uncompressed data - On par for compressed COGs - 100% pixel-exact across all tested formats Write performance: - On par for uncompressed - 2-4x slower for compressed (zlib is C-native in GDAL) Tested against Landsat 8, Copernicus DEM, USGS 1-arc-second, and USGS 1-meter DEMs. 154 tests cover codecs, header parsing, geo metadata, round-trips across 8 dtypes x 3 compressions, edge cases (corrupt files, NaN/Inf, extreme shapes, PixelIsPoint), and the public API.
…rite Six features filling the main gaps for real-world use: 1. Multi-band write: 3D arrays (height, width, bands) now write as multi-band GeoTIFFs with correct BitsPerSample, SampleFormat, and PhotometricInterpretation (RGB for 3+ bands). Overviews work for multi-band too. read_geotiff returns all bands by default (band=None) with a 'band' dimension. 2. Integer nodata masking: uint8/uint16/int16 arrays with nodata values are promoted to float64 and masked with NaN on read, matching rioxarray behavior. Previously only float arrays were masked. 3. PackBits compression (tag 32773): simple RLE codec, both read and write. Common in older TIFF files. 4. JPEG decompression (tag 7): read support via Pillow for JPEG-compressed tiles/strips. Import is optional and lazy. 5. BigTIFF write: auto-detects when output exceeds ~4GB and switches to BigTIFF format (16-byte header, 20-byte IFD entries, 8-byte offsets). Prevents silent offset overflow corruption on large files. 6. Dask lazy reads: read_geotiff_dask() returns a dask-backed DataArray using windowed reads per chunk. Works for single-band and multi-band files with nodata masking per chunk. 178 tests passing.
Strip-based windowed reads now only decompress strips that overlap the requested row range. Previously, all strips were decompressed into a full image buffer and then sliced. For a 4096x512 deflate file with 256-row strips, reading a 10x10 window from the top-left goes from 31 ms to 1.9 ms (16x). On a 100,000-row file the savings scale linearly with the number of strips skipped.
Read and write Zstandard-compressed GeoTIFFs using the zstandard package (lazy import, clear error if missing). On a 2048x2048 float32 raster, ZSTD vs deflate: - Write: 39 ms vs 420 ms (10.7x faster) - Read: 14 ms vs 66 ms (4.7x faster) - Size: 15.5 MB vs 15.5 MB (comparable) 9 new tests covering codec round-trips, stripped/tiled layouts, uint16, predictor, multi-band, and the public API.
TIFF PlanarConfiguration=2 stores each band as a separate set of strips or tiles (RRR...GGG...BBB) instead of interleaved (RGBRGB...). The reader now detects this from the IFD and iterates band-by-band through the strip/tile offset array, placing each single-band chunk into the correct slice of the output. Both strip and tile layouts are handled. Windowed reads and single- band selection work correctly with planar files. 6 new tests: planar strips (RGB, 2-band), planar tiles, windowed read, band selection, and public API.
Adds read support for non-byte-aligned pixel data, common in bilevel masks (1-bit), palette images (4-bit), and medical/scientific sensors (12-bit). Changes: - _dtypes.py: Map 1/2/4-bit to uint8 and 12-bit to uint16 - _compression.py: Add unpack_bits() for MSB-first bit unpacking at 1, 2, 4, and 12 bits per sample - _reader.py: Add _decode_strip_or_tile() helper that handles the full decompress -> predictor -> unpack -> reshape pipeline, detecting sub-byte depths automatically. Both strip and tile readers refactored to use it. 7 new tests: 1-bit bilevel, 1-bit non-byte-aligned width, 4-bit, 4-bit odd width, 12-bit, and direct codec tests.
Reads TIFF files with Photometric=3 (Palette) and ColorMap tag (320).
The TIFF color table (uint16 R/G/B arrays) is converted to a
matplotlib ListedColormap and stored in da.attrs['cmap'].
New plot_geotiff() convenience function uses the embedded colormap
with BoundaryNorm so that integer class indices map to the correct
palette colors when plotted. Works out of the box:
da = read_geotiff('landcover.tif')
plot_geotiff(da) # colors match the TIFF's palette
Also stores raw RGBA tuples in attrs['colormap_rgba'] for custom use.
Supports both 8-bit (256-color) and 4-bit (16-color) palettes.
5 new tests: 8-bit palette read, 4-bit palette, colormap object
verification, plot_geotiff smoke test, and non-palette attr check.
The .xrs accessor (registered on all DataArrays by xrspatial) now has
a plot() method that checks for an embedded TIFF colormap in attrs.
If present, it applies BoundaryNorm with the ListedColormap so that
integer class indices map to the correct palette colors.
da = read_geotiff('landcover.tif')
da.xrs.plot() # palette colors applied automatically
For non-palette DataArrays, falls through to the standard da.plot().
The old plot_geotiff() function is kept as a thin wrapper.
Multiple threads reading the same file now share a single read-only mmap instead of each opening their own. A module-level _MmapCache protected by a threading.Lock manages reference counts per file path. The mmap is closed when the last reader releases it. Read-only mmap slicing (which is what the strip/tile readers do) is thread-safe at the OS level -- no seek or file position involved. Tested with 16 concurrent threads reading different windows from the same deflate+tiled file, and a stress test of 400 reads across 8 threads. Zero errors, cache drains properly. For dask lazy reads, this means all chunk-read tasks for the same file share one mmap instead of opening/closing the file per chunk.
Writes now go to a temporary file in the same directory, then os.replace() atomically swaps it over the target path. This gives: - No interleaved output when multiple threads write the same path - Readers never see a half-written file - No corrupt file left behind if the process crashes mid-write - Temp file cleaned up on any exception os.replace is atomic on POSIX (single rename syscall) and near-atomic on Windows (ReplaceFile).
_make_overview() now accepts a method parameter instead of hardcoding 2x2 block averaging. Available methods: - mean (default): nanmean of each 2x2 block, same as before - nearest: top-left pixel of each block (no interpolation) - min/max: nanmin/nanmax of each block - median: nanmedian of each block - mode: most frequent value per block (for classified rasters) - cubic: scipy.ndimage.zoom with order=3 (requires scipy) All methods work on both 2D and 3D (multi-band) arrays. Exposed via overview_resampling= parameter on write() and write_geotiff(). 12 new tests covering each method, NaN handling, multi-band, COG round-trips with nearest and mode, the public API, and error on invalid method names.
Adds support for TIFF resolution metadata used in print and cartographic workflows: - Tag 282 (XResolution): pixels per unit, stored as RATIONAL - Tag 283 (YResolution): pixels per unit, stored as RATIONAL - Tag 296 (ResolutionUnit): 1=none, 2=inch, 3=centimeter Read: resolution values are stored in DataArray attrs as x_resolution, y_resolution (float), and resolution_unit (string: 'none', 'inch', or 'centimeter'). Write: accepted as keyword args on write() and write_geotiff(), or extracted automatically from DataArray attrs. Written as RATIONAL tags (numerator/denominator pairs). Also adds RATIONAL type serialization to the writer's tag encoder. 5 new tests: DPI round-trip, centimeter unit, no-resolution check, DataArray attrs preservation, unit='none'.
…ical CRS GeoInfo and DataArray attrs now include all commonly-used GeoKeys parsed from the GeoKeyDirectory: - crs_name: full CRS name (e.g. "NAD83 / UTM zone 18N") - geog_citation: geographic CRS name (e.g. "WGS 84", "NAD83") - datum_code: EPSG geodetic datum code - angular_units / angular_units_code: e.g. "degree" (9102) - linear_units / linear_units_code: e.g. "metre" (9001) - semi_major_axis, inv_flattening: ellipsoid parameters - projection_code: EPSG projection method code - vertical_crs, vertical_citation, vertical_units: for compound CRS EPSG unit codes are resolved to human-readable names via lookup tables (ANGULAR_UNITS, LINEAR_UNITS). Raw geokeys dict is still available for anything not covered by a named field. 6 new tests covering geographic and projected CRS extraction, real-file verification, no-CRS baseline, and unit lookups.
_HTTPSource now uses a module-level urllib3.PoolManager that reuses TCP connections (and TLS sessions) across range requests to the same host. For a COG with 64 tiles, this eliminates 63 TCP handshakes. On localhost: 1.7x faster for 16 range requests. Over real networks the benefit is much larger since each TLS handshake adds 50-200ms. Falls back to stdlib urllib.request if urllib3 is not installed. The pool is created lazily on first use with retry support (2 retries, 0.1s backoff).
CRS can now be specified as WKT strings, PROJ strings, or EPSG integers. pyproj (lazy import) resolves between them: Read side: - crs_wkt attr is populated by resolving EPSG -> WKT via pyproj - Falls back gracefully if pyproj is not installed (EPSG still works) Write side: - crs= parameter on write_geotiff accepts int (EPSG), WKT string, or PROJ string. String inputs are resolved to EPSG via pyproj.CRS.from_user_input().to_epsg(). - DataArray with crs_wkt attr (no integer crs) is also handled: the WKT is resolved to EPSG for the GeoKeyDirectory. This means files with user-defined CRS no longer lose their spatial reference when round-tripped, as long as pyproj can resolve the WKT/PROJ string to an EPSG code. 5 new tests: WKT from EPSG, write with WKT string, write with PROJ string, crs_wkt attr round-trip, and no-CRS baseline.
Band names, statistics, and other GDAL-specific metadata stored in
the GDALMetadata XML tag (42112) are now read, exposed, and written
back.
Read: the XML is parsed into a dict stored in attrs['gdal_metadata'].
Dataset-level items use string keys ('DataType'), per-band items use
(name, band_int) tuple keys (('STATISTICS_MAXIMUM', 0)). The raw XML
is also available in attrs['gdal_metadata_xml'].
Write: accepts gdal_metadata_xml on write(), or extracts from
DataArray attrs on write_geotiff(). If attrs has a gdal_metadata
dict but no raw XML, it's re-serialized automatically.
Round-trip verified on the USGS 1-meter DEM which has statistics,
layer type, and data type metadata -- all items survive intact.
7 new tests: XML parsing, dict serialization, file round-trip,
DataArray attrs preservation, no-metadata baseline, real-file read,
and real-file round-trip.
Any IFD tag that the writer doesn't explicitly manage (Software, DateTime, ImageDescription, Copyright, custom private tags, etc.) is now collected on read, stored in attrs['extra_tags'], and re-emitted on write. Read: extract_geo_info collects (tag_id, type_id, count, value) tuples for all tags not in the _MANAGED_TAGS set (structural tags that the writer builds from scratch: dimensions, compression, offsets, geo tags, etc.). Stored in attrs['extra_tags']. Write: extra_tags are appended to the IFD, skipping any tag_id that was already written to avoid duplicates. The tag values are serialized using the same type-aware encoder as built-in tags. Tested with a hand-crafted TIFF containing Software (305) and DateTime (306) tags. Both survive read -> write -> read intact. 3 new tests: read detection, round-trip preservation, and no-extra-tags baseline.
The auto-detection now estimates total file size (header + IFDs + overflow + pixel data) instead of only checking compressed pixel data size, and compares against UINT32_MAX (4,294,967,295) instead of a hardcoded 3.9 GB threshold. Also adds a bigtiff= parameter to write() and write_geotiff(): - bigtiff=None (default): auto-detect based on estimated file size - bigtiff=True: force BigTIFF even for small files - bigtiff=False: force classic TIFF (user's responsibility if >4GB) 3 new tests: force BigTIFF via public API, small file stays classic, force classic via bigtiff=False.
Big-endian TIFFs (byte order marker 'MM') now byte-swap pixel data to native order after decompression. Previously, the reader did .view(dtype) with a native-order dtype, producing garbage values for multi-byte types (uint16, int32, float32, float64). Fix: _decode_strip_or_tile uses dtype.newbyteorder(file_byte_order) for the view, then .astype(native_dtype) if a swap is needed. Single-byte types (uint8) need no swap. The COG HTTP reader path has the same fix. Also fixed the test conftest: make_minimal_tiff(big_endian=True) now actually writes pixel bytes in big-endian order. 7 new tests: float32, uint16, int32, float64, uint8 (no swap), windowed read, and public API -- all with big-endian TIFFs.
Read and write GeoTIFFs to/from cloud storage using fsspec as the filesystem abstraction layer. Any URI with a :// scheme (that isn't http/https) is routed through fsspec, which delegates to the appropriate backend: - s3://bucket/key.tif (requires s3fs) - gs://bucket/key.tif (requires gcsfs) - az://container/blob.tif (requires adlfs) - abfs://container/blob.tif (requires adlfs) - Any other fsspec-supported scheme (memory://, ftp://, etc.) Read: _CloudSource uses fsspec.core.url_to_fs() then fs.open() for full reads and range reads. Falls through to the same TIFF parsing pipeline as local files. Write: _write_bytes detects fsspec URIs and writes via fs.open() instead of the local atomic-rename path (which doesn't apply to cloud storage). If fsspec or the backend library isn't installed, a clear ImportError is raised with install instructions. 5 new tests using fsspec's memory:// filesystem for integration testing without real cloud credentials.
Reads GDAL .vrt files by parsing the XML and assembling pixel data
from the referenced source GeoTIFFs using windowed reads.
Supported VRT features:
- SimpleSource: direct pixel copy with source/destination rects
- ComplexSource: scaling (ScaleRatio) and offset (ScaleOffset)
- Source nodata masking
- Multiple bands
- GeoTransform and SRS/CRS propagation
- Relative and absolute source file paths
- Windowed reads (only fetches overlapping source regions)
Usage:
da = read_geotiff('mosaic.vrt') # auto-detected by extension
da = read_vrt('mosaic.vrt') # explicit function
da = read_vrt('mosaic.vrt', window=(0, 0, 100, 100)) # windowed
read_geotiff auto-detects .vrt files and routes them through the
VRT reader. The DataArray gets coordinates from the VRT's
GeoTransform and CRS from the SRS tag.
New module: _vrt.py with parse_vrt() and read_vrt() functions.
8 new tests: single tile, 2x1 mosaic, 2x2 mosaic, windowed read,
CRS propagation, nodata, read_vrt API, and XML parser unit test.
1. Band-first DataArray (CRITICAL): write_geotiff now detects (band, y, x) dimension order and transposes to (y, x, band). Prevents silent data corruption from rasterio-style arrays. 2. HTTP COG sub-byte support (CRITICAL): the COG HTTP reader now routes through _decode_strip_or_tile like the local readers, so 1-bit/4-bit/12-bit COGs over HTTP work correctly. 3. Dask VRT support (USEFUL): read_geotiff_dask detects .vrt files and reads eagerly then chunks, since VRT windowed reads need the virtual dataset's source layout. 4. VRT writer (USEFUL): write_vrt() generates a VRT XML file from multiple source GeoTIFFs, computing the mosaic layout from their geo transforms. Supports relative paths and CRS/nodata. 5. ExtraSamples tag (USEFUL): RGBA writes now include tag 338 with value 2 (unassociated alpha). Multi-band with >3 bands also gets ExtraSamples for bands beyond RGB. 6. MinIsWhite (USEFUL): photometric=0 (MinIsWhite) single-band files are now inverted on read so 0=black, 255=white. Integer values are inverted via max-value, floats via negation. 7. Post-write validation (POLISH): after writing, the header bytes are parsed to verify the output is a valid TIFF. Emits a warning if the header is corrupt. 8. Float16/bool auto-promotion (POLISH): float16 arrays are promoted to float32, bool arrays to uint8, instead of raising ValueError. 275 tests passing. 7 new tests for the fixes plus updated edge case tests.
Removes the rioxarray dependency from all example notebooks: - multispectral.ipynb: rioxarray.open_rasterio -> read_geotiff - classification-methods.ipynb: same - viewshed_gpu.ipynb: same - 25_GLCM_Texture.ipynb: rioxarray.open_rasterio COG read -> read_geotiff with window= and band= parameters. Also removes GDAL-specific env vars (AWS_NO_SIGN_REQUEST, etc.) since our reader doesn't use GDAL. Also updates reproject/_crs_utils.py to check attrs['crs'] and attrs['crs_wkt'] (xrspatial.geotiff convention) before falling back to .rio.crs (rioxarray). This means DataArrays from read_geotiff work directly with xrspatial.reproject without needing rioxarray installed. The rioxarray fallback is kept in _crs_utils.py for backwards compatibility with users who pass rioxarray-decorated DataArrays.
Both are now required (not optional): - matplotlib: needed for palette colormap (ListedColormap) and da.xrs.plot() with palette TIFFs - zstandard: needed for ZSTD compression (tag 50000), increasingly common in modern COGs This fixes the CI failures where these packages weren't installed.
read_geotiff_gpu() decodes tiled GeoTIFFs on the GPU and returns a CuPy-backed DataArray that stays on device memory. No CPU->GPU transfer needed for downstream xrspatial GPU operations (slope, aspect, hillshade, etc.). CUDA kernels implemented: - LZW decode: one thread block per tile, LZW table in shared memory (20KB per block, fast on-chip SRAM) - Predictor decode (pred=2): one thread per row, horizontal cumsum - Float predictor (pred=3): one thread per row, byte-lane undiff + un-transpose - Tile assembly: one thread per pixel, copies from decompressed tile buffer to output image Supports LZW and uncompressed tiled TIFFs. Falls back to CPU for unsupported compression types or stripped files. 100% pixel-exact match with CPU reader on all tested files (USGS LZW+pred3 3612x3612, synthetic LZW tiled). Performance: GPU LZW is comparable to CPU (~330ms vs 270ms for 3612x3612) because LZW is inherently sequential per-stream. The value is in keeping data on GPU for end-to-end pipelines without CPU->GPU transfer overhead. Future work: CUDA inflate (deflate) kernel would unlock the parallel decompression win since deflate tiles are much more common in COGs.
Implements RFC 1951 deflate decompression as a Numba @cuda.jit kernel for GPU-accelerated TIFF tile decoding. One thread block per tile, all tiles decompress in parallel. Supports all three deflate block types: - BTYPE=0: stored (no compression) - BTYPE=1: fixed Huffman codes - BTYPE=2: dynamic Huffman codes (most common in real files) Uses a two-level Huffman decode: - Fast path: 10-bit shared-memory lookup table (1024 entries) - Slow path: overflow array scan for codes > 10 bits (up to 15) Fixes the infinite loop bug where 14-bit lit/len codes exceeded the original 10-bit table size. Tested: 100% pixel-exact match on Copernicus deflate+pred3 COG (3600x3600, 16 tiles) vs CPU zlib. Performance: GPU inflate is ~20x slower than CPU zlib for this file size (16 tiles). Deflate is inherently sequential per-stream, so each thread block runs a long serial loop while most SMs sit idle. The value is keeping data on GPU for end-to-end pipelines. For files with hundreds of tiles, the parallelism would help more.
gpu_decode_tiles() now tries kvikio.nvcomp.DeflateManager for batch deflate decompression before falling back to the Numba CUDA inflate kernel. nvCOMP is NVIDIA's optimized batched compression library that decompresses all tiles in a single GPU API call. Fallback chain for GPU decompression: 1. nvCOMP via kvikio (if installed) -- optimized CUDA kernels 2. Numba @cuda.jit inflate kernel -- pure Python/Numba implementation 3. CPU zlib fallback -- if GPU decode raises any error kvikio is an optional dependency (pip install kvikio-cu12 or conda install -c rapidsai kvikio). When not installed, the Numba kernels are used transparently.
Fixed the nvCOMP C API ctypes binding to pass opts structs by value using proper ctypes.Structure subclasses. The previous byte-array approach caused the struct to be misinterpreted by nvCOMP. Working: nvCOMP ZSTD batch decompress (nvcompBatchedZstdDecompressAsync) - 100% pixel-exact match on all tested files - 1.5x end-to-end speedup on 8192x8192 ZSTD with 1024 tiles (GPU pipeline: 404ms vs CPU+transfer: 620ms) Not working on Ampere: nvCOMP deflate returns nvcompErrorNotSupported (status 11). Deflate GPU decompression requires Ada Lovelace or newer GPU with HW decompression engine. Falls back to the Numba CUDA inflate kernel on Ampere. nvCOMP is auto-detected by searching for libnvcomp.so in CONDA_PREFIX and sibling conda environments. When found, ZSTD tiles are batch-decompressed in a single GPU API call.
When kvikio is installed, read_geotiff_gpu() can read compressed tile bytes directly from NVMe SSD to GPU VRAM via GPUDirect Storage, bypassing CPU memory entirely: Normal: SSD -> CPU (mmap) -> cupy.asarray (CPU->GPU copy) With GDS: SSD -> GPU VRAM (direct DMA, no CPU involved) The full pipeline for a ZSTD COG with GDS + nvCOMP: SSD --(GDS)--> GPU compressed tiles --(nvCOMP)--> GPU decompressed --> GPU predictor decode --> GPU tile assembly --> CuPy DataArray Fallback chain in read_geotiff_gpu: 1. KvikIO GDS file read + nvCOMP batch decompress (fastest) 2. CPU mmap tile extract + nvCOMP batch decompress 3. CPU mmap tile extract + Numba CUDA kernels 4. CPU read_to_array + cupy.asarray transfer (slowest) Also adds: - gpu_decode_tiles_from_file(): accepts file path + offsets instead of pre-extracted bytes, enabling the GDS path - _try_nvcomp_from_device_bufs(): nvCOMP on tiles already in GPU memory (from GDS), avoiding a device-to-host round-trip - _apply_predictor_and_assemble(): shared GPU post-processing used by both GDS and mmap paths KvikIO is optional: conda install -c rapidsai kvikio GDS requires: NVMe SSD + NVIDIA kernel module (nvidia-fs)
- GDS tile read: added sync + verification after each pread to catch partial reads and CUDA errors early. Catches exception and tries to reset CUDA state before falling back. - gpu_decode_tiles: unsupported GPU codecs (ZSTD without nvCOMP, etc.) now decompress on CPU then transfer to GPU instead of raising ValueError. This keeps the predictor + assembly on GPU. - Fixes cudaErrorIllegalAddress from kvikio version mismatch (26.02 C lib vs 26.06 Python bindings) by catching the error gracefully instead of poisoning the GPU state.
nvCOMP deflate decompression now works on all CUDA GPUs by using backend=2 (CUDA software implementation) instead of backend=0 (DEFAULT, which tries hardware decompression first and fails on pre-Ada GPUs). Benchmarks (read + slope, A6000 GPU, nvCOMP via libnvcomp.so): Deflate: 8192x8192 (1024 tiles): GPU 769ms vs CPU 1364ms = 1.8x 16384x16384 (4096 tiles): GPU 2417ms vs CPU 5788ms = 2.4x ZSTD: 8192x8192 (1024 tiles): GPU 349ms vs CPU 404ms = 1.2x 16384x16384 (4096 tiles): GPU 1325ms vs CPU 2087ms = 1.6x Both codecs decompress entirely on GPU via nvCOMP batch API. No CPU decompression fallback needed when nvCOMP is available. 100% pixel-exact match verified.
Adds a GeoTIFF / COG I/O section to the feature matrix covering: - read_geotiff, write_geotiff, read_geotiff_gpu, VRT, open_cog - Compression codecs (deflate, LZW, ZSTD, PackBits, JPEG) - GPU decompression via nvCOMP (2.4x speedup at 16K x 16K) - Cloud storage, GDS, metadata preservation, sub-byte support - Overview resampling modes Updates Quick Start to use read_geotiff instead of synthetic data. Updates Notes on GDAL to reflect native reader capabilities. Updates Dependencies to list core and optional packages.
Sections now go from most-used to least-used in a typical workflow: 1. GeoTIFF / COG I/O (read your data first) 2. Surface (slope, aspect, hillshade -- the basics) 3. Hydrology (flow direction, accumulation, watersheds) 4. Flood (downstream from hydrology) 5. Multispectral (satellite imagery) 6. Classification (binning results) 7. Focal (neighborhood analysis) 8. Proximity (distance) 9. Zonal (zonal stats) 10. Reproject / Merge 11. Interpolation 12. Morphological 13. Fire 14. Raster / Vector Conversion 15. Utilities 16. Multivariate, Pathfinding, Diffusion, Dasymetric (specialized) Previously alphabetical, which put Classification first and buried Surface and Hydrology in the middle.
write_geotiff_gpu() compresses tiles on the GPU and writes a valid GeoTIFF. The CuPy array stays on device throughout -- only the compressed bytes transfer to CPU for file assembly. GPU pipeline: CuPy array → tile extraction (CUDA kernel) → predictor encode (CUDA kernel) → nvCOMP batch compress → CPU file assembly CUDA kernels added: - _extract_tiles_kernel: image → per-tile buffers (1 thread/pixel) - _predictor_encode_kernel: horizontal differencing (1 thread/row) - _fp_predictor_encode_kernel: float predictor (1 thread/row) - _nvcomp_batch_compress: deflate + ZSTD via nvCOMP C API Deflate write performance (tiled 256, A6000): 2048x2048: GPU 135ms vs CPU 424ms = 3.1x faster 4096x4096: GPU 302ms vs CPU 1678ms = 5.6x faster 8192x8192: GPU 1114ms vs CPU 6837ms = 6.1x faster GPU deflate is also 1.5-1.8x faster than rioxarray/GDAL at 4K+. All round-trips verified pixel-exact (deflate, ZSTD, uncompressed).
README: - Updated feature matrix: write_geotiff now shows Dask ✅ and CuPy 🔄 (fallback). Added write_geotiff_gpu and read_geotiff_dask rows. Updated VRT to show Dask support. - Added comprehensive benchmark tables for read (real-world + synthetic) and write (CPU vs GPU vs rioxarray) across all sizes and codecs. - 100% consistency verified across all tested files. Backend support for write_geotiff: - NumPy: direct write (existing) - Dask DataArray: .compute() then write (existing, now documented) - CuPy raw array: .get() to numpy then write (new) - CuPy DataArray: .data.get() then write (new) - Dask+CuPy: .compute().get() then write (new) - Python list: np.asarray() then write (existing) For GPU-native compression (no CPU transfer), use write_geotiff_gpu.
read_geotiff_gpu:
- New chunks= parameter returns a Dask+CuPy DataArray
- read_geotiff_gpu('dem.tif', chunks=512) decompresses on GPU
then chunks the result for out-of-core GPU pipelines
write_geotiff_gpu:
- Accepts Dask+CuPy DataArrays (.compute() then compress on GPU)
- Accepts Dask+NumPy DataArrays (.compute() then transfer to GPU)
- Accepts raw CuPy, numpy, or list inputs
All 7 input combinations verified:
read_geotiff_gpu -> CuPy DataArray (existing)
read_geotiff_gpu(chunks=N) -> Dask+CuPy DataArray (new)
write_geotiff_gpu(cupy_array) (existing)
write_geotiff_gpu(cupy_DataArray) (existing)
write_geotiff_gpu(dask_cupy_DataArray) (new)
write_geotiff_gpu(numpy_array) (auto-transfer)
write_geotiff_gpu(dask_numpy_DataArray) (auto-compute+transfer)
Also fixed write_geotiff CuPy fallback for raw arrays and
Dask+CuPy DataArrays (compute then .get() to numpy).
read_geotiff and write_geotiff now dispatch to the correct backend
automatically:
read_geotiff('dem.tif') # NumPy (default)
read_geotiff('dem.tif', gpu=True) # CuPy via nvCOMP
read_geotiff('dem.tif', chunks=512) # Dask lazy
read_geotiff('dem.tif', gpu=True, chunks=512) # Dask+CuPy
write_geotiff(numpy_arr, 'out.tif') # CPU write
write_geotiff(cupy_arr, 'out.tif') # auto-detects CuPy -> GPU write
write_geotiff(data, 'out.tif', gpu=True) # force GPU write
Auto-detection: write_geotiff checks isinstance(data, cupy.ndarray)
to decide whether to use GPU compression. Falls back to CPU if
cupy is not installed or nvCOMP fails.
read_vrt also supports gpu= and chunks= parameters for all four
backend combinations.
Users no longer need to call read_geotiff_gpu/write_geotiff_gpu
directly -- the main functions handle everything.
read_geotiff already accepts HTTP URLs, cloud URIs (s3://, gs://,
az://), local files, and VRT files. open_cog is now a thin
deprecated wrapper. Users just use read_geotiff for everything:
read_geotiff('https://example.com/cog.tif')
read_geotiff('s3://bucket/cog.tif')
read_geotiff('/local/dem.tif')
read_geotiff('mosaic.vrt')
All with gpu=, chunks=, window=, band= options.
Removed open_cog from the README feature matrix.
The public API is now: read_geotiff(source, ...) # read anything: file, URL, cloud, VRT write_geotiff(data, path, ...) # write any backend write_vrt(path, sources) # generate VRT mosaic XML read_geotiff auto-detects: - .vrt extension -> VRT reader - http:// / https:// -> COG range-request reader - s3:// / gs:// / az:// -> cloud via fsspec - gpu=True -> nvCOMP GPU decompression - chunks=N -> Dask lazy windowed reads Removed from __all__: open_cog (deprecated wrapper), read_vrt (called internally), read_geotiff_dask (use chunks=), read_geotiff_gpu / write_geotiff_gpu (use gpu=True). All these functions still exist for backwards compatibility but are no longer the recommended entry points.
* Add JPEG 2000 codec with optional nvJPEG2000 GPU acceleration (#1048) CPU path via glymur (same pattern as JPEG/Pillow and ZSTD/zstandard). GPU path via nvJPEG2000 ctypes bindings (same pattern as nvCOMP). Both are optional -- graceful fallback if libraries aren't installed. * Add JPEG 2000 test coverage and fix glymur numres constraint (#1048) - 14 new tests: codec roundtrip, TIFF write/read roundtrip, public API, availability checks, and ImportError fallback - Fix jpeg2000_compress: calculate numres from tile dimensions, remove pre-existing temp file before glymur write - Update test_edge_cases: use 'webp' as unsupported compression example since 'jpeg2000' is now supported * Add JPEG 2000 compression user guide notebook (#1048) Covers write/read roundtrip, file size comparison with deflate, multi-band RGB example, and GPU acceleration notes. * Update README: add JPEG 2000 / nvJPEG2000 to codec lists (#1048)
* Rename GeoTIFF API to xarray conventions (#1047) open_geotiff replaces read_geotiff, to_geotiff replaces write_geotiff. Adds .xrs.to_geotiff() accessor on DataArray and Dataset, and .xrs.open_geotiff() on Dataset for spatially-windowed reads. * Add tests for .xrs.to_geotiff and .xrs.open_geotiff accessors (#1047) * Update README and notebooks to open_geotiff / to_geotiff (#1047) * Add GeoTIFF I/O user guide notebook (#1047) * Rename old API refs in JPEG 2000 tests and notebook (#1047)
* Update README narrative to match current scope (#1045) The intro, GDAL caveat, and docs note were written when the library was much smaller. Now there are 100+ functions, native GeoTIFF I/O, 4-backend dispatch, and 33+ user guide notebooks. Updated the prose to say what the library actually does instead of underselling it. * Fix GPU emoji shortcode, add read/reproject/write example (#1045) :gpu: doesn't render on GitHub -- switched to :desktop_computer:. Added a quick start snippet showing the read_geotiff -> reproject -> write_geotiff workflow. * Use import xrspatial as xrs, prioritize accessor methods (#1045) Switched quick start to 'import xrspatial as xrs', moved the read/reproject/write flow into the main example using the .xrs accessor, and simplified the standalone function example to use the xrs namespace. * Fix output grid computation for reproject (#1045) Three fixes to _grid.py: - Resolution estimation now transforms each axis independently and uses the geometric mean for square pixels (was diagonal step) - Edge sampling increased from 21 to 101 points plus a 21x21 interior grid for better bounds coverage - ceil() replaced with round() for dimension calculation to prevent floating-point noise from adding spurious rows/columns * Add Numba JIT and CUDA projection kernels for reproject (#1045) Ports six projections from the PROJ library to Numba (CPU) and Numba CUDA (GPU), bypassing pyproj for common CRS pairs: - Web Mercator (EPSG:3857) -- spherical, 3 lines per direction - Transverse Mercator / UTM (326xx, 327xx, 269xx) -- 6th-order Krueger series (Karney 2011), closed-form forward and inverse - Ellipsoidal Mercator (EPSG:3395) -- Newton inverse - Lambert Conformal Conic (e.g. EPSG:2154) -- Newton inverse - Albers Equal Area (e.g. EPSG:5070) -- authalic latitude series - Cylindrical Equal Area (e.g. EPSG:6933) -- authalic latitude series CPU Numba kernels are 6-9x faster than pyproj. CUDA kernels are 40-165x faster. Unsupported CRS pairs fall back to pyproj. _transform_coords now tries Numba first, then pyproj. The CuPy chunk worker tries CUDA first, keeping coordinates on-device. * Add reproject benchmark vs rioxarray (#1045) Compares xrspatial (approx, exact, numba) against rioxarray on synthetic grids and real-world GeoTIFFs. Measures both performance (median of 5 runs) and pixel-level consistency (RMSE, R^2, NaN agreement) by forcing both libraries onto the same output grid. * Update README with Numba/CUDA projection table (#1045) Updated Reproject description and added a table listing the six projections with native Numba CPU and CUDA GPU kernels. * Add Sinusoidal and Polar Stereographic Numba kernels (#1045) Three new CPU Numba projections for bathymetric/ocean use cases: - Sinusoidal (ellipsoidal): MODIS ocean/land products. Uses pj_mlfn meridional arc length with 5th-order series. Forward: sub-micrometer vs pyproj. Roundtrip: exact. - Polar Stereographic (N/S): IBCAO Arctic (3996/3413), IBCSO Antarctic (3031), UPS. Iterative inverse (15 iter max). Forward: sub-nanometer. Roundtrip: exact. LAEA implemented but disabled in dispatch pending investigation of ~940m oblique-mode error (kernels are in the file for future fix, just not wired into the dispatch table). * Fix LAEA xmf/ymf swap, re-enable in dispatch (#1045) The oblique-mode LAEA had xmf and ymf swapped (rq/dd vs rq*dd). Research agent found the bug by comparing against PROJ's laea.cpp source. Forward accuracy is now sub-millimeter vs pyproj. * Add generic tmerc dispatch for State Plane zones (#1045) State Plane zones that use Transverse Mercator (Maine, New Hampshire, Wisconsin, etc.) now hit the Numba fast path. Uses the same Krueger 6th-order series as UTM, with a Zb offset for non-zero lat_0. Meter-based zones only; US survey foot zones fall back to pyproj. * Support US survey foot and international foot units (#1045) State Plane zones in us-ft (136 zones) and ft (28 zones) now use the Numba fast path. The Krueger/LCC kernels compute in metres internally, then divide by the unit conversion factor (0.3048006 for us-ft, 0.3048 for ft) on output. x_0/y_0 from PROJ4 dicts are already in metres regardless of the units parameter. * Add projection benchmark table to README, fix dispatch for custom CRS (#1045) - Benchmark table showing Numba vs pyproj for all 12 supported projections - Fixed dispatch to work with custom PROJ strings (no EPSG code), needed for Sinusoidal and other non-registered CRS definitions - Fixed _utm_params to handle None epsg_code * Add GPU column and inline speedups to projection benchmark table (#1045) * Add CUDA kernels for Sinusoidal, LAEA, Polar Stere, State Plane (#1045) All 12 projections now have GPU CUDA kernels. Performance on A6000: - Sinusoidal: 18ms (56x vs pyproj) - LAEA Europe: 18ms (92x) - Polar Stere: 57ms (64-67x) - State Plane tmerc: 23ms (88x) - State Plane lcc ftUS: 36ms (124x) - LCC France: 39ms (78x) All bit-exact against CPU Numba kernels. Updated README benchmark table and projection support matrix. * Guard Numba dispatch against non-WGS84 datums (#1045) CRS on non-WGS84/GRS80 ellipsoids (Airy for BNG, Clarke 1866 for NAD27, Bessel for Tokyo, etc.) now fall back to pyproj instead of silently skipping the datum transformation. Without this, BNG had ~100m error, NAD27 ~24m, Tokyo ~900m. Each _*_params() function now checks _is_wgs84_compatible_ellipsoid() before returning parameters. EPSG-specific paths (UTM, Web Mercator) were already safe since they only match WGS84/NAD83 codes. * Add NAD27 datum support via geocentric Helmert shift (#1045) NAD27 (EPSG:4267) sources and targets now go through the Numba fast path. After the projection kernel runs in WGS84, a 3-parameter geocentric Helmert transform (dx=-8, dy=160, dz=176 for Clarke 1866) shifts coordinates to/from NAD27. Accuracy: mean 2.9m, p95 5.8m vs pyproj across CONUS. This matches PROJ's own accuracy when NADCON grids aren't installed. The framework is extensible to other datums by adding entries to _DATUM_PARAMS. Also broadened geographic CRS detection from WGS84/NAD83-only to include NAD27, so NAD27 -> UTM and NAD27 -> State Plane dispatch correctly. * Add CUDA resampling kernels for end-to-end GPU reproject (#1045) Native CUDA nearest, bilinear, and cubic (Catmull-Rom) resampling kernels replace cupyx.scipy.ndimage.map_coordinates. When the CUDA projection path produces on-device coordinates, the entire pipeline now stays on GPU with no CPU roundtrip. Full reproject pipeline (2048x2048, bilinear, 4326->UTM): GPU end-to-end: 78ms CPU Numba: 161ms Speedup: 2.1x * Add merge benchmark: xrspatial vs rioxarray (#1045) Merges 4 overlapping WGS84 tiles and compares timing and pixel-level consistency against rioxarray.merge_arrays. Baseline results (xrspatial is currently slower on merge): 512x512 tiles: 59ms vs 56ms (1.1x) 1024x1024: 293ms vs 114ms (2.6x) 2048x2048: 2.52s vs 656ms (3.8x) Consistency: RMSE < 0.012, NaN agreement > 99.8%. * Fast same-CRS merge and early-exit in numba dispatch (#1045) Two optimizations that make merge 4.5-7.3x faster: 1. Same-CRS tiles skip reprojection entirely. When source and target CRS match, tiles are placed into the output grid by direct coordinate alignment (array slicing). Falls back to full reprojection if resolutions differ by >1%. 2. try_numba_transform now bails before allocating coordinate arrays when neither CRS side is a supported geographic system. This saved ~100ms per tile for unsupported pairs. Merge benchmark (4 overlapping WGS84 tiles): 512x512: 13ms (was 59ms, now 2.3x faster than rioxarray) 1024x1024: 48ms (was 293ms, now 2.6x faster than rioxarray) 2048x2048: 344ms (was 2520ms, now 1.6x faster than rioxarray) * Replace coordinate-only benchmarks with end-to-end reproject/merge tables (#1045) README now shows full pipeline times (transform + resampling) and merge times, both compared against rioxarray. More useful than the previous coordinate-transform-only table since users care about total wall time. * Dask+CuPy reproject: single-pass GPU instead of per-chunk (#1045) For dask+cupy inputs, eagerly compute the source to GPU memory and run the in-memory CuPy reproject in a single pass. This avoids the per-chunk overhead of 64+ dask.delayed calls, each creating a pyproj Transformer and launching small CUDA kernels. Before: 958ms (64 delayed chunks, 512x512 each) After: 43ms (single CuPy pass, pixel-exact same output) Speedup: 22x The output is a plain CuPy array. For truly out-of-core GPU data that doesn't fit in GPU memory, the old dask.delayed path remains available by passing the data as dask+numpy. * Chunked dask+cupy reproject without full-source eager compute (#1045) Replaces the eager .compute() approach with a chunked GPU pipeline that fetches only the needed source window per output chunk. This handles sources larger than GPU memory while still being 8-20x faster than the old dask.delayed path. The key optimizations vs dask.delayed: - CRS objects and transformer created once (not per chunk) - CUDA projection + native CUDA resampling per chunk - Default 2048x2048 GPU chunks (not 512x512) - Sequential loop avoids dask scheduler overhead Performance (4096x4096 WGS84 -> UTM, bilinear): CuPy single pass: 34ms Dask+CuPy (2048): 49ms (was 958ms) Dask+CuPy (512): 71ms Dask+CuPy (256): 124ms All chunk sizes are pixel-exact vs plain CuPy (max_err < 1e-11). * Add NADCON grid-based datum shift for sub-meter NAD27 accuracy (#1045) Vendored two NOAA shift grids into the package (306KB total): - us_noaa_conus.tif: NADCON classic (121x273, 0.25° resolution) - us_noaa_nadcon5_nad27_nad83_1986_conus.tif: NADCON5 (105x237) The grid loader checks the vendored directory first, then a user cache, then downloads from the PROJ CDN as a last resort. Numba JIT bilinear interpolation applies the lat/lon arc-second offsets per pixel, with an iterative inverse for target->source direction. When a grid covers the data, it replaces the Helmert shift (which had ~3-5m accuracy). The grid gives sub-meter accuracy matching PROJ with NADCON grids installed. Points outside grid coverage fall back to Helmert automatically. * Vendor 14 datum shift grids for worldwide sub-meter accuracy (#1045) Shipped 8.4MB of NOAA/NTv2 shift grids covering: US: NAD27 CONUS (NADCON + NADCON5), Alaska, Hawaii, Puerto Rico UK: OSGB36 -> ETRS89 (Ordnance Survey OSTN15) Australia: AGD66 -> GDA94 (NT region) Europe: Germany (DHDN), Austria (MGI), Spain (ED50, eastern coast), Netherlands (RD), Belgium (BD72), Switzerland (CH1903), Portugal (D73) Added Helmert fallback parameters for all 14 datums plus Tokyo. Grid lookup automatically selects the best grid covering a point, falling back to Helmert for regions without grid coverage. All grids are Public Domain or Open Government Licence. * Add vertical datum support with vendored EGM96 geoid (#1045) New public API for ellipsoidal <-> orthometric height conversion: geoid_height(lon, lat) # geoid undulation N (metres) geoid_height_raster(da) # N for every pixel ellipsoidal_to_orthometric(h, ...) # GPS height -> map height orthometric_to_ellipsoidal(H, ...) # map height -> GPS height depth_to_ellipsoidal(depth, ...) # chart datum depth -> ellipsoidal ellipsoidal_to_depth(h, ...) # ellipsoidal -> chart datum depth Vendored EGM96 global geoid model (2.6MB, 15-arcmin / ~28km grid, 721x1440 pixels). EGM2008 (77MB, 2.5-arcmin) available via CDN download on first use. Numba JIT bilinear interpolation with longitude wrapping for the global grid. Performance: 80 Mpix/s on CPU (1M points in 12ms). * Add time-dependent ITRF frame transforms (#1045) 14-parameter Helmert (7 static + 7 rates) for converting between ITRF frames at a given observation epoch. Parameters parsed from the PROJ data files shipped with pyproj. Available frames: ITRF2000, ITRF2008, ITRF2014, ITRF2020 (and all intermediate frames back to ITRF89). Usage: lon2, lat2, h2 = itrf_transform( -74.0, 40.7, 10.0, src='ITRF2014', tgt='ITRF2020', epoch=2024.0, ) Typical shifts: 2-4mm horizontal, 1-3mm vertical between ITRF2014 and ITRF2020 at epoch 2024. The rates capture tectonic plate motion (~mm/year) which accumulates over years. Numba JIT parallelized for batch transforms. * 7-parameter Helmert and 6-term authalic latitude series (#1045) Helmert upgrade (3-param -> 7-param Bursa-Wolf): - Added rx/ry/rz rotations (arcsec) and ds scale (ppm) to the geocentric datum shift kernel - Updated OSGB36, DHDN, MGI, ED50, BD72 with published 7-param values from the EPSG dataset - OSGB36 Helmert fallback improved from 15.73m to 1.44m vs pyproj - Same kernel handles 3-param datums (rx=ry=rz=ds=0) Authalic latitude series (3-term -> 6-term): - Extended _authalic_apa to 6 coefficients (10th-order in e²) - Updated _authalic_inv and CUDA _d_authalic_inv to evaluate 5 terms - Theoretical improvement: sub-mm for the series itself, though the q->beta->phi roundtrip error is dominated by the asin(q/qp) step at high latitudes (~4.8m at 89°, <0.1m at mid-latitudes) * Add oblique stereographic and oblique Mercator kernels (disabled) (#1045) Numba kernels for two additional projections: - Oblique Stereographic: Gauss conformal sphere + stereographic. Kernel roundtrips perfectly but forward differs from PROJ's specific Gauss-Schreiber conformal mapping by ~50km. Needs alignment with PROJ's sterea.cpp double-projection approach. - Oblique Mercator (Hotine): rotation matrix + Mercator on the oblique cylinder. Forward has errors in the u/v -> x/y rotation. Needs closer alignment with PROJ's omerc.cpp variant handling. Both kernels are implemented and compile but disabled in the dispatch table pending accuracy fixes. They fall through to pyproj. Also: Equidistant Conic skipped -- has zero EPSG definitions in the PROJ database, essentially unused in practice. * Fix oblique stereographic with Gauss conformal sphere (#1045) The oblique stereographic now uses the correct PROJ double-projection: 1. Gauss conformal mapping: phi -> chi via scaling factor C = sqrt(1 + e²cos⁴φ₀/(1-e²)) and normalization constant K 2. Standard stereographic on the conformal sphere Forward accuracy: sub-nanometre vs pyproj. Roundtrip: exact (1.4e-14 degrees). Also fixed R scaling: R_metric = a * k0 * R_conformal, where R_conformal is the dimensionless conformal radius from PROJ. Oblique Mercator (Hotine) remains disabled -- the u/v rotation and variant handling need more work. * 2D Numba kernels for LCC/tmerc: eliminate tile/repeat + fuse unit conv (#1045) New lcc_inverse_2d and tmerc_inverse_2d kernels take 1D x/y arrays and write directly to 2D output, avoiding: - np.tile (199ms for 4096x4096) - np.repeat (357ms for 4096x4096) - Separate unit division pass (115ms for ftUS) The unit conversion (feet -> metres) is fused into the inner loop, operating on scalars instead of 16.8M-element arrays. California zone 5 ftUS (4096x4096, bilinear): Before: 915ms (0.9x vs rioxarray) After: 712ms (1.2x vs rioxarray) * Fix longitude wrapping in all projection inverses (#1045) Added _norm_lon_rad() and applied it to all inverse functions that compute lam + lon0. Without normalization, projections with non-zero lon0 (e.g. EPSG:3413 Arctic Stere with lon0=-45°) could return longitudes outside [-180, 180], causing 360° errors and false NaN pixels in the source lookup. Polar Stere N (EPSG:3413) consistency: Before: RMSE=8.29, NaN agree=90.4%, 1.1M extra NaN After: RMSE=0.008, NaN agree=100%, 79 extra NaN (edge pixels) * Relax resampling boundary check to match GDAL behavior (#1045) Changed out-of-bounds threshold from -0.5/sh-0.5 to -1.0/sh in all resampling kernels (nearest, bilinear, cubic, and CUDA). Pixels within one pixel of the source edge are now clamped to the nearest valid pixel instead of being set to nodata. This matches GDAL/rasterio's boundary convention, fixing 2568 false-NaN pixels at the edges of LAEA Europe reprojections. LAEA NaN agreement: 99.8% -> 100.0% All other projections: unchanged or improved to 100.0% * Match GDAL's bilinear weight renormalization (#1045) Changed bilinear resampling (CPU Numba + CUDA) from clamp-and-use-all to skip-and-renormalize, matching GDAL's GWKBilinearResample4Sample: - Out-of-bounds neighbors: skipped, weight redistributed to valid ones (was: clamped to edge pixel, counted at full weight) - NaN neighbors: skipped, interpolated from remaining valid pixels (was: any NaN contaminates the whole output pixel) This eliminates the one-pixel NaN halo around nodata regions that the old approach produced, and gives mathematically exact results on linear surfaces at raster boundaries. The ~0.017 RMSE vs rioxarray on gradient rasters is unchanged -- it comes from sub-pixel floating-point coordinate differences in the interior, not edge handling. * Harden reproject: thread safety, NaN crash, uint8 overflow (#1045) Thread safety: - Added threading.Lock to global caches in _datum_grids.py, _vertical.py, and _itrf.py. Concurrent callers no longer race on grid loading or ITRF parameter parsing. All-NaN raster crash: - np.nanmin on an all-NaN array returns NaN; int(NaN) is undefined. Added np.isfinite guards in both numpy and cupy chunk workers. uint8 cubic overflow: - Cubic resampling can ring outside [0, 255] on sharp edges. Added np.clip clamping before the dtype cast for all integer source types (uint8, int16, etc.) across nearest/bilinear/cubic. Geoid poles: - _interp_geoid_point now returns NaN (not 0.0) outside the grid's latitude range, preventing silent wrong values at the poles. Exception narrowing: - Bare except Exception:pass blocks around Numba/CUDA fast paths narrowed to except (ImportError, ModuleNotFoundError). Actual bugs now propagate instead of silently falling back to pyproj. New tests: - test_reproject_1x1_raster: single-pixel source - test_reproject_all_nan: 100% NaN source - test_reproject_uint8_cubic_no_overflow: cubic on uint8 sharp edge * Fix cubic NaN, add merge tests, validate grids, improve docs (#1045) Cubic NaN handling: - When any of the 16 Catmull-Rom neighbors is NaN, falls back to bilinear with weight renormalization instead of returning nodata. Eliminates the one-pixel nodata halo around NaN regions that cubic resampling previously produced. Merge strategy tests: - Added end-to-end tests for last, max, min strategies (were only tested at the internal _merge_arrays_numpy level). Datum grid validation: - load_grid() now validates band shapes match, grid is >= 2x2, and bounds are sensible. Invalid grids return None (Helmert fallback) instead of producing garbage. Documentation: - reproject() and merge() docstrings now note output CRS is WKT format in attrs['crs'], and merge() documents CRS selection when target_crs=None. * Integrate vertical datum shift into reproject() (#1045) New parameters src_vertical_crs and tgt_vertical_crs on reproject() enable automatic vertical datum transformation during reprojection: reproject(dem, 'EPSG:32633', src_vertical_crs='EGM96', # input is MSL heights tgt_vertical_crs='ellipsoidal') # want GPS heights Supported vertical CRS: 'EGM96', 'EGM2008', 'ellipsoidal'. Implementation: - After horizontal reprojection, output pixel coordinates are inverse-projected to geographic (lon/lat) for the geoid lookup - Geoid undulation N is interpolated from the vendored EGM96 grid - Heights are adjusted: h = H + N (ortho->ell) or H = h - N (ell->ortho) - Processes in 128-row strips to bound memory (12MB peak vs 768MB for a 4096x4096 raster) - Zero roundtrip error (ortho -> ell -> ortho recovers exact input) Also handles cross-geoid transforms (EGM96 -> EGM2008) by applying both shifts: H2 = H1 + N1 - N2. * Skip pyproj Transformer in chunk worker when Numba handles transform (#1045) For supported CRS pairs, the chunk worker now tries the Numba fast path BEFORE creating a pyproj.Transformer. If Numba succeeds (which it does for all 10+ supported projections), the Transformer is never created, saving ~15ms of pyproj setup per chunk. Before: 2 Transformer.from_crs calls per reproject (grid + chunk) After: 1 call (grid only, ~500 points for boundary estimation) The grid computation still uses pyproj for boundary transforms since it's only ~500 points and runs once. The per-pixel transform (millions of points) is handled entirely by Numba. * Button up reproject: README, benchmarks, write_geotiff WKT fix (#1045) README Reproject section updated: - All 11 projection families listed (added oblique stereographic) - Full pipeline benchmark table (read+reproject+write, all backends) - Datum shift coverage (14 grids, 10 Helmert fallbacks) - Vertical datum support (EGM96/EGM2008, integrated into reproject) - ITRF time-dependent frame transforms - pyproj usage documented (metadata only, Numba does the math) - Merge performance table updated benchmarks/reproject_benchmark.md: - 254-line comprehensive benchmark document with 6 sections - Full pipeline: NumPy 2.7s, CuPy 348ms, rioxarray 418ms - 13 projections tested for accuracy vs pyproj - Datum, vertical, ITRF coverage documented - All numbers from actual benchmark runs write_geotiff WKT fix: - reproject() stores CRS as WKT string in attrs['crs'] - write_geotiff assumed integer EPSG code, crashed with TypeError - Added isinstance check to parse WKT via _wkt_to_epsg() * Update benchmark with warm-kernel numbers (#1045) Separated reproject-only from full-pipeline timing. With warm Numba/CUDA kernels: - CuPy reproject: 73ms (2.0x faster than rioxarray) - rioxarray reproject: 144ms - NumPy reproject: 413ms Full pipeline (read+reproject+write) is dominated by I/O for compressed GeoTIFFs, where rioxarray's C-level rasterio beats our Python/Numba reader. Added note about ~4.5s JIT warmup on first call. * Parallel tile compression + ZSTD default: 13x faster writes (#1045) Three optimizations to the GeoTIFF writer: 1. Default compression changed from deflate to ZSTD: Same file size (40MB), 6x faster single-threaded compression. ZSTD is the modern standard; deflate still available via parameter. 2. Parallel tile compression via ThreadPoolExecutor: Tiles are independent, and zlib/zstd/LZW all release the GIL. Uses os.cpu_count() threads. Falls back to sequential for uncompressed or very few tiles (< 4). 3. Optimized uncompressed path: Pre-allocates contiguous buffer for all tiles. Combined results (3600x3600 float32): Write with new default (zstd parallel): 101ms (was 1388ms deflate sequential) Write deflate (parallel): 155ms (was 1388ms) vs rasterio: zstd 2.0x faster, deflate 3.0x faster Full pipeline (read + reproject + write): NumPy: 890ms (was 2907ms) Also fixed write_geotiff crash when attrs['crs'] contains a WKT string (produced by reproject()) -- added isinstance check to parse WKT via _wkt_to_epsg(). * Parallel tile decompression in GeoTIFF reader (#1045) Tile decompression (deflate, LZW, ZSTD) now runs in parallel using ThreadPoolExecutor, same approach as the writer. zlib, zstandard, and Numba LZW all release the GIL. Read performance (Copernicus 3600x3600 deflate): Before: 291ms (sequential) After: 101ms (parallel) -- 2.9x faster rasterio: 189ms -- we're now 1.9x FASTER than rasterio Full pipeline improvement (read + reproject + write): NumPy: 2907ms -> 697ms (4.2x faster total) * Support multi-band (RGB/RGBA) raster reprojection (#1045) Multi-band rasters (y, x, band) now reproject correctly: - Each band is reprojected independently using shared coordinates (coordinate transform computed once, reused for all bands) - Output preserves the band dimension name and coordinates - Works with any dtype (float32, uint8 with clamping, etc.) - Custom band dim names (e.g. 'channel') preserved Also fixed spatial dimension detection to use name-based lookup (_find_spatial_dims) instead of hardcoded dims[-2]/dims[-1], which failed for 3D rasters where the band dim was last. Previously crashed with TypingError on 3D input. * Add 14 edge case tests for reproject (#1045) New TestEdgeCases class covering: - Multi-band RGB/RGBA reprojection (float32, uint8) - Antimeridian crossing - Y-ascending coordinate order - Checkerboard NaN pattern (bilinear renormalization) - UTM -> geographic (reverse projection direction) - Projected -> projected (LCC -> UTM) - Sentinel nodata (-9999) - Integer EPSG as target CRS - Explicit resolution and width/height parameters - Merge with non-overlapping tiles - Merge with single tile * Prevent OOM on large datasets in reproject and merge (#1045) Three safeguards for datasets that exceed available RAM: 1. Auto-chunk large non-dask inputs (reproject): If the source array exceeds 512MB, automatically wrap it in dask.array with the configured chunk_size (default 512x512). This routes it through the chunked dask path instead of the in-memory path that would call .values and OOM. 2. Auto-promote merge to dask path: If the combined output size (output_shape * n_tiles * 8 bytes) exceeds 512MB, use the dask merge path even if no input is dask. This prevents _merge_inmemory from calling .values on each tile. 3. Cap source window size in chunk workers: If a single output chunk maps to a source window larger than 64 Mpixels (~512MB for float64), return nodata instead of materializing the window. This prevents extreme projections (e.g. polar stereographic edge pixels mapping to the entire source hemisphere) from OOMing individual chunk workers. A 30TB dataset with 16GB RAM would now: - Auto-chunk into 512x512 tiles - Process each tile independently (~2MB working memory per tile) - Never materialize more than 512MB in a single operation * Streaming reproject for datasets that exceed dask graph limits (#1045) For a 30TB raster at 2048x2048 chunks, dask's task graph would be 1.9GB -- larger than many machines' RAM. The streaming path bypasses dask entirely and processes output tiles in a sequential loop: for each output tile: compute source coordinates (Numba) read source window (lazy slice, no full materialization) resample write tile to output array free tile Memory usage: O(tile_size^2) per tile, ~16MB at 2048x2048. No graph overhead. No scheduler overhead. The routing logic: - Source < 512MB: in-memory (fastest) - Source > 512MB, graph < 1GB: auto-chunk to dask (parallel) - Source > 512MB, graph > 1GB: streaming (bounded memory) The streaming path produces results identical to the in-memory path (max error ~5e-13, floating-point noise only). * Add max_memory parameter for parallel streaming reproject (#1045) The streaming path (for datasets > ~1TB) now uses ThreadPoolExecutor with bounded concurrency based on the max_memory budget: reproject(huge_raster, 'EPSG:3857', max_memory='4GB') The budget controls how many output tiles can be processed in parallel. Numba kernels release the GIL, so threads give real parallelism. Memory stays bounded: max_memory='4GB', tile=2048x2048: ~42 concurrent tiles max_memory='16GB', tile=2048x2048: ~170 concurrent tiles Accepts int (bytes) or strings: '512MB', '4GB', '1TB'. Default is 1GB when not specified. On a 512x512 test with 256x256 tiles: Sequential (32MB budget): 233ms Parallel (4GB budget): 24ms -- 10x faster, identical output * Distributed streaming via dask.bag for multi-worker clusters (#1045) When dask.distributed is active, the streaming path uses dask.bag to distribute tile batches across workers instead of processing everything in one process: Local (no cluster): ThreadPoolExecutor within one process, max_memory bounded Distributed (dask cluster active): 1. Partition 2M tiles into N batches (one per worker) 2. dask.bag.from_sequence(batches, npartitions=N) 3. bag.map(process_batch) -- each worker gets its batch 4. Within each worker, ThreadPoolExecutor for intra-worker parallelism (Numba releases GIL) 5. Assemble results Graph size comparison for 30TB: Old dask.array approach: 1,968,409 nodes (1.9GB graph, OOM) New dask.bag approach: 4-64 nodes (one per worker) Each worker's memory bounded by max_memory parameter. Auto-detects distributed via get_client(). * Rename read_geotiff/write_geotiff to open_geotiff/to_geotiff (#1045) Aligns with the base branch rename (PR #1056) to match xarray conventions (open_* for readers, to_* for writers). Updated in: - xrspatial/geotiff/__init__.py (function defs, __all__) - xrspatial/reproject/_datum_grids.py, _vertical.py (grid loading) - README.md (all code examples and feature matrix) - benchmarks/reproject_benchmark.md - All geotiff test files (test_cog.py, test_edge_cases.py, test_features.py, bench_vs_rioxarray.py) All 343 tests pass (271 geotiff + 72 reproject).
JPEG is now a supported codec, so passing it with float32 data hits the dtype validation before the unsupported-compression check. Use bzip2 (actually unsupported) to test the right error path.
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.
Summary
open_geotiff(),to_geotiff(),da.xrs.plot()accessorTest plan
pytest xrspatial/geotiff/tests/(reader, writer, compression, JPEG 2000, edge cases, accessor)pytest xrspatial/tests/test_reproject.py(datum transforms, edge cases, multiband)