rwarp

rwarp is GDAL's warp pipeline in Rust: coordinate transforms, approximate transforms, source window planning and resampling kernels, with no runtime GDAL dependency. It follows GDAL's algorithms closely, written with reference to the GDAL source and checked against gdalwarp output, but it is an independent implementation, not affiliated with or endorsed by the GDAL project.

Try it in the browser: the demo warps web map tiles into a projection of your choice, client-side, compiled to WebAssembly.

Why another warper

GDAL's warper bundles planning, reading, coordinate mapping and resampling into one operation. rwarp separates them, so you can see and control what will be read before any pixel is touched:

target grid
    |
    v
[plan]  plan_chunks / compute_source_window
    |     which chunks, and which source pixels does each need?
    |     inspectable, serialisable, resumable
    v
[read]  the caller fetches bytes: COG tiles, Zarr chunks, byte-reference stores
    |
    v
[warp]  warp(): nearest, bilinear, cubic, lanczos, average, mode

That makes it a building block for caches, tile servers, browser viewers and planners that need to know the I/O cost of a reprojection up front.

What it covers

CRS transforms go through a small trait with two backends: libproj (the default, full PROJ) or proj4rs (pure Rust, builds for wasm32-unknown-unknown). See Getting started.

Status

Early development: the API will change. Each component is validated against GDAL output; the GDAL parity page lists what has been checked, how, and where rwarp deliberately differs. Feedback and testing are very welcome, at github.com/hypertidy/rwarp.

Getting started

Add the crate

[dependencies]
rwarp = "0.3"

For unreleased changes on main, depend on the repository:

[dependencies]
rwarp = { git = "https://github.com/hypertidy/rwarp" }

Choose a CRS backend

cargo test                                            # libproj
cargo test --no-default-features --features proj4rs   # pure Rust
cargo build --no-default-features --features proj4rs --target wasm32-unknown-unknown

Or bring your own: GenImgProjTransformer::with_backend takes any pair of CrsTransform implementations (one each way), for example a JavaScript PROJ build reached from WebAssembly.

A first warp

This plans and warps one 256 x 256 tile in a local Lambert azimuthal equal-area projection from a Web Mercator tile pyramid. It is examples/tile_laea.rs; run it with cargo run --example tile_laea.

use rwarp::source_window::compute_source_window;
use rwarp::transform::GenImgProjTransformer;
use rwarp::warp::{warp, Nodata, ResampleAlg, SourceBuffer, WarpOptions};

const WEBMERC_HALF: f64 = 20037508.342789244;
const LAEA_TAS: &str = "+proj=laea +lat_0=-42 +lon_0=147 +datum=WGS84 +units=m +no_defs";

fn main() -> Result<(), String> {
    // Source: Web Mercator at zoom 8, as one virtual raster.
    let n = 256 * (1 << 8);
    let px = 2.0 * WEBMERC_HALF / n as f64;
    let src_gt = [-WEBMERC_HALF, px, 0.0, WEBMERC_HALF, 0.0, -px];
    let src_size = [n, n];

    // Target: a 256 x 256 tile, 400 km across, centred on Tasmania.
    let extent = 400_000.0;
    let dst_gt = [-extent / 2.0, extent / 256.0, 0.0, extent / 2.0, 0.0, -extent / 256.0];

    // 1. Transformer: target pixel -> source pixel, approximated as gdalwarp does.
    let exact = GenImgProjTransformer::new("EPSG:3857", src_gt, LAEA_TAS, dst_gt)?;
    let t = ApproxTransformer::new(exact, 0.125);

    // 2. Plan: which source pixels does the tile need? No pixels are read.
    let win = compute_source_window(&t, [0, 0], [256, 256], src_size, 1)
        .ok_or("tile does not map onto the source")?;
    println!(
        "source window: offset ({}, {}), size {} x {} pixels",
        win.xoff, win.yoff, win.xsize, win.ysize
    );

    // 3. Read: fetch that window (here, synthesise it).
    let (w, h) = (win.xsize as usize, win.ysize as usize);
    let data: Vec<f32> = (0..w * h).map(|i| (win.xoff as usize + i % w) as f32).collect();
    let src = SourceBuffer {
        data: &data,
        ncol: w,
        nrow: h,
        nbands: 1,
        col_off: win.xoff as usize,
        row_off: win.yoff as usize,
    };

    // 4. Warp.
    let opts = WarpOptions::new(ResampleAlg::Bilinear)
        .with_src_nodata(Nodata::Value(f64::NAN));
    let out = warp(&t, &src, 256, 256, &opts);
    let valid = out.iter().filter(|v| !v.is_nan()).count();
    println!("warped {} of {} target pixels", valid, out.len());
    Ok(())
}

The four steps are the whole pipeline:

1. Transformer. Maps target pixels to source pixels. Wrapping it in ApproxTransformer with a 0.125 pixel error, as gdalwarp does by default, transforms a few points per scanline exactly and interpolates the rest.
2. Plan. compute_source_window finds the source pixels the tile needs, without reading any.
3. Read. Fetch that window however your data is stored.
4. Warp. warp resamples into the target, with options for the method, kernel scaling, nodata and more (see Resampling).

Transformers

A transformer maps coordinates between target pixels and source pixels. Every part of rwarp that needs coordinates (planning, kernels) takes anything that implements the Transformer trait, so they compose freely.

pub trait Transformer {
    /// dst_to_src = true: input is target pixel/line, output is source pixel/line.
    fn transform(&self, dst_to_src: bool, x: &mut [f64], y: &mut [f64]) -> Vec<bool>;
}

It follows GDAL's transformer callback: coordinate arrays are modified in place, and one flag per point says whether it transformed. The warp kernels always call it target to source.

The transformers

