geoloc backmap interpolation only doing one real iteration.

[warmerdam]

On Tue, 22 Jul 2008 17:18:18 +0100, Frank Warmerdam <warmerdam@…> wrote:

I imagine this is a limitation of the "hole filling" code:

http://trac.osgeo.org/gdal/browser/trunk/gdal/alg/gdalgeoloc.cpp#L267

You might want to try changing this code to go more than 3 iterations and see if it helps.

The loop serves no purpose at the moment (other than wasting some time!) because iterations after the first simply do the same thing again. I think you want to mark a newly-filled pixel as valid, for example something like

	pabyValidFlag[iBMX+iBMY*nBMXSize] = 1 ??

Obviously I don't understand the code because I found that increasing the size of the backmap (currently 1.3 times bigger than something else) made the holes worse. Using a factor of 0.5 actually gave me an image without any holes at all! Maybe the output image wasn't actually correct but at first glance it seemed fine :-)

Andrew

[warmerdam]

I have tried the obvious change to mark the interpolated values as valid, but this produces bad smear-out effects at image edges. I think there needs to be some more logic in the code using the backmap.

I'm shelving this for a bit but it needs to be fixed before 1.6 release.

(ps the smearing effect is clear in the autotest/gcore/geoloc.py output images)

[arb]

The attached code exhibits less obvious smearing effects, although if too many iterations are performed then it will start smearing at the edges. However it should be better at upto 3 iterations. It also checks diagonal neighbours which may or may not be a sensible thing to do. It breaks early if nothing needs to be done. (The code still has debugging statements included.)

[arb]

(first attempt to fix loop)

[arb]

It turns out that the code I attached seems to solve the problems I was experiencing.

There seems to be no need to increase the number of iterations beyond 3.

There seems to be no need to change the backmap size factor from 1.3, although doing so does produce noticable changes in the output.

[arb]

Final patch which fixes the problem

[Mateusz Łoskot]

Frank,

I was about to apply Andrew's patch but I think you may want to review it first. Is the fix OK?

[Markus Neteler]

As co-sponsor of the TPS method I am still interested after years to get AVHRR/NOAA data properly geocoded. According to ​http://lists.osgeo.org/pipermail/gdal-dev/2009-May/020696.html this patch helps.

[Even Rouault]

I've applied Andrew's patch in trunk r17015, which just stylistic changes, and after a quick reviewing to check that it is not fundamentally wrong. But I'm not too familiar with geoloc, so it would be great if other people (Markus ?) could confirm that it, at least, doesn't make things worse in situations where it used to work before. I've just visually checked with the test in autotest/gcore/geoloc.py that it didn't change things too much (expected checksum changed in r17016).

I'm attaching the result of autotest/gcore/geoloc.py before and after.

[Markus Neteler]

I have made tests with my long-term test AVHRR scene:

