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.