Opened 9 years ago

Closed 9 years ago

Last modified 8 years ago

#2388 closed defect (duplicate)

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

Reported by: arb Owned by: dron
Priority: normal Milestone:
Component: default Version: unspecified
Severity: critical Keywords:
Cc: dron, Markus Neteler

Description

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.

Attachments (5)

geo-1b-a.jpg (254.5 KB) - added by arb 9 years ago.
First geotiff (small)
geo-1b-b.jpg (423.2 KB) - added by arb 9 years ago.
Second geotiff (larger area)
geo-1b-a-dron.jpg (110.8 KB) - added by dron 9 years ago.
Warping of the sample NOAA L1B dataset using default algorithm
lagrange.c (3.4 KB) - added by arb 9 years ago.
Lagrange interpolation for NOAA Level 1B geo-location
geo-1b-a-tps-r14530.jpg (36.8 KB) - added by Even Rouault 9 years ago.
Result of "gdalwarp -tps -t_srs EPSG:23031 -r bilinear -tr 1000 1000 geo-1b-a.l1b geo-1b-a-tps.tif" with r14530

Download all attachments as: .zip

Change History (19)

Changed 9 years ago by arb

Attachment: geo-1b-a.jpg added

First geotiff (small)

Changed 9 years ago by arb

Attachment: geo-1b-b.jpg added

Second geotiff (larger area)

comment:1 Changed 9 years ago by arb

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

comment:2 Changed 9 years ago by Even Rouault

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

comment:3 Changed 9 years ago by dron

Cc: dron added

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

comment:4 Changed 9 years ago by 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?

comment:5 Changed 9 years ago by 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

Changed 9 years ago by dron

Attachment: geo-1b-a-dron.jpg added

Warping of the sample NOAA L1B dataset using default algorithm

comment:6 Changed 9 years ago by 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

comment:7 Changed 9 years ago by 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?

Changed 9 years ago by arb

Attachment: lagrange.c added

Lagrange interpolation for NOAA Level 1B geo-location

comment:8 in reply to:  7 Changed 9 years ago by dron

Owner: changed from warmerdam to dron
Status: newassigned

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.

comment:9 Changed 9 years ago by warmerdam

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

comment:10 Changed 9 years ago by 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.

Changed 9 years ago by Even Rouault

Attachment: geo-1b-a-tps-r14530.jpg added

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

comment:11 Changed 9 years ago by 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.

comment:12 Changed 9 years ago by 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.

comment:13 Changed 9 years ago by Even Rouault

Resolution: duplicate
Status: assignedclosed

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

#2403 tracks the enhancement request for geolocation array.

comment:14 Changed 8 years ago by Markus Neteler

Cc: Markus Neteler added
Note: See TracTickets for help on using tickets.