Opened 16 years ago

Closed 13 years ago

#2501 closed defect (fixed)

geoloc backmap interpolation only doing one real iteration.

Reported by: warmerdam Owned by:
Priority: normal Milestone:
Component: GDAL_Raster Version: svn-trunk
Severity: normal Keywords: geolocation
Cc: Markus Neteler, dron

Description

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

Attachments (4)

gdalgeoloc-arb.cpp (37.9 KB ) - added by arb 16 years ago.
(first attempt to fix loop)
gdalgeoloc.diff (5.6 KB ) - added by arb 16 years ago.
Final patch which fixes the problem
result_of_autotest_before_patch.tif (6.3 KB ) - added by Even Rouault 15 years ago.
result_of_autotest_after_patch.tif (6.3 KB ) - added by Even Rouault 15 years ago.

Download all attachments as: .zip

Change History (15)

comment:1 by warmerdam, 16 years ago

Milestone: 1.6.0
Status: newassigned

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)

comment:2 by arb, 16 years ago

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

by arb, 16 years ago

Attachment: gdalgeoloc-arb.cpp added

(first attempt to fix loop)

comment:3 by arb, 16 years ago

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.

by arb, 16 years ago

Attachment: gdalgeoloc.diff added

Final patch which fixes the problem

comment:4 by Mateusz Łoskot, 16 years ago

Frank,

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

comment:5 by Markus Neteler, 15 years ago

Cc: Markus Neteler added

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.

comment:6 by Even Rouault, 15 years ago

Milestone: 1.6.11.7.0
Owner: changed from warmerdam to Even Rouault
Status: assignednew

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.

by Even Rouault, 15 years ago

by Even Rouault, 15 years ago

comment:7 by Markus Neteler, 15 years ago

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.

comment:8 by Markus Neteler, 15 years ago

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

comment:9 by Even Rouault, 15 years ago

Cc: dron added
Owner: Even Rouault removed

comment:10 by Markus Neteler, 15 years ago

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

comment:11 by Even Rouault, 13 years ago

Milestone: 1.8.1
Resolution: fixed
Status: newclosed

Patch has been applied. Closing

Note: See TracTickets for help on using tickets.