Recommended GDAL parity work

Initial assessment date: 2026-10-24. Based on the checkout before PR #239 and official GDAL documentation. Priorities and effort are engineering judgments, not measured estimates. Implementation is tracked in #238 and the roadmap.

Existing foundation

Terrio already implements common raster/vector formats, cloud and archive access, COGs, window reads, streaming vector conversion, streaming warp/mosaic/terrain outputs, GCP/RPC transformers, the full warp resampler list, STAC/STACIT/STACTA, GRIB, netCDF/HDF5/Zarr, Arrow, topology, contours and zonal statistics. These are documented in the roadmap, formats and CLI documentation. Recommendations below extend these implementations. The optional GDAL backend provides additional format access; it does not constitute native implementation of the proposed operations.

Deliberate TIFF mask difference

Issue #748 requests that a TIFF image's own reduced-resolution flag select its transparency mask, including the first image of an external overview file. Terrio honors a first-image/mask pair marked 1/5; GDAL 3.13.3 ignores that first mask. This applies to existing GDAL-written external pyramids as well as TIFFs opened directly and intentionally changes their masked reads.

The frozen masked-external cubic 26×41 reference pins 488 differing values out of 3,198 against unmodified GDAL output. Terrio is checked against a second GDAL control whose first image/mask flags are normalized to 0/4 without changing stored samples. Both references remain checked; this behavior is not described as GDAL pixel parity.

Prioritized additions

1. Fill nodata

Implemented in PR #239: a shared Rust operation exposed through CLI/MCP, Python and JavaScript, with GDAL comparison tests. The following describes the gap identified in the initial assessment.

Before PR #239, the CLI command inventory and core operations contained no fill-nodata operation. The rasterio-compatible CLI explicitly rejects fillnodata in rio/cli.py. The xarray accessor's interpolate_na is an existing interpolation surface, but does not supply this GDAL-compatible core operation.

The shared Rust operation repairs holes in DEMs and continuous imagery before contours, terrain analysis or reprojection. GDAL supports inverse-distance or nearest filling, an explicit mask, search distance and smoothing. GDAL fill-nodata.

Acceptance: compare pixels with GDAL for interior holes, edge holes, disconnected valid regions, NaN nodata and explicit masks. Valid source pixels must be preserved.

2. Overview-aware warping

Implemented in PR #242: the warper selects AUTO, AUTO-n, NONE or a forced source level after planning the output grid from full-resolution metadata. CLI/MCP, Python paths/lazy handles and JavaScript bytes/datasets retain stored pyramids. The conformance suite checks GDAL's selected pixels and measures actual HTTP response bytes for a USGS COG. Mosaics and cutlines currently load each selected level before compositing; single-source warps read only the required source windows. GDAL warp.

Acceptance: compare selected levels and output pixels with GDAL under matching settings; measure transferred bytes for remote COG downsampling. Preserve an explicit full-resolution path because overview generation changes sampled values.

3. Geolocation-array warping

Implemented in PR #241: the warper uses indexed bilinear geolocation cells for scientific swaths, including subsampled arrays, pixel-center conventions, invalid coordinates and antimeridian crossings. Rust, CLI/MCP, Python and JavaScript expose the capability. The supported subset and in-memory indexing are documented. Sentinel-5P and full-resolution valid arrays match GDAL exactly; subsampled/invalid-node cases retain 9–32 boundary-mask differences per 160×100 band, with common-valid values identical. GDAL warp.

Acceptance: use real swath products already represented in the scientific-format fixtures; compare bounds, valid-data coverage and resampled values. Include subsampled coordinate arrays, pixel-center conventions, invalid coordinates and antimeridian crossings.

4. Richer VRT interoperability

The VirtualRaster reader now extends existing simple/complex/averaged sources with serialized warped VRTs (affine/reprojection, GCP/TPS/homography/RPC), selected scalar derived functions, full/separable kernel filters and lazy masks. Rust parses/serializes the descriptions; CLI/MCP and Python write warped VRT files without materializing pixels, and Python/JavaScript accept uploaded source maps. Whole/cropped GDAL cross-reads are registered in the vrt conformance suite; a counted source regression checks bounded pixel reads. See the supported subset for remaining serialized cutline/alpha/geolocation and external-function limitations. Pansharpened VRTs follow the pansharpening implementation. GDAL VRT.

Acceptance: open GDAL-generated VRT fixtures and compare window reads, masks and pixels; verify Terrio-generated files in GDAL. Keep evaluation lazy and implement built-in operations before considering arbitrary user-supplied pixel functions.

5. Richer vector SQL

