Lacking interpolation of geolocation arrays

[knutfrode]

Warping datasets with geolocation arrays works fine if the arrays have the same size as the other raster bands. However, if the geolocation arrays are much smaller, the quality is very bad.

Say I have a 2000x2000 px dataset with 100x100 px sized geolocation arrays. Then a warped image have only 200x200 blocks of identical values, even if the size of the warped image is much larger. If I change the interpolation algorithm (neareast, cubic...) the values change slightly, but the real output resolution remains low. This is the same whether using Python API or gdalwarp. As a workaround I convert the geolocation arrays to GCPs. This gives decent warping quality, but is a bit awkward.

(from gdal mailing list: ​http://thread.gmane.org/gmane.comp.gis.gdal.devel/33336/focus=33339)

[knutfrode]

A Tiff-image and a Python script which illustrates this bug is attached.

[Even Rouault]

trunk r25403 "geoloc: add bilinear interpolation of coordinates from backmap (#4907)"

[Even Rouault]

trunk r25404 "Adjust expected results after r25403 (#4907)"

[knutfrode]

This interpolation is a big improvement, but I believe there must be a small bug in the interpolation code.

Running the attached example script (above), one can see some artefacts at some locations in the image, probably around the lon/lat geolocation pixels. Increasing LINE_STEP and PIXEL_STEP in the script to 30 or higher, the artefacts are dominating, as also seen with some real case satellite images which cannot be attached here.

I am not sure, but have a feeling that the interpolation algorithm sometimes thinks it is on the boundary of the image when it is not.

[knutfrode]

I have attached two real-case examples of AVHRR images containing geolocation arrays, which have been warped (to mercator projection) with -geoloc option.

In one case the geolocation array was defined every 40th pixel, but for every line (all that is available). In the second case the line sampling was also reduced to every 40th line - so that the geolocation array is more evenly spaced.

This gives a hint what is wrong with the interpolation: as longitudes and latitudes do not make an orthogonal coordinate system (except at equator), the bilinear interpolation gives stretching, which becomes worse at higher latitudes. This should not be a problem if the geolocation coordinates were in an equidistant projected system, but I guess that longitudes and latitudes will be the most common by far.

Thus the interpolation should probably rather be done in an equidistant coordinate system. Alternatively, polynomial fitting should probably work better. That should also allow for extrapolation, for the cases where the geolocation arrays do not cover the boundaries (which they do not for AVHRR at least).

[Even Rouault]

A more recent commit probably fixed the issue. So closing

Geoloc transformer: fix wrong bilinear interpolation in GDALGeoLocTransform() (#5305)