gdalwarp AVHRR in NOAA level 1b format looks like Dali painting

[arb]

I used this command to convert AVHRR in NOAA Level-1B format into a GeoTIFF (in the UTM31 projection with ED50 datum):

gdalwarp -tps -t_srs EPSG:23031 -r bilinear -tr 1000 1000 level1b.dat geo.tif

I believe that -tps is the appropriate option when using GCPs from the 1B input file.

The first result looked like a Dali painting, especially at the corners. I enlarged the input file hoping that the deformations would be less near my region of interest but the result is worse.

[arb]

First geotiff (small)

[arb]

Second geotiff (larger area)

[arb]

The sample files can be downloaded from ​ftp://ftp.sat.dundee.ac.uk/misc/testdata/geotiff/

[Even Rouault]

The ftp directory is not accessible. The access rights are missing for 'other'

[dron]

I believe this is related to the wrong handling of GCP points from the southern hemisphere (see the number of FIXME notes related to ascending/descending flag processing in l1bdataset.cpp). Sample file needed to fix the problem.

Best regards, Andrey

[arb]

Sorry, ftp access rights now corrected.

Additional info on interpolation (and extrapolation) of GCPs using the Lagrandian method as suggested by NOAA can be found here: ​http://www2.ncdc.noaa.gov/docs/klm/html/c2/sec2-4.htm Maybe something worth considering implementing in GDAL?

[dron]

My first guess was wrong: the GCP point are being extracted well and there is no problem here. The actual problem could be related to TPS algorithm. When trying to do the transformation of the first sample file I am getting the error:

$ gdalwarp -tps -t_srs EPSG:23031 -r bilinear -tr 1000 1000 geo-1b-a.l1b geo-1b-a.tif

GDAL: GDALOpen(geo-1b-a.l1b) succeeds as L1B.
ERROR 1: Unable to compute a TPS based transformation between pixel/line
and georeferenced coordinates for geo-1b-a.l1b.

Everything works fine without "-tps" switch and I do not have any distortions on my image. Strange that I can't reproduce it with my copy of GDAL (taken from trunk). What is your GDAL version?

Best regards,
Andrey

[dron]

Warping of the sample NOAA L1B dataset using default algorithm

[dron]

Talking on warping and interpolation methods I should say that I am realizing that our current approaches are not quite suitable for NOAA. It is desirable to have implementation of the NOAA specific algorithm. But it should be either contributed by someone using this data or funded, because I just do not have a free time for that (and I am not using these data myself anymore).

Best regards,
Andrey

[arb]

I had trouble installing gdal so I'm using the latest FWtools (which is great, by the way, because it's easy to install, so you should make nightly builds of it on all platforms!). Version 2.0.6 for linux seems to contain libgdal-1.6.0 if that helps.

The gdal manual warned me to stay away from the default interpolation and use "-tps".

You say you have no distortions, well no Dali-esque ones, but look at Spain, it is clearly distorted! Did you specify the same projection?

Attached is some code to perform the interpolation. Does that help?

[arb]

Lagrange interpolation for NOAA Level 1B geo-location

[dron]

Replying to arb:

You say you have no distortions, well no Dali-esque ones, but look at Spain, it is clearly distorted! Did you specify the same projection?

You can find my options in the message above. It is the same command line as in your initial message with except of "-tps" option.

The Spain is slightly shrinked, yes. That is because it is located at the border of the scene where pixels became larger. L1B driver reports just a subset of the total amount of GCPs. Look at gdalinfo output: the first GCP appears at the 22nd pixel. That should be changed: we need more GCPs near the scene edge and lesser at the image center. I will fix this problem.

Attached is some code to perform the interpolation. Does that help?

Well, actually I would like to have a patch against the GDAL code with some preliminary discussion on what part of GDAL should be improved in which way.

[warmerdam]

Please ensure that the reporter isn't just running into #2300, now fixed in trunk and 1.5.

[Even Rouault]

Yes, absolutely. 'arb', you're are certainly hitting #2300 as I've no problem with my GDAL tree. See the screenshot attached.

Andrey, I've just fixed , in r14530, the '-tps' switch in gdalwarp.cpp that was broken recently in trunk only due to changes in gdaltransform.cpp.

[Even Rouault]

Result of "gdalwarp -tps -t_srs EPSG:23031 -r bilinear -tr 1000 1000 geo-1b-a.l1b geo-1b-a-tps.tif" with r14530

[arb]

There is a difference between the two latest images, but if they used the same cmd line (except -tps) then one of them must be wrong! To me this one simply looks wrong: "geo-1b-a-dron.jpg - added by dron on 05/26/08 11:44:09". This one looks correct: "geo-1b-a-tps-r14530.jpg - added by rouault on 05/26/08 14:40:31". If the only difference is that the latter uses "-tps" then the former must be doing some broken interpolation.

I realise the code I sent is not a patch, that is because I'm not a GDAL developer, just a user! But I'm sure a GDAL developer could easily use it to return a geolocation array instead of a set of GCPs.

Other datasets (notably HDF) contain their own geolocation arrays. In fact some HDF datasets contain a sparse array, eg. MODIS 250m bands only have 1000m geolocation. Can you just copy the code for HDF? Then you can use my code to dilate the spare geolocation array using the Lagrange method.

[dron]

r14530 really fixes TPS algorithm for me. I am having the same results as Even now. The image, attached by me looks different from the Even one because default warping method just not suitable for AVHRR. Use TPS, it works now.

Also I started adding a geolocation array extraction to L1B driver a long time ago but this work wasn't finished due to the lack of time, as usual (and it is not as simple as just copy the code from another driver). Could you create a separate enhancement request for that and assign it to me, so it will not be forgotten anymore? AFAIK HDF driver is the only driver using this feature yet.

Interpolation of geolocation points is not an option for GDAL drivers. GDAL reports the dataset content without any processing inside the driver. If you need to process this data in some way it should be done outside the driver. You can obtain the geolocation array from GDAL and apply interpolation algorithm somwhere in your processing toolchain.

[Even Rouault]

Closing now. The issue was with the broken TPS that has been solved.

#2403 tracks the enhancement request for geolocation array.