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

[Brendan Heberlein]

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


Hi,


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


Even


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),
                      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



_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/gdal-dev


--
http://www.spatialys.com
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>