GDAL parity
rwarp is checked against GDAL output, component by component. The test
fixtures are generated by scripts in tests/fixtures/ from a stated GDAL
build, and the tests compare rwarp's results to them.
What has been checked
| Component | Test data | Result |
|---|---|---|
GenImgProjTransformer | Sentinel-2 B04 | 65536 of 65536 pixels identical to pixel_extract |
ApproxTransformer | same | bit-identical to GDALApproxTransform |
compute_source_window | GEBCO global, Fiji LCC | matches GDAL's source windows |
plan_chunks (with ChunkOptions::gdalwarp) | whole lon/lat sources to Mercator, south polar stereographic (pole in the target) and an overhanging LAEA across 180 degrees, -wm 100000 or the default, MEM and 32x32-tiled targets (16 cases, 1 to 28 chunks) | chunk list (target and source windows) identical; kernel scale per chunk within GDAL's logged precision |
Whole-raster chunked warp (warp_chunks) | same, near, bilinear, cubic, lanczos, nodata outside the source | identical to gdalwarp, every pixel, with both CRS backends |
Longitude wrapping (CENTER_LONG), GenImgProjTransformer::with_source_size | regional lon/lat sources across 180 degrees stored as 170..190, -190..-170 and 140..200, a global source in 0..360 and one 361 degrees wide under a -te of 90..270; to LAEA centred on 170 and 180, Mercator windows across 180 and whole-world, south polar stereographic and lon/lat targets (11 cases, near, bilinear, cubic; chunked and in one chunk) | centre longitude as GDAL's (%g); 162 target-to-source and source-to-target points per case through the exact transformer within 7e-13 px of GDAL's (libproj; proj4rs the same away from the poles, where it fails and PROJ returns a point); chunk lists and every output pixel identical to gdalwarp with both CRS backends. Without the wrap, 17% to 61% of each target differs |
| Retry of points just off the source window | same (LAEA, near) | 2 pixels per 10000 differ without it |
| Nearest neighbour | Sentinel-2, IBCSO Mawson | 65536 of 65536 identical to gdalwarp -r near |
| Bilinear, cubic, lanczos | GEBCO at about 1:1 | 65536 of 65536 identical |
| Kernel widening, scale estimate | synthetic, 1.5:1 to 4:1, anisotropic, nodata | within 1e-6 relative of gdalwarp |
| Point kernels at any ratio | same, against -wo XSCALE=1 -wo YSCALE=1 | within 1e-6 relative |
| Average, GDAL box | plus upsampling, overhanging and 30 degree rotated targets | within 1e-6 relative of -r average |
| Average, exact footprint | axis-aligned cases against the box; a 45 degree diamond against the analytic mean | within 1e-5 |
Mode, all MODE_TIES | categorical Byte/Int16/Int32/Float32, 4:1 and 2.78:1, nodata, rotation | identical to -r mode (26 cases) |
| Nodata, per band and unified, average thresholds | single- and three-band Float32, every method | within 1e-6 relative (30 cases) |
| Rounding, clamping, nodata avoidance, source bounds | Byte/Int16/Int32/Float32 | identical (12 cases) |
| Antimeridian split read | GEBCO, Fiji tile | 70 times fewer source pixels than a naive read |
| Geolocation arrays | synthetic swath, both conventions, curvilinear source and target, both inverse methods | coordinates within 1e-6 px; nearest warps identical to gdalwarp -et 0 |
Geolocation backmap inverse (GeolocInverse::Backmap) | 9 synthetic cases: a curved swath by cell centres and by corners to lon/lat (fine and coarse) and south polar stereographic, a polar ring and a swath across the antimeridian, holes in the arrays with source nodata, overlapping scans; -et 0, and -et 0.125 for two | success identical to GDAL's transformer at every target pixel centre (31,356), positions within 5e-10 px; the backmap's estimate (GDAL_GEOLOC_USE_MAX_ACCURACY=NO) the same; 24 warps (near, bilinear, cubic) identical to gdalwarp, every pixel, with proj4rs and the fixtures' PROJ; with an older libproj the nearest warps are identical and interpolated values within 2e-8 relative |
Neighbour on the edge of the transformer's domain (GWKAdjustSrcOffsetOnEdge) | the holes case above | 12 of 4560 pixels differ without it |
| Cutline: centre rule, all touched, blend distance | synthetic Float32; edges and vertices through pixel centres, holes, multipolygons, source nodata, every method but mode, 1:1, 2.78:1 and upsampled, 30 degree rotated, LAEA to lon/lat with the cutline in either CRS | within 1e-6 relative of gdalwarp (67 cases; the reprojected ones identical with libproj and proj4rs) |
| Cutline tie rules | GDAL's own rasterizer on centre-tie polygons | identical, including centres exactly on diagonal edges |
The fixtures use GDAL 3.14.0dev (739fcf03bf). The first geolocation
fixture was made with GDAL 3.8.4; regenerated with 3.14 only the CRS
step's last digits moved (1e-8 px), and GDAL's geolocation code is the
same from 3.12 to 3.14.
Longitude wrapping (CENTER_LONG)
For a source CRS that is geographic in degrees and has a geotransform,
GDALCreateGenImgProjTransformer2 (and so gdalwarp) wraps source
longitudes around the source raster's centre: InsertCenterLong takes the
middle of the smallest and largest x of the raster's four corners, writes
it with %g (six significant digits: a raster centred on 147.1284567 is
wrapped at 147.128), and gives it to the coordinate transformation as its
source centre. OGRProjCT then adds 360 to a source longitude below
centre - 180 and takes 360 from one above centre + 180, once, on the
way into the forward transform and on the way out of the reverse one that
the warper uses. A raster wider than 360 degrees gets no centre, unless the
target is geographic too and gdalwarp's -te crosses 180 degrees. gdalwarp
turns it off (INSERT_CENTER_LONG=FALSE) only to bring a cutline into
source pixels.
rwarp does the same in GenImgProjTransformer once it knows the source
raster's size (with_source_size; with_target_extent for -te).
CenterLong::Off and CenterLong::At choose otherwise, and
source_center_long gives GDAL's value for a transformer built with
with_backend, which does not know its CRSs. Whether a CRS is geographic
in degrees is decided as GDAL does: with libproj, a bound CRS by its base
and a compound CRS by its horizontal part; NTF (Paris), in grads, is not.
Where rwarp deliberately differs
| Behaviour | GDAL | rwarp | Why |
|---|---|---|---|
| Kernel scale estimate in a chunked warp | taken at the top-left of the whole output for every chunk | taken over the window being warped (ScaleOrigin::Chunk, the default; ScaleOrigin::Output or ChunkOptions::gdalwarp for GDAL's) | GDAL's sample points omit the chunk offset; chunked and single-chunk gdalwarp output differ. On the Mercator lanczos fixture, half the pixels differ |
Source window padding under XSCALE/YSCALE | padded for the estimated scale | padded for the smaller of the given and estimated scales (ChunkOptions::pad_for_fixed_scale, default on; off in ChunkOptions::gdalwarp) | a kernel stretched by the given scale is otherwise clipped at the window edge, so tiles warped on their own differ there |
MODE_TIES for 32/64-bit and float data | Min/Max skipped for a value's first occurrence | the 8/16-bit behaviour for every type | the same data should not give different modes by type |
| Geolocation inverse, by default | backmap plus local search | index of all quadrilaterals (GeolocInverse::QuadIndex); GDAL's is GeolocInverse::Backmap, identical to gdalwarp | the index does not depend on a sampled raster. Same answer inside the raster on grids without overlapping quadrilaterals; differs at folds and at the seam of a grid wrapping 360 degrees, where GDAL leaves pixels empty; see Curvilinear grids |
| Cutline all touched | centre rule plus each boundary line walked pixel to pixel, with tolerances | pixels with any covered area (above 1e-9) | GDAL also burns pixels touched only at a corner (below-right of it), skips edges within 1e-4 of a pixel wall, and burns one column for an edge within 0.01 px of vertical even where the polygon is in the column to its left |
Cutline::from_transformer densification | skipped when -cblend is given (unless GDALWARP_DENSIFY_CUTLINE=YES) | always | the helper does not know about blending; densifying only makes the transformed cutline more accurate |
| Longitude wrapping of a geographic source | on whenever the source has a geotransform (the dataset gives its size) | on once GenImgProjTransformer::with_source_size is given | GenImgProjTransformer::new takes no source size; a transformer built as in 0.3.0 wraps nothing, as before |
Not modelled
- The
CENTER_LONGconfiguration option (--config CENTER_LONG), which GDAL applies to the target CRS as well as the source;CenterLong::Atwraps the source only, as aCENTER_LONGin the source CRS does. CHECK_WITH_INVERT_PROJ,STREAMABLE_OUTPUT, and GDAL's shortcut for a translation by whole pixels in chunk planning.- Geolocation:
SWAP_XY, GDAL's experimental quadtree inverse (GDAL_GEOLOC_INVERSE_METHOD=QUADTREE), and its temporary-file backmap for arrays of 24 megapixels or more (same arithmetic as the in-memory one rwarp ports). GDAL's alpha and cutline-density variant of the domain-edge neighbour (GWKAdjustSrcOffsetOnEdgeUnifiedSrcDensity).
GDAL's retry of a sample point just off the source window with the exact
transformer is done: warp() transforms such a point again on its own,
which through ApproxTransformer (exact for five points or fewer, as
GDAL's) is exact, so the kernel stays independent of transformer type.
GDAL issues found along the way
Matching gdalwarp closely turned up behaviour in GDAL itself that looks
unintended: the chunk-offset scale estimate, MODE_TIES by data type, the
MODE_TIES=FIRST documentation, edge weights above 1 in average when a
target pixel overhangs the source, destination nodata avoidance with a float
working type, and two different rounding rules into integer outputs. They
are written up with repros in
mdsumner/osgeo-upstream.
The CENTER_LONG work found that GDAL computes its special points (the
poles, which force grid sampling of a chunk that contains one) without
reprojecting them when the source CRS is geographic:
GDALTransformLonLatToDestGenImgProjTransformer takes (0, +-89.9999) as
target coordinates, so for a projected target they land at its false
origin. On the chunk fixtures this changes no chunk list.
The cutline work found more, not yet in that write-up: all-touched burns
pixels the boundary touches only at a corner, and a whole column along an
edge within 0.01 pixels of vertical; blend distance dilates the cutline for the area methods instead
of feathering it; and -cblend turns off cutline densification.
The geolocation backmap port found these: the inverse's spiral search
admits a square when iX >= -1 || iY >= -1 (probably meant &&), so it
reaches extrapolated squares well beyond the top and left edges but only
one beyond the right and bottom; on arrays wrapping 360 degrees the
backmap averages positions from both ends of the rows across the
antimeridian, and the search radius of 3 cannot recover, leaving a line of
empty pixels along the antimeridian (51 of 4900 in the polar test case);
next to a missing sample PixelLineToXY extrapolates past the last line
along the line only, folding the extrapolated square back over the grid;
isPointOnRingBoundary tests for an exactly zero cross product (an OGR
FIXME), so only vertices and axis-aligned edges count as on the boundary;
GenerateBackMap casts a negative sample position to unsigned; and its
per-block sample ranges drop the last samples beyond the arrays' edge
when the edge falls within two samples of a multiple of 256 (they land in
a block that is never processed). In the warp kernel,
GWKAdjustSrcOffsetOnEdge gives a source pixel that is nodata the value
of a neighbour whenever a corner near it fails to transform, so source
nodata next to holes in geolocation arrays is filled.