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

Brendan Heberlein bheberlein at wisc.edu
Sat Jun 12 21:15:44 PDT 2021

Hi Even,

thanks for your response. I will try to test this with the updated code, although I have yet to learn the compilation process.

I will let you know if I am able to verify this is fixed.

All the best,
     — Brendan Heberlein
            Research Intern, BS
            Townsend Lab
            University of Wisconsin - Madison

From: Even Rouault <even.rouault at spatialys.com>
Sent: Friday, June 11, 2021 2:49 PM
To: Brendan Heberlein <bheberlein at wisc.edu>; gdal-dev at lists.osgeo.org <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Edge Artifacts with GDAL Warp & Geolocation Arrays


I believe https://github.com/OSGeo/gdal/pull/3964 should improve this


Le 11/06/2021 à 21:30, Brendan Heberlein via gdal-dev a écrit :
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),

Any insight would be much appreciated. Cheers!

     — Brendan Heberlein
            Research Intern, BS
            Townsend Lab
            University of Wisconsin - Madison

gdal-dev mailing list
gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210613/19437f15/attachment-0001.html>

More information about the gdal-dev mailing list