Geotransforms are in GDAL's order, [x0, dx, rx, y0, ry, dy]. Pixel coordinates are continuous: pixel i covers [i, i + 1), so its centre is i + 0.5.

Sources across 180 degrees

A lon/lat source stored as 170..190, as -190..-170, or as a 0..360 global grid has longitudes the CRS step never returns: inverse projections give -180..180. gdalwarp handles this by wrapping the source's longitudes around its centre (GDAL's CENTER_LONG), and GenImgProjTransformer does the same once it knows the source raster's size:

use rwarp::transform::{CenterLong, GenImgProjTransformer};

let t = GenImgProjTransformer::new("EPSG:4326", src_gt, dst_crs, dst_gt)?
    .with_source_size([src_width, src_height]);
// the middle of the source's corner longitudes, to six digits as GDAL
let centre = t.center_long();

Target to source, a longitude more than 180 degrees from the centre is brought back by 360 (once), so a target crossing 180 finds both halves of the source; source to target, the same is done before the CRS step. Only a source CRS that is geographic in degrees is wrapped, as in GDAL. A source wider than 360 degrees is not, unless the target is lon/lat too and its extent crosses 180 (give it with with_target_extent, gdalwarp's -te).

with_center_long(CenterLong::Off) turns it off (gdalwarp does, for cutlines) and CenterLong::At(lon) sets the centre. A transformer from with_backend does not know its CRSs, so CenterLong::Auto does nothing there; use CenterLong::At(source_center_long(&src_gt, size, None).unwrap()) for GDAL's value. Without with_source_size, nothing is wrapped, as in rwarp 0.3.0.

CRS backends

The CRS step goes through the CrsTransform trait:

pub trait CrsTransform {
    fn convert(&self, x: f64, y: f64) -> Option<(f64, f64)>;
}

Its direction is fixed at construction. Geographic coordinates are always longitude, latitude in degrees, whichever backend is used. crs_pair builds the two directions with the backend compiled in; with_backend constructors accept your own.

Vector quantities

For velocity fields and other (u, v) quantities stored in grid coordinates, reprojecting the scalar values is not enough: the vectors also rotate with the local frame. Transformer::jacobian gives the local 2 x 2 derivative of a transform at a point (by finite differences unless an implementation overrides it), the piece a vector-preserving warp needs.

Planning reads

Planning answers "which source pixels does this target need?" before any pixels are read. It needs only a transformer and the source raster's size.

One window

let win = compute_source_window(&t, dst_off, dst_size, src_raster_size, resample_padding)?;
// win.xoff, win.yoff, win.xsize, win.ysize: the source window in pixels
// win.fill_ratio: how much of the unclamped window lies inside the source

This follows GDAL's ComputeSourceWindow: it samples the target window's edges (and, if that fails, a grid over its interior), transforms them to the source, and takes their bounding box, padded for the resampling kernel (0 for nearest, 1 for bilinear, 2 for cubic, 3 for lanczos). It returns None if too few points transform.

A low fill_ratio means much of the window falls outside the source, or the target straddles a discontinuity such as the antimeridian.

compute_source_window is a quick estimate that pads one pixel beyond the kernel radius and skips GDAL's forward probe and special points. For the window exactly as gdalwarp computes it, use ChunkPlanner::new(&t, src_size, dst_size, opts).source_window(dst_off, dst_size).

Splitting into chunks

use rwarp::plan::{plan_chunks, ChunkOptions};

let opts = ChunkOptions::gdalwarp(ResampleAlg::Bilinear)
    .with_memory_limit(64.0 * 1024.0 * 1024.0)            // gdalwarp -wm
    .with_pixel_cost(DataType::Float32, 1, true, true);   // with nodata
let plan = plan_chunks(&approx, src_size, dst_size, &opts);
for c in plan.live() {
    // c.dst_off, c.dst_size: a piece of the target
    // c.src_off, c.src_size: what it reads from the source
    // c.fill_ratio, c.scale: GDAL's fill ratio and kernel scale
    let tile = warp(&approx, &read(c.src_off, c.src_size), c.dst_size[0], c.dst_size[1],
                    &c.warp_options(&WarpOptions::new(ResampleAlg::Bilinear)));
}

This is gdalwarp's planning (CollectChunkList): the source window of the whole target is computed, and while it would take more working memory than the limit, or its fill ratio is below 0.5 (much of it outside the source), the target is halved along its longer side and each half planned again. Each chunk is then warped on its own, with the kernel scale GDAL would estimate for it. Chunks are sorted top to bottom, then left to right, as GDAL warps them.

The chunk boundaries matter to the output, not only to memory: they set the kernel scale, the edges where kernels are clipped, and the scanlines the approximate transformer interpolates over. With ChunkOptions::gdalwarp and the same memory limit, the plan and the output equal gdalwarp's for a whole raster (see GDAL parity). ChunkOptions::new differs in two documented places: each chunk's kernel scale is estimated over the chunk (GDAL estimates every chunk at the output's top-left), and source windows are padded for a fixed kernel scale.

A few inputs need the caller: GDAL's special points (longitude 0 at latitude +-89.9999, in target pixels, which make every chunk sample a grid when one falls in the target) need the CRSs, and the target's block size matters only for a tiled target (OPTIMIZE_SIZE).

warp_chunks runs a whole plan from a source buffer in memory, which is how the parity tests use it.

At the antimeridian

with_split_antimeridian(true) adds the split of 0.2's collect_chunk_list: a chunk whose source window spans the whole source width (as at 180 degrees) is split where source x jumps, so a tile straddling 180 degrees becomes tiles each side with compact source windows rather than one spanning the globe (a 70 times reduction in pixels read for a GEBCO tile over Fiji). GDAL reaches a similar result only through its memory limit.

Why planning is separate

Because the plan is just data (target windows and source windows), it can be inspected, costed, cached, serialised (the serde feature) and resumed. You can count the COG tiles or Zarr chunks a reprojection touches, pick an overview level from the scale before reading, or batch requests across many targets, all before any I/O.

Resampling

warp resamples a source window into a target, for any sample type (u8, i8, u16, i16, u32, i32, f32, f64) and any number of pixel-interleaved bands, in one pass:

let out = warp(&transformer, &SourceBuffer { .. }, dst_ncol, dst_nrow, &opts);

SourceBuffer is the window you read, laid out (row, col, band), with its offset in the full source image. WarpOptions holds everything else; WarpOptions::default() is nearest neighbour with GDAL's defaults. Build it with WarpOptions::new(alg) and the with_ methods, one per option:

let opts = WarpOptions::new(ResampleAlg::Average)
    .with_footprint(Footprint::Exact { densify: 1 })
    .with_src_nodata(Nodata::Value(-9999.0));

WarpOptions and the option enums are #[non_exhaustive], so new options can arrive in minor releases. Outside rwarp this means no struct literals for WarpOptions, GeolocOptions or AverageThresholds, and a _ arm when matching on the enums.

Kernel widening when downsampling

When the output is coarser than the source, GDAL widens the bilinear, cubic and lanczos kernels to cover each output pixel's source area. rwarp does the same by default:

• WidenWhenDownsampling estimates the target/source scale over the target window, as GDAL does per warp chunk, then widens the kernel support to ceil(radius / scale) on any axis where the scale is below 1. Upsampling never narrows the kernel.
• Fixed { x, y } gives the scale, like GDAL's XSCALE and YSCALE. Use it to make results independent of how the target is tiled.
• Point keeps the fixed 2 x 2, 4 x 4 and 6 x 6 support at every ratio, which is gdalwarp -wo XSCALE=1 -wo YSCALE=1.

As in GDAL, bilinear and cubic widen only once the scale reaches 0.5 or less; lanczos widens at any scale below 1.

One deliberate difference: rwarp estimates the scale over the window actually being warped. GDAL takes its estimate at the top-left of the whole output for every chunk, so a chunked gdalwarp can widen by the wrong amount away from the origin (see GDAL parity).

Average and mode

Both weigh the source pixels a target pixel's footprint covers, measuring area in source pixel space:

• Footprint::GdalBox is GDAL's footprint: the axis-aligned box between the transformed top-left and bottom-right corners of the target pixel, with fractional weights at its edges. It matches gdalwarp, quirks included; under rotation or strong curvature the box can be far from the true footprint.
• Footprint::Exact { densify } transforms the target pixel's outline and intersects it exactly with the source pixels. With renormalisation over valid pixels, Average is then first-order conservative remapping. It costs about 25 times the box. The same overlaps are available as a weight table; see Conservative weights.

Mode counts covered area per value and breaks ties with ModeTies. Counts are re-compared every time a value's count grows, in scan order, so First is the tied value that reached the winning count first. GDAL does this for 8- and 16-bit data; rwarp does it for every type.

Output values

• Integer outputs are clamped to the type and rounded half up, floor(v + 0.5), as GDAL does with an integer working type. Nearest and mode copy source values and never round.
• Nodata avoidance: a resampled value equal to the nodata value is moved off it by one step, so a valid pixel never reads back as missing.
• Source bounds: sample points outside [0, n) of the source are off the source for every method except average and mode, which use footprints.

Conservative weights

Footprint::Exact measures how much of each source pixel every target pixel's outline covers, then averages and throws the measurements away. rwarp::weights keeps them: a sparse table with one entry per overlapping (target cell, source cell) pair. That table is the weight matrix of first-order conservative remapping, ESMF's conserve method. Compute it once, then apply it to every band, time step or variable on the same grids.

use rwarp::weights::{exact_weights, ApplyOptions, GridWindow, Normalization, WeightOptions};

// The source grid needs no data, only its size (and offset, for a window).
let src = GridWindow::new(src_ncol, src_nrow);
let w = exact_weights(&transformer, &src, dst_ncol, dst_nrow, &WeightOptions::new(4));

// Apply to any pixel-interleaved buffer on that source grid.
let opts = ApplyOptions::new().with_src_nodata(Nodata::Value(-9999.0));
let out: Vec<f32> = w.apply(&source_buffer, Normalization::FracArea, &opts);

The transformer is the same one warp() takes (target pixel to source pixel), so anything rwarp can warp it can make weights for.

What's in a Weights

Areas are measured in source pixel space, where a source cell has area 1. For a source on a regular grid in its own CRS that is area in the source CRS, up to a constant. Target cells have edges straight in target pixel space, each split into densify segments whose end points go through the transformer exactly.

frac is relative to the source cell, not the target, so it doesn't depend on how the target is normalised. dst_area gives either normalisation:

• Normalization::FracArea: S = frac / (sum of frac over the target cell), the area-weighted mean over the part of the cell that has source. This is what warp() computes. In ESMF it's fracarea; in xESMF, conserve_normed.
• Normalization::DestArea: S = frac / dst_area. A cell half off the source gets half the weight, so sum(value * dst_area) over the target equals the integral of the source over the region the target covers. ESMF destarea, xESMF conserve.

Weights::normalised(norm) returns S per entry. Weights::apply also skips source nodata: FracArea renormalises over the valid entries, and DestArea doesn't (a missing source cell counts as zero).

Conservation

Neighbouring target cells build their outlines from the same transformed points, so together they tile source pixel space with no gaps or overlaps. For every source cell the target grid covers completely, its frac sums to 1 over the target cells. Where the grid's edge cuts a source cell, the sum is the covered fraction. This needs no reference to check, and tests/weights.rs checks it against the coverage of the target grid's outer boundary:

Densify

Both settings conserve: conservation comes from the shared edges, not from how straight they are. They define slightly different cells. With densify 1 a target cell is the quadrilateral through its four transformed corners in source pixel space; with more, its edges follow the curve the target's straight edges map to. Under an affine transform the two agree to rounding (4e-14). For lon/lat to LAEA at 20 km, frac differs by up to 7e-5 and dst_area by up to 6e-7 relative between densify 1 and 4. Cost grows with densify (the transform dominates): 4 is about 4 times 1 when target cells are smaller than source cells.

Same numbers as warp()

warp() with Footprint::Exact and ResampleAlg::Average uses the same lattice and the same per-pixel coverage step. For compatibility with rwarp 0.2 it snaps coverage within 1e-6 of 1 to 1, drops coverage at or below 1e-6 and rounds the rest to f32. exact_weights keeps every positive overlap, unrounded, by default. With WeightOptions::warp_compatible(densify) it does what warp() does, and apply(FracArea) then gives warp()'s output bit for bit (tested on the exact-footprint kernel fixtures and on lon/lat to LAEA with two bands, nodata, NaN, every BandNodata rule, a mask band and Int16 to Byte).

warp() reduces directly and never stores weights. Use it for a one-off. Use weights when the same grids are applied many times, when you need DestArea, or when the weights are the product.

Tiling and streaming

WeightOptions::with_dst_offset(col, row) computes the weights of a target window, with full-grid indices. A window's lattice is the whole grid's, so tiles give exactly the whole grid's entries: row bands concatenate, and column offsets agree too when densify is a power of two.

exact_weights_with(.., |cell| ..) calls back once per target cell with its entries, storing nothing, for problems too big to hold. Entries cost 24 bytes each (four u32 and an f64). Over a 0.25° global source:

(densify 1, release build, one thread, PROJ backend.)

Weight files

rwarp has no I/O. Writing ESMF/SCRIP netCDF, xESMF's format or Parquet is wxylib's job, and the columns map directly. With the full grids' widths:

Weights::esmf_row_col(dst_ncol_full, src_ncol_full) returns the two index columns.

Limits

• Planar areas in source pixel space. ESMF measures on the sphere with great-circle edges. A spherical measure is planned (design/conservative.md, stage 3).
• No antimeridian handling. A target cell whose outline crosses the source's wrap column folds across the whole source width, and its weights are wrong. Keep target grids clear of the source's seam for now, or use a source that is continuous across the target.
• Destination cells are what a transformer describes. Curvilinear grids work through GenImgGeolocTransformer with its bilinear cell definition. Cells from corner bounds and meshes are later stages.

Curvilinear grids

Many datasets have no geotransform: satellite swaths, ocean models on rotated or tripolar grids, and anything else whose cell positions are given by 2D longitude and latitude arrays. rwarp handles these with geolocation arrays, following GDAL's geolocation transformer (alg/gdalgeoloc.cpp).

The model

Each geolocation sample is a point. Between four neighbouring samples, positions are interpolated bilinearly; beyond the outermost samples, they are extrapolated from the closest grid square. Going the other way, from a coordinate to a pixel, rwarp finds the quadrilateral of samples that contains the point and solves the exact inverse bilinear interpolation (two ways to find it, below).

Where a sample sits within its pixel is set by GeolocConvention:

• PixelCenter: samples are cell centres, as for CF coordinate variables. Cell edges then fall at the midpoint between neighbouring centres (in the arrays' own coordinates), the outer edges half a spacing beyond the outermost centres, and each corner at the mean of four centres.
• TopLeftCorner: samples are top-left pixel corners. This is GDAL's default when GEOREFERENCING_CONVENTION is not set, so centre coordinates read through the default land half a cell off. Say which one you have.

GeolocOptions also carries pixel and line offsets and steps (for arrays that sample every Nth pixel), a nodata value for the X array, and whether longitudes wrap at the antimeridian.

Warping from a swath

GenImgGeolocTransformer puts geolocation arrays on either side of the usual CRS step. This is examples/curvilinear_swath.rs (cargo run --example curvilinear_swath):

use rwarp::geoloc_transform::{
    GenImgGeolocTransformer, GeolocArrays, GeolocConvention, GeolocOptions, PixelGeo,
};
use rwarp::warp::{warp, Nodata, ResampleAlg, SourceBuffer, WarpOptions};

fn main() -> Result<(), String> {
    let (ncol, nrow) = (24, 18);
    let mut lon = Vec::with_capacity(ncol * nrow);
    let mut lat = Vec::with_capacity(ncol * nrow);
    for j in 0..nrow {
        for i in 0..ncol {
            let (i, j) = (i as f64, j as f64);
            lon.push(142.0 + 0.30 * i + 0.08 * j + 0.004 * j * j);
            lat.push(-42.0 - 0.22 * j + 0.05 * i - 0.002 * i * i);
        }
    }

    // Coordinates are cell centres, as CF coordinate variables are. GDAL's
    // default convention is TopLeftCorner, so say so explicitly.
    let opts = GeolocOptions::default()
        .with_convention(GeolocConvention::PixelCenter)
        .with_lon_wrap(true);
    let swath = GeolocArrays::new(lon, lat, ncol, nrow, opts)?;

    // Target: a 40 x 32 grid of 10 km pixels in a local LAEA.
    let dst_crs = "+proj=laea +lat_0=-45 +lon_0=146 +datum=WGS84 +units=m +no_defs";
    let dst_gt = [-420000.0, 10000.0, 0.0, 330000.0, 0.0, -10000.0];
    let t = GenImgGeolocTransformer::new(
        PixelGeo::Geoloc(swath),
        "+proj=longlat +datum=WGS84 +no_defs",
        PixelGeo::affine(dst_gt)?,
        dst_crs,
    )?;

    let data: Vec<i32> = (0..nrow)
        .flat_map(|r| (0..ncol).map(move |c| (c + 1000 * r) as i32))
        .collect();
    let src = SourceBuffer { data: &data, ncol, nrow, nbands: 1, col_off: 0, row_off: 0 };
    let wopts = WarpOptions::new(ResampleAlg::NearestNeighbour)
        .with_src_nodata(Nodata::Value(-1.0));
    let out = warp(&t, &src, 40, 32, &wopts);

    let covered = out.iter().filter(|&&v| v != -1).count();
    println!("{covered} of {} target pixels fall on the swath", out.len());
    // The centre row, as source (col, row) pairs.
    let row: Vec<String> = out[16 * 40..17 * 40]
        .iter()
        .map(|&v| if v < 0 { ".".to_string() } else { format!("{},{}", v % 1000, v / 1000) })
        .collect();
    println!("centre row: {}", row.join(" "));
    Ok(())
}

823 of the 1280 target pixels fall on the swath, the same count gdalwarp gives for this grid.

Warping onto a curvilinear grid

The same transformer works the other way round: put the geolocation arrays on the target side to resample a regular grid onto a swath or model grid. This is GDAL's DST_METHOD=GEOLOC_ARRAY, and rwarp matches gdalwarp there too.

Two inverse methods

From a coordinate to a pixel, both methods find a quadrilateral of samples that contains the point and solve the same exact inverse bilinear interpolation in it. They differ in how they find it, chosen with GeolocOptions::with_inverse:

• GeolocInverse::QuadIndex (the default) indexes every quadrilateral, plus one ring of extrapolated ones outside the arrays, and takes the first that contains the point in row-major order. Its domain is the union of the quadrilaterals.
• GeolocInverse::Backmap is GDAL's method (GEOLOC_INVERSE_METHOD BACKMAP, GDAL's default): a float32 raster over the arrays' extent, filled by forward-projecting the arrays and patching holes as GDAL does, gives an approximate pixel, and a spiral search of up to 3 squares around it (5 near the poles) finds the quadrilateral. It is ported in GDAL's arithmetic order, so warps through it are identical to gdalwarp's, pixel for pixel. GeolocInverse::BackmapEstimate stops at the backmap's estimate (GDAL_GEOLOC_USE_MAX_ACCURACY=NO). The backmap's oversampling is with_backmap_oversample (GEOLOC_BACKMAP_OVERSAMPLE_FACTOR, 1.3).

let opts = GeolocOptions::default()
    .with_convention(GeolocConvention::PixelCenter)
    .with_lon_wrap(true) // what GDAL does for lon/lat in [-180, 180]
    .with_inverse(GeolocInverse::Backmap);

Use Backmap to reproduce gdalwarp, and QuadIndex where the answer should not depend on a sampled raster. On a grid whose quadrilaterals do not overlap, the two agree wherever both succeed inside the raster. On the test cases (tests/geoloc_backmap.rs, 2600 to 4900 target pixels each, nearest neighbour):

Outside the raster GDAL's search also succeeds at more points (squares further than one beyond the top and left edges, within its search radius), which no kernel reads.

Other differences from GDAL

• rwarp treats NaN samples as missing, as well as the X array's nodata value; GDAL compares with the nodata value only.
• Whether longitudes wrap at the antimeridian is the explicit with_lon_wrap. GDAL decides it from the SRS (geographic, every longitude in [-180, 180]) or forces it with GEOLOC_NORMALIZE_LONGITUDE_MINUS_180_PLUS_180.
• SWAP_XY and GDAL's experimental quadtree inverse (GDAL_GEOLOC_INVERSE_METHOD=QUADTREE) are not implemented. GDAL keeps the backmap of arrays of 24 megapixels or more in temporary files, with the same arithmetic; rwarp keeps it in memory.

Cutlines

A cutline limits a warp to the source pixels inside a polygon, as gdalwarp -cutline does. rwarp takes the polygon in source pixel/line coordinates of the full source image, which is where gdalwarp puts it before warping:

use rwarp::cutline::{Cutline, CutlineRule};

// Polygons: an exterior ring and its holes, each ring a list of [x, y].
let cut = Cutline::new(vec![vec![exterior, hole]]);
let opts = WarpOptions::new(ResampleAlg::Bilinear).with_cutline(cut);
let out = warp(&transformer, &source_buffer, dst_ncol, dst_nrow, &opts);

Several polygons are a multipolygon (their union).

From world coordinates

• Cutline::from_geotransform(polygons, src_gt): the polygons are in the source raster's CRS. The map to pixels is affine, so nothing is densified, as in gdalwarp.

• Cutline::from_transformer(polygons, &transformer): the polygons are in whatever a Transformer maps to source pixels. For a cutline in another CRS, build a GenImgProjTransformer from that CRS to the source with the identity geotransform [0, 1, 0, 0, 0, 1] on the cutline side, so its "pixels" are CRS coordinates:

  let t = GenImgProjTransformer::new(src_crs, src_gt, "EPSG:4326", [0.0, 1.0, 0.0, 0.0, 0.0, 1.0])?;
  let cut = Cutline::from_transformer(lonlat_polygons, &t)?;

  Edges are densified as gdalwarp densifies them, so that each transformed edge is about one source pixel long.

Rules

with_blend_dist(d) is gdalwarp's -cblend d: within d source pixels of the boundary the density ramps from 1 inside to 0 outside, 0.5 on the line.

Centre ties. A pixel centre exactly on the boundary follows GDAL: on a left edge it is out, on a right edge in, on a top edge in, on a bottom edge out, and a hole whose top edge runs along a line of centres does not cut that line. Ties are decided exactly, and match GDAL on diagonal edges through centres too.

All touched counts a pixel the polygon covers by more than 1e-9 of its area. GDAL walks each boundary line instead, which differs in three places: GDAL also burns a pixel the boundary only touches at its corner (below-right of the corner), skips an edge within 1e-4 of a pixel wall, and burns a single column for an edge less than 0.01 pixels from vertical.

What a cutline does to each method

The cutline becomes a density per source pixel, computed over the source buffer only.

• Point methods (nearest, bilinear, cubic, cubic spline, lanczos): a target pixel whose sample point lands outside the cutline gets no value; pixels outside are dropped from every kernel, like nodata, so cubic falls back to bilinear near the boundary.
• Area methods (average, mode, min, max, quantiles, RMS): pixels outside the cutline are left out, as nodata pixels are.
• Blend distance: point methods interpolate the density along with the value. With a destination nodata the value is kept wherever the density is at least 0.0001; without one it is multiplied by the density, fading to the zero fill, which is what gdalwarp does when warping into a new file. Area methods only use the density as validity, so for them the blend distance widens the cutline.
• Exact: area methods weight each pixel by footprint times coverage (so a pixel half inside counts half in Average); point kernels weight each tap by its coverage, and a target pixel gets a value when its sample point's pixel is at least half covered.

Cost

The cutline is rasterized once per warp call, over the source buffer: the centre rule with controlledburn's centre-rule burn, all touched and exact with its sparse coverage burn, and the blend ramp only within the blend distance of each edge. For a 14-vertex star over a 4096 x 4096 buffer that is about 15 ms (centre), 20 ms (all touched, exact) and 20 ms (blend distance 10) in a release build.

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

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

Not modelled

• The CENTER_LONG configuration option (--config CENTER_LONG), which GDAL applies to the target CRS as well as the source; CenterLong::At wraps the source only, as a CENTER_LONG in 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.

API documentation

• This version (built from main with this site): api/rwarp/
• Released versions: docs.rs/rwarp

The crate documentation goes into more detail than this guide on each module, including the GDAL source each part follows, with file and line references.

Changelog

All notable changes to rwarp. Dates are when the change reached main.

Unreleased

Added

• A spherical measure for conservative weights (stage 3a of design/conservative.md): WeightOptions::with_measure with weights::lonlat_measure(src_gt) measures areas on the sphere for a source on a regular longitude/latitude grid, through controlledburn's Measure::CosLatitude (closed-form edge integrals of cos(lat); the cell partition and conservation are unchanged). Weights records the measure, Weights::src_area(k) gives an entry's source cell measure, and normalised, dst_covered and apply weight overlaps by it; for the default planar measure every cell is 1 and results are unchanged, bit for bit. Exact against closed-form lon/lat box overlaps on the sphere (to 1e-14). Needs controlledburn 0.2.4. (2026-10-29)
• ESMF oracle suite (stage 2 of design/conservative.md, tests/esmf.rs): six cases, each against ESMPy 8.9.1's first-order conservative weights in Cartesian source pixel space (rwarp's cell definition: equal to about 1e-12) and on the sphere with great-circle edges (differences recorded as ceilings to shrink). Fixtures from tests/fixtures/esmf/gen.py, run in ghcr.io/hypertidy/gdal-r-esmf by .github/workflows/esmf-fixtures.yml. (2026-10-29)

Changed

• The serde feature also enables controlledburn's, for Measure in WeightOptions and Weights. (2026-10-29)

0.3.1 (2026-10-28)

Added

• Longitude wrapping of a geographic source, GDAL's CENTER_LONG: GenImgProjTransformer::with_source_size lets the transformer wrap source longitudes around the source raster's centre as gdalwarp does (InsertCenterLong, OGRProjCT), so a source stored as 170..190, -190..-170 or 0..360 is found by a target crossing 180 degrees. CenterLong (Auto, the default; Off; At(lon), #[non_exhaustive]) with with_center_long, with_target_extent for gdalwarp's -te, center_long(), source_center_long for GDAL's centre, and crs::is_geographic_degrees. A transformer built without with_source_size wraps nothing, as before. Checked against GDAL 3.14.0dev on 11 cases across 180 degrees (LAEA, Mercator, polar stereographic, lon/lat targets; chunked and whole): centres, exact transformer points, chunk lists and every output pixel, with both CRS backends. The LAEA chunk fixtures are back across 180 degrees. The proj feature now also depends on proj-sys (already a dependency of proj), to ask PROJ whether a CRS is geographic.
• weights module: conservative remapping weights from the exact footprint, as values (stage 1 of design/conservative.md). exact_weights returns Weights, one entry per overlapping (target cell, source cell) pair with full-grid indices and the overlap as an unrounded f64 fraction of the source cell, plus each target cell's outline area in source pixel space (dst_area, NaN where it failed). Weights::normalised gives ESMF's S under Normalization::FracArea or DestArea; Weights::apply applies the weights to one or more bands, skipping source nodata; Weights::esmf_row_col gives ESMF/SCRIP's 1-based flattened indices. exact_weights_with streams them per target cell. WeightOptions (densify, target offset for tiling, drop and snap tolerances, f32 rounding; WeightOptions::warp_compatible), ApplyOptions and GridWindow, all #[non_exhaustive] with builders; serde derives behind serde. apply(FracArea) with warp-compatible weights reproduces warp() with Footprint::Exact and Average bit for bit. warp() shares the lattice and per-pixel coverage step and its output is unchanged. (2026-10-28)
• GDAL's backmap geolocation inverse, as an option beside the quad index: GeolocOptions::with_inverse(GeolocInverse::Backmap) (GeolocInverse: QuadIndex, the default; Backmap; BackmapEstimate for GDAL_GEOLOC_USE_MAX_ACCURACY=NO; #[non_exhaustive]), with with_backmap_oversample (GEOLOC_BACKMAP_OVERSAMPLE_FACTOR) and GeolocArrays::backmap_grid. Ported from GDAL 3.14 in its arithmetic order (float32 backmap, GDALFillNodata hole filling, the spiral search), by way of wxylib's wxy-warp (geoloc.rs, GDAL 3.12; the geolocation code is unchanged since), adding PIXEL_OFFSET / STEP. Checked against GDAL 3.14.0dev on 9 new cases (swath by centres and by corners, lon/lat and polar stereographic targets, a polar ring and a swath across the antimeridian, holes in the arrays with source nodata, overlapping scans): success and position at every target pixel, the backmap's estimate, and 24 warps (near, bilinear, cubic; -et 0 and 0.125) identical to gdalwarp, every pixel, with proj4rs and the fixtures' PROJ (with an older libproj, nearest identical and interpolated values within 2e-8 relative). Where the two inverses differ, and why, is in the curvilinear guide. (2026-10-28)

Changed

• The chunk fixtures' special points are computed as GDAL computes them for a geographic source (not reprojected; see ChunkOptions::special_points). No chunk list changes.
• warp() does GDAL's GWKAdjustSrcOffsetOnEdge: when a corner of the source buffer fails to transform forward and a sample point lands on an invalid source pixel whose corners fail too, a valid neighbour (right, below, left, above) is used instead. This changes output only where the transformer has holes or edges inside the source, as geolocation arrays with missing samples do (12 of 4560 pixels in the holes test case, now identical to gdalwarp). warp() therefore calls the transformer source to destination too, for such pixels only. Point methods; GDAL's variant for alpha and cutline densities is not reproduced. (2026-10-28)
• The geolocation swath fixture is regenerated with GDAL 3.14.0dev (was 3.8.4): positions moved by up to 3e-8 pixels (PROJ), nothing else.
• controlledburn 0.2.3 (corner passes, near-coincident edges, exact centre ties). The cutline centre rule now matches gdalwarp on centres exactly on diagonal edges: 5 pixels of the mask_diamond_tie fixture that differed now agree. They came from controlledburn 0.2.2's rounded intercepts; the parity page had put them down to GDAL.

0.3.0 (2026-10-28)

This release changes the API: options types are #[non_exhaustive].

Changed

• warp() retries a sample point that lands less than one pixel outside the source buffer by transforming it again on its own, as GDAL does (GWKCheckAndComputeSrcOffsets). Through ApproxTransformer that is exact, so interpolation error no longer drops pixels at the edge of a chunk's source window; with an exact transformer nothing changes. Point methods only (nearest, bilinear, cubic, cubic spline, lanczos), as in GDAL. (2026-10-28)
• WarpOptions, GeolocOptions and AverageThresholds are #[non_exhaustive], as are ResampleAlg, KernelScaling, Footprint, BandNodata, Nodata, DstNodata, DataType and GeolocConvention. New options and methods can then come in minor releases. Outside rwarp, build options with WarpOptions::new(alg) and the with_ methods (GeolocOptions::default().with_convention(..), and so on) instead of struct literals, and add a _ arm when matching on the enums. Fields stay public for reading and assignment.

Added

• WarpOptions::new and with_ builder methods on WarpOptions, GeolocOptions and AverageThresholds.
• plan module: plan_chunks and ChunkPlanner give gdalwarp's chunk list as values (WarpPlan of Chunks: target window, source window, fill ratio, kernel scale), following GDAL 3.14's CollectChunkList and ComputeSourceWindow exactly: warp memory limit (-wm) and per-pixel cost, fill ratio splits, OPTIMIZE_SIZE on target blocks, special points, the forward probe, SAMPLE_GRID, SAMPLE_STEPS, SOURCE_EXTRA. ChunkOptions (builders, #[non_exhaustive]); ChunkOptions::gdalwarp for gdalwarp's planning exactly; an optional round-trip check for backends that do not fail at a projection's antipode; the antimeridian split of 0.2's collect_chunk_list as split_antimeridian. warp_chunks warps a plan from a source buffer in memory. Ported from wxylib's wxy-warp (chunked.rs). (2026-10-28)
• serde feature: Serialize/Deserialize for WarpPlan and Chunk. (2026-10-28)
• Parity fixtures for whole rasters warped in many chunks (tests/chunks.rs, tests/fixtures/chunks/): seven gdalwarp cases, chunk lists and output identical. (2026-10-28)
• Cutlines (gdalwarp -cutline): WarpOptions::cutline / with_cutline, and the cutline module with Cutline (polygons in source pixel/line), CutlineRule::Centre (GDAL's default, with its tie rules), CutlineRule::AllTouched (CUTLINE_ALL_TOUCHED), blend distance (-cblend), and CutlineRule::Exact, which weights each source pixel by the fraction of it the cutline covers (not a GDAL mode). Cutline::from_geotransform and Cutline::from_transformer (densified as gdalwarp does) map world geometry to source pixels; Cutline::density is the per-pixel density the warp uses. The cutline is rasterized over the source buffer only, with controlledburn.

Removed

• source_window::collect_chunk_list and ChunkPlan: use plan::plan_chunks (with ChunkOptions::split_antimeridian for the same split at the antimeridian). compute_source_window stays. (2026-10-28)

0.2.0 (2026-10-28)

This release changes the API: warp() with WarpOptions and SourceBuffer replaces the warp entry points of 0.1, and ResampleAlg gains variants.

Added

• geoloc_transform: geolocation-array transformers, for a curvilinear source (gdalwarp -geoloc), a curvilinear target (DST_METHOD=GEOLOC_ARRAY), or both. Both of GDAL's georeferencing conventions; offsets and steps; nodata; antimeridian handling. (2026-10-28)

• warp(), WarpOptions, SourceBuffer: one entry point for every sample type, method and option. WarpOptions carries a destination offset for warping one tile of a larger grid. (2026-10-28)

• KernelScaling: GDAL's kernel widening when downsampling, now the default; Fixed and Point alternatives. estimate_scale is public. (2026-10-28)

• ResampleAlg::Average, with Footprint::GdalBox (gdalwarp's) and Footprint::Exact (first-order conservative). controlledburn is a new dependency. (2026-10-28)

• ResampleAlg::Mode with ModeTies (MODE_TIES). (2026-10-28)

• BandNodata (UNIFIED_SRC_NODATA) and AverageThresholds (NODATA_VALUES_PCT_THRESHOLD, EXCLUDED_VALUES). (2026-10-28)

• ResampleAlg::CubicSpline, Min, Max, Med, Q1, Q3 and Rms, and ResampleAlg::is_area(). (2026-10-28)

• Sample for i64, u64, and half::f16 behind a new f16 feature. (2026-10-28)

• warp_to::<S, D>: a different output type from the source's, through GDAL's working type (DataType::working) and GDALCopyWords conversion. (2026-10-28)

• Nodata and DstNodata: source and destination nodata per band, as gdalwarp's -srcnodata and -dstnodata. (2026-10-28)

• Examples: tile_laea, curvilinear_swath.

• A documentation site, with the browser demo under /demo/.

Changed

• Integer outputs round half up and clamp, and a resampled value equal to nodata is moved off it, as gdalwarp does with an integer working type.
• Sample points outside [0, n) of the source are off the source for every point method, as in GDAL.
• Lanczos clips its window and renormalises at edges and nodata, and a target pixel whose centre lands on nodata gets no value, for every interpolating kernel. Some edge and nodata pixels change compared with 0.1.
• Multi-band nodata follows GDAL's default (a pixel is nodata only where every band is); the 0.1 rule is BandNodata::AnyBand.
• Sample requires PartialOrd.
• Footprint::Exact is 2 to 6x faster. The destination pixel lattice is transformed once per row band (shared corners and edge points) instead of four points per pixel, and coverage comes from controlledburn::CoverageWindow (0.2.2) instead of one burn call per pixel. With an exact transformer the output is unchanged, except for pixels straddling the antimeridian, which the exact footprint does not handle (it gave a different wrong answer before). The lattice's row lines go to the transformer as scanlines, so an ApproxTransformer now approximates them, as for the box footprint; before, it transformed the exact footprint's points exactly.
• WarpOptions is not generic: nodata: Option<T> is replaced by src_nodata: Nodata and dst_nodata: DstNodata.
• Bilinear uses GDAL's general formula for types and cases where GDAL does, which can change the last bit (and so integer rounding near 2^40 for 64-bit values).

Deprecated

• warp_resample_t and warp_resample_u8, kept as shims with the 0.1 point-kernel behaviour. Use warp().

Fixed

• Nearest neighbour mapped source coordinates in (-1, 0) to pixel 0.
• ApproxTransformer left one point untransformed and reported it failed when subdivision reached a left half of two points, which happens on strongly curved transforms or with a small error threshold (whole columns of a polar warp came out as nodata). GDAL transforms it.

0.1.2 (2026-03-11, not on crates.io)

• Use the system libproj through proj 0.31.

Not released to crates.io. On main after it (2026-09 to 2026-10):

• rwarp-wasm: browser shim and demo (OSM, WMTS, ArcGIS tile services, COG sources), Warper.withTransforms for JavaScript CRS transforms.
• A generic Sample kernel over u8 to f64, interleaved bands, nodata as a value or an alpha band.
• The CrsTransform trait with proj and proj4rs backends; builds for wasm32-unknown-unknown.

0.1.1 (2026-03-06)

• GcpTransformer: polynomial transformer from ground control points, order 1 to 3.
• Use the published vaster crate; crate documentation.

0.1.0 (2026-03-06)

• First release, extracted from cogcache: GenImgProjTransformer, ApproxTransformer, compute_source_window, collect_chunk_list and the nearest, bilinear, cubic and lanczos kernels.

Applications

Browser demo

Open the demo. A Leaflet map in a projection of your choice (a local Lambert azimuthal equal-area by default, centred wherever you like) whose tiles are warped client-side from web map tiles, WMTS, ArcGIS tile services or a Cloud Optimized GeoTIFF. rwarp runs in a Web Worker, compiled to WebAssembly with the pure-Rust proj4rs backend: no server, no libproj.

The shim is rwarp-wasm/ in the repository. Its JavaScript API plans a source window and warps RGBA or float tiles:

const w = new Warper(srcCrs, srcGt, dstCrs, dstGt, maxError);
const win = w.source_window(dstW, dstH, srcW, srcH, padding);
const rgba = w.warp_rgba(srcRgba, srcW, srcH, xoff, yoff, dstW, dstH, "nearest");

Warper.withTransforms accepts CRS transforms as JavaScript functions, so a page can lend rwarp a projection library (proj4js, a PROJ build) when proj4rs declines a CRS. See rwarp-wasm/README.md.

cogcache

rwarp was extracted from cogcache, which plans warps from a target tile back to the source chunks it needs, and caches both the raw encoded source chunks and the warped target tiles. rwarp is the planning and resampling core.

Reference for array planning

The ArrIR design for portable array operation plans defines resampling independently of GDAL, and uses rwarp and gdalwarp together to pin down what gdalwarp actually does: how ties are broken, how edges are placed for curvilinear grids, where kernels widen. rwarp's fixtures are the evidence for those answers.

Using rwarp in your project

If you build on rwarp, an issue or pull request adding your project here is welcome.