Opened 9 years ago

Closed 4 years ago

#2403 closed defect (fixed)

[PATCH] NOAA level 1b geolocation interpolation

Reported by: arb Owned by: dron
Priority: normal Milestone: 1.10.2
Component: GDAL_Raster Version: svn-trunk
Severity: normal Keywords: AVHRR L1B
Cc: neteler, jurien, spareeth, mmetz

Description

The current noaa level 1b geolocation relies on sparse GCPs. This might be fine when reprojecting a small area but when reprojecting a whole swath width it produces significant errors.

It might be better if the GCPs were interpolated inside the driver and returned as a geolocation array, or if a new Lagrangian interpolation method was added to GDAL.

See ticket 2388: http://trac.osgeo.org/gdal/ticket/2388 although the sample code there doesn't handle the case where longitude wraps around (360 back to 0), or switches sign (180 to -180). When copying the 51 longitudes into a 2048-element array (prior to interpolation) it would be advisable to check for a difference > 180 compared to the previous longitude and just continue with the same magnitude/sign, then perform some modulo function after interpolation. Hope that helps!

Attachments (2)

noaa_l1b_gcp_equidistant_subsampling.patch (1.7 KB) - added by ramiro 6 years ago.
Modification to L1BDataset::ProcessRecordHeaders? to return equidistant points.
l1bdataset.diff (3.8 KB) - added by mmetz 4 years ago.
patch for frmts/l1b/l1bdataset.cpp to read more GCPs

Download all attachments as: .zip

Change History (15)

comment:1 Changed 9 years ago by dron

Status: newassigned

comment:2 Changed 8 years ago by neteler

Cc: neteler added

comment:3 Changed 8 years ago by neteler

Component: defaultGDAL_Raster
Milestone: 1.7.0
Version: unspecifiedsvn-trunk

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.

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

Changed 6 years ago by ramiro

Modification to L1BDataset::ProcessRecordHeaders? to return equidistant points.

comment:4 Changed 6 years ago by ramiro

Function L1BDataset::ProcessRecordHeaders subsamples the GCPs of the l1b image, but the points returned are not equidistant.
The attached patch modifies this behaviour.

I tested "gdalwarp -tps" on the AVHRR image and seems to work ok(no error "There is a problem to invert the interpolation matrix" produced).

gdalwarp -tps -t_srs EPSG:4322 A1726874.L2460388 out_noaa_l1b.tif

To apply the patch execute from gdal:

patch -p0 < noaa_l1b_gcp_equidistant_subsampling.patch

comment:5 Changed 6 years ago by jurien

Cc: jurien added

comment:6 Changed 6 years ago by Even Rouault

Milestone: 1.8.1
Summary: NOAA level 1b geolocation interpolation[PATCH] NOAA level 1b geolocation interpolation

comment:7 in reply to:  6 Changed 6 years ago by jurien

Replying to rouault:

Does this mean there are no plans to include it on a future release?

Regards

comment:8 Changed 6 years ago by Even Rouault

No, it just means that as 1.8.1 has been released and didn't include this patch, it was no longer the appropriate milestone. And I didn't assign a new one because - to the extent of my knowledge - there's no one actively working on this bug. But the ticket is still open

comment:9 Changed 4 years ago by neteler

Cc: spareeth added
Keywords: AVHRR L1B added
Milestone: 1.10.0

See also (un)related tickets #2848 and #1569

Changed 4 years ago by mmetz

Attachment: l1bdataset.diff added

patch for frmts/l1b/l1bdataset.cpp to read more GCPs

comment:10 Changed 4 years ago by mmetz

The attached patch l1bdataset.diff reads all GCPs per line and attempts even spacing across X and Y. It also fixes a bug for the index of the last line. The patch has been tested with AVHRR data, using tps for georeferencing. The patch improves the georeferencing accuracy considerably.

comment:11 Changed 4 years ago by Even Rouault

MarkusM,

I'm reviewing the patch. I've scratched my head a bit before realizing that the part between line 707 and 738 (after patched applied) is now dead code, particularly the if(nDesiredGCPsPerLine < nGCPsOnThisLine && nGCPStep > 1 test that is always FALSE. So the intent is to have all available GCPs per line (51), and a sub-sampling of 51 GCP lines maximum, or every 40 (=2048 / 51) rows if nRasterYSize > nRasterXSize, or nRasterYSize if nRasterYSize < 51 ? It'll be interested (for me) to test with datasets where nRasterYSize > nRasterXSize since the sample datasets in the GDAL autotest have a very low number of lines. I'll take care of cleaning things up.

comment:12 Changed 4 years ago by neteler

Cc: mmetz added

comment:13 Changed 4 years ago by Even Rouault

Milestone: 1.10.2
Resolution: fixed
Status: assignedclosed

trunk r26710 "L1B: report correct values for GCP (#2403); report more GCPS than before; implement geolocation array"

branches/1.10 r26711 "L1B: report correct values for GCP (#2403); report more GCPS than before"

Note: See TracTickets for help on using tickets.