gdalinfo A1726874.L2460388
Driver: L1B/NOAA Polar Orbiter Level 1b Data Set 
Files: A1726874.L2460388    
Size is 2048, 247    
Coordinate System is `'
GCP Projection = GEOGCS["WGS 72",DATUM["WGS_1972",SPHEROID["WGS 72",6378135,298.26,AUTHORITY["EPSG",7043]],
 TOWGS84[0,0,4.5,0,0,0.554,0.2263],AUTHORITY["EPSG",6322]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
 UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",9108]],AUTHORITY["EPSG",4322]]   
GCP[  0]: Id=, Info=   
   (25.5,0.5) -> (-5.7777,48.9182,0) 
GCP[  1]: Id=, Info=   
   (225.5,0.5) -> (1.5265,48.7581,0) 
GCP[  2]: Id=, Info=   
   (105.5,0.5) -> (-2.1543,48.8966,0)
GCP[  3]: Id=, Info=   
   (625.5,0.5) -> (8.4579,48.1722,0) 
GCP[  4]: Id=, Info=   
   (185.5,0.5) -> (0.4539,48.8106,0) 
GCP[  5]: Id=, Info=   
   (225.5,0.5) -> (1.5265,48.7581,0) 
GCP[  6]: Id=, Info=   
   (1225.5,0.5) -> (15.0325,47.2071,0)
GCP[  7]: Id=, Info=   
   (305.5,0.5) -> (3.359,48.6449,0)  
GCP[  8]: Id=, Info=   
   (345.5,0.5) -> (4.1565,48.5863,0) 
GCP[  9]: Id=, Info=   
   (385.5,0.5) -> (4.8925,48.5273,0) 
GCP[ 10]: Id=, Info=   
   (2025.5,0.5) -> (29.7619,43.4221,0)
GCP[ 11]: Id=, Info=   
   (25.5,13.5) -> (-5.7829,48.7936,0)
GCP[ 12]: Id=, Info=   
   (225.5,13.5) -> (1.5025,48.6319,0)
GCP[ 13]: Id=, Info=   
   (105.5,13.5) -> (-2.1688,48.771,0)
...
GCP[216]: Id=, Info=   
   (305.5,2047.5) -> (2.851,46.3723,0)
GCP[217]: Id=, Info=   
   (345.5,2047.5) -> (3.6142,46.3139,0)
GCP[218]: Id=, Info=   
   (385.5,2047.5) -> (4.3188,46.2554,0)
GCP[219]: Id=, Info=   
   (2025.5,2047.5) -> (28.3301,41.3809,0)   
Metadata:
  DATASET_NAME=NSS.LHRR.NL.D01326.S0139.E0145.B0601717.WI 
  SATELLITE=NOAA-16(L) 
  DATA_TYPE=AVHRR LAC  
  REVOLUTION=06017
  SOURCE=Wallops Island, Virginia, USA
  PROCESSING_CENTER=NOAA/NESDIS - Suitland, Maryland, USA 
  START=year: 2001, day: 326, millisecond: 5993085 
  STOP=year: 2001, day: 326, millisecond: 6034085  
  LOCATION=Descending  
Corner Coordinates:    
Upper Left  (    0.0,    0.0) 
Lower Left  (    0.0,  247.0) 
Upper Right ( 2048.0,    0.0) 
Lower Right ( 2048.0,  247.0) 
Center( 1024.0,  123.5) 
Band 1 Block=2048x1 Type=UInt16, ColorInterp=Undefined    
  Description = AVHRR Channel 1:  0.58  micrometers -- 0.68 micrometers 
Band 2 Block=2048x1 Type=UInt16, ColorInterp=Undefined    
  Description = AVHRR Channel 2:  0.725 micrometers -- 1.10 micrometers 
Band 3 Block=2048x1 Type=UInt16, ColorInterp=Undefined    
  Description = AVHRR Channel 3B: 3.55  micrometers -- 3.93 micrometers 
Band 4 Block=2048x1 Type=UInt16, ColorInterp=Undefined    
  Description = AVHRR Channel 4:  10.3  micrometers -- 11.3 micrometers 
Band 5 Block=2048x1 Type=UInt16, ColorInterp=Undefined    
  Description = AVHRR Channel 5:  11.5  micrometers -- 12.5 micrometers 

# bilinear
gdalwarp -t_srs EPSG:4326 -r bilinear -tr 1000 1000 A1726874.L2460388 A1726874.L2460388_LL.tif
Creating output file that is 0P x 0L.
ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

# TPS
gdalwarp -tps -t_srs EPSG:4326 -tr 1000 1000 A1726874.L2460388 A1726874.L2460388_LL.tif
 There is a problem to invert the interpolation matrix
 There is a problem to invert the interpolation matrix
ERROR 1: Too many points (441 out of 441) failed to transform,
unable to compute output bounds.

# default
gdalwarp A1726874.L2460388 A1726874.L2460388_LL.tif
Creating output file that is 2085P x 489L.
Processing input file A1726874.L2460388.
0...10...20...30...40...50...60...70...80...90...100 - done.
# -> Result looks quite "funky" and is completely unusable

In ticket #2388 is lagrange.c attached (Lagrange interpolation for NOAA Level 1B geo-location) which might solve it for sparse GCPs as appear in this image.

[Markus Neteler]

This AVHRR can be downloaded here (1.8MB): ​http://gis.fem-environment.eu/outgoing/A1726874.L2460388.gz

[Markus Neteler]

Sorry for posting to the wrong ticket (I have GCPs not geolocation), shifted report to ticket #2403

[Even Rouault]

Patch has been applied. Closing