[gdal-dev] Edge Artifacts with GDAL Warp & Geolocation Arrays

[Brendan Heberlein]

Hello all,

I have been using gdalwarp​ to georeference sensor data with geolocation arrays.

I have observed that GDAL assigns values to pixels that clearly fall outside the geolocation footprint, producing some funny artifacts at the image edges. For example, see the image here<https://pasteboard.co/K69DUni.png>.

(Sensor data over Google satellite basemap; red dots are geolocation points.)

Notice how the data values are 'extruded' from the edges.

This occurs with various interpolation methods (all the ones I tried), including nearest neighbor, average & bilinear. Why does this happen? How can I go about fixing this?

Ideally, any pixel in the output grid which lies more than 1 pixel width from any of the geolocation points (red dots) would be assigned to no-data.

My code is very straightforward, basically I open the raw sensor image & geolocation array (a.k.a. input geometry map or IGM) with GDAL & then assign geolocation metadata to the sensor image before warping:
# Geolocation domain metadata
geoloc_metadata = {
    'SRS': WGS84_UTM16N.ExportToWkt(),
    'X_BAND': '1',
    'X_DATASET': str(igm_img_file),
    'Y_BAND': '2',
    'Y_DATASET': str(igm_img_file),
    'PIXEL_OFFSET': '0',
    'LINE_OFFSET': '0',
    'PIXEL_STEP': '1',
    'LINE_STEP': '1',
}

# Set geolocation domain metadata on dataset to be warped
dn_img_ds.SetMetadata(geoloc_metadata, 'GEOLOCATION')

x_min, y_min = igm[:2].min(axis=(1,2))
x_max, y_max = igm[:2].max(axis=(1,2))

# warp_options = gdal.WarpOptions(warpMemoryLimit=)
warped_ds = gdal.Warp(f'/warped.tif', dn_img_ds,
                      outputBounds=(x_min, y_min, x_max, y_max),
                      xRes=1.0,
                      yRes=1.0,
                      resampleAlg='near',
                      dstSRS=WGS84_UTM16N,
                      srcNodata=None,
                      dstNodata=-9999,
                      multithread=True,
                      targetAlignedPixels=True,
                      outputType=get_gdal_dtype(np.int16),
                      geoloc=True)

Any insight would be much appreciated. Cheers!

     — Brendan Heberlein
            Research Intern, BS
            Townsend Lab
            University of Wisconsin - Madison
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210611/9c552243/attachment.html>