Implemented in PR #245: the documented SQL subset supports computed SELECT columns, single-field DISTINCT, COUNT/SUM/AVG/MIN/MAX summaries and same-dataset attribute LEFT JOINs. Duplicate join keys use OGR's first matching secondary row; null keys and unmatched rows retain the primary feature with null joined fields. Rust, CLI/MCP, Python and JavaScript share the subset. Queries retain primary geometry except DISTINCT and summaries. The grouped/spatial extension now supports attribute GROUP BY, selected output ordering and the six planar/EPSG spatial functions documented in the SQL guide, with GDAL SQLite conformance comparisons. Explicit geometry expressions replace the primary geometry; GROUP BY discards it. GDAL OGR SQL.

Acceptance: compare values, nulls, field types and geometry retention with OGR SQL, including mixed-case fields and duplicate join keys. Preserve streaming for operations that do not require materialization.

6. Viewshed

Implemented in PR #240: shared Rust single-observer binary visibility and minimum visible height, with observer/target heights, maximum distance, curvature/refraction and projected-unit handling. CLI/MCP, Python and JavaScript expose the operation. Real projected DEM comparisons are registered in conformance; numeric input nodata is treated as elevation, as in GDAL, and nonfinite elevations must be filled or cropped first. Cumulative visibility is available with --om ACCUM --os N, matching GDAL 3.13.3 normalization and its disabled curvature in cumulative mode. GDAL viewshed.

Acceptance: compare with GDAL on projected DEMs, hills/ridges, maximum distances and curvature settings. Define input nodata behavior explicitly: GDAL documents that it does not specially process input nodata cells.

7. Indexed raster catalogs

Terrio has mosaics and STAC mosaics, but its native driver registry has no GTI driver. Add a spatially indexed raster catalog using existing GeoPackage or FlatGeobuf readers. Store source paths and footprints, query only intersecting sources, and expose ordering and cache limits. GDAL's GTI driver uses vector indexes for large virtual mosaics and supports mixed source CRSs and compositing order. This is a scale-oriented extension to mosaics. GDAL GTI.

Acceptance: cross-read tile indexes with GDAL; compare overlap, alpha and nodata behavior. Benchmark catalog startup and a small window over a large source list, including files opened and bytes fetched, rather than only whole-mosaic throughput.

8. Pansharpening

No native pansharpen command appears in the inventory. Build on the existing resamplers and raster I/O. Start with GDAL-compatible weighted Brovey processing, explicit spectral weights, nodata, bit depth and extent handling. Add virtual output after the VRT extensions. GDAL produces either materialized rasters or pansharpened VRTs. GDAL pansharpen, pansharpened VRT algorithm.

Acceptance: compare against GDAL on real panchromatic/multispectral pairs, unequal resolutions, masked borders, explicit weights and integer output clipping.

9. Native multidimensional API

The native multidimensional API now exposes groups, named dimensions, coordinate arrays and lazy positive-stride slicing/transposition across numeric netCDF/HDF5 and Zarr arrays. Rust, CLI/MCP, Python and JavaScript bindings retain packed values and native integer precision. Existing raster views remain available. General string/compound arrays are listed as unsupported; netCDF snapshots currently materialize numeric arrays, while Zarr snapshots iterate chunks.

Conformance compares real NOAA wind, JPSS geolocation, Sentinel-5P and SRTM arrays with GDAL MD and xarray. A wind-derived four-dimensional cube checks time/depth slices, coordinate/chunk/unlimited-axis preservation across netCDF–Zarr snapshots, and reads that skip damaged unselected chunks. A counting-storage regression measures selected hyperslab I/O. The documented numeric scope follows GDAL's multidimensional model.

Suggested sequence

Start with overview-aware warping. Next implement geolocation warping and richer VRT support to complete existing reader workflows. SQL, viewshed and pansharpening can ship independently. Schedule indexed catalogs and the multidimensional model according to demand for continental mosaics and scientific data cubes, respectively.

Weighted/fractional zonal statistics now extend the existing zonal operation, with aligned weighting rasters and exactextract-style coverage fractions (see the CLI for the supported statistics and numeric allowances). Read-only OGC API Features access is tracked in #522, with supported scope and reference comparisons in the CLI and testing documentation. Overview refresh tooling remains a follow-up. Avoid prioritizing proprietary drivers and core database connections: the current roadmap explicitly excludes them and provides an optional GDAL backend for the long tail.

The initial recommendation-only assessment used source and documentation inspection. Fill-nodata implementation validation now includes Rust regression tests, Python comparisons with GDAL 3.12.2, real USGS terrain comparisons in the dem-calc conformance suite, and JavaScript tests of the WebAssembly build. The separate terrio.rio compatibility CLI still rejects fillnodata; the new operation is available through terrio.fill_nodata and terrio fill-nodata. Effort, demand and performance benefits for the remaining recommendations are estimates.

Cumulative viewshed uses optional observer coordinates in the CLI. Python and JavaScript retain the existing required positional coordinates; pass 0, 0 for ACCUM, which ignores them.