Ticket #772 (assigned defect)

Opened 8 years ago

Last modified 3 years ago

HDF - Georeferencing for ASTER HDF_EOS products

Reported by: warmerdam Owned by: warmerdam
Priority: normal Milestone:
Component: GDAL_Raster Version: svn-trunk
Severity: major Keywords: hdf, aster
Cc: neteler, dnadeau, schut@…, mloskot

Description (last modified by mloskot) (diff)

Some aster HDF EOS products lack any georeferencing.

For example AST_09_003011420030252520000000.hdf1

The file level metadata contains lots of potentially useful metadata but none of it is reflected on the HDF_EOS:EOS_SWATH subdatasets. Metadata includes: 

 ASTERMapProjection=Universal Transverse Mercator
  UPPERLEFT=-2.309363, 113.508649
  UPPERRIGHT=-2.406042, 114.172772
  LOWERLEFT=-2.872922, 113.427496
  LOWERRIGHT=-2.969329, 114.092004
  SCENECENTER=-2.639387, 113.800183
  MAPORIENTATIONANGLE=8.351934
  BANDSCALEFACTORS=0.000000, 0.000000, 0.000000, 0.000000, 0.217400, 0.069600, 0
.062500, 0.059700, 0.041700, 0.031800, 0.000000, 0.000000, 0.000000, 0.000000, 0
.000000
  PROJECTIONPARAMETERS1=6378137.000000, 6356752.300000, 0.999600, 0.000000, 1.93
7315, 0.000000, 500000.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
 0.000000, 0.000000, 0.000000, 0.000000
  UTMZONECODE1=49

File lives in Frank's HDF collection.

Attachments

ASTER_prdat011.hdr Download (20.0 KB) - added by kuszinger 4 years ago.
Aster product header acquired using gdalinfo
prdat011_dat.txt Download (26.0 KB) - added by neteler 3 years ago.
AST3A1 HDF metadata

Change History

Changed 8 years ago by warmerdam 

Please ignore, just testing bugzilla.

Changed 5 years ago by neteler

  • cc neteler added; neteler@… removed
  • description modified (diff)

Changed 5 years ago by mloskot

  • cc mloskot added
  • keywords hdf added
  • description modified (diff)

Changed 4 years ago by kuszinger

  • keywords hdf, aster added; hdf removed
  • version changed from unspecified to 1.5.3
  • severity changed from normal to major

Now years later I faced the same problem. Using gdalwarp everything seems fine, but when looking the converted image, it seems that it is ~10 kms away.

gdlainfo header of original aster product is attached. 

	E	N
UL	104,368197	14,087236
UR	105,075064	14,088052
LL	104,369926	13,442014
LR	105,074859	13,442791

CC	104,719631	13,765897  (image center)

in utm wgs84 zone 48N it is: 

proj transform UTM WGS84 48N	
431795,691705	1557465,741863
508103,125870	1557465,748215
431795,622932	1486104,404283
508103,181400	1486104,383859

469692,227378 1521853,570515 center..

GDALWARP command is here: 

gdalwarp -t_srs '+proj=utm +no_defs +zone=48 +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000' HDF4_EOS:EOS_SWATH:"/home/kuszi/grass/etc/ASTER_prdat011.dat":VNIR_Swath:ImageData1 a5.tif

After running gdalinfo on the resulting tif, it shows that converted boundaries ~10km to the south (eg. 1 476 502 instead of 1 486 104) 

UL	431761,36115951	1547454,27111329
UR	508114,65089312	1547454,27111329
LL	431761,36115951	1476502,41927774
LR	508114,65089312	1476502,41927774

GDAL version: 1.5.3 OS: Linux tpgentoo 2.6.27-gentoo-r8 #1 SMP PREEMPT Tue Mar 10 17:11:25 CET 2009 i686 Intel(R) Pentium(R) M processor 1.50GHz GenuineIntel? GNU/Linux

Changed 4 years ago by kuszinger

Aster product header acquired using gdalinfo

Changed 4 years ago by kuszinger

  • version changed from 1.5.3 to unspecified
  • type changed from defect to enhancement

Status change - enhancement may be necessary:

An option would be nice to use corner data instead of GCPs in a HDF file. For rothorectified data products it is OK.

Reason:

The error in my case is in HDF-EOS aster file. Meanwhile header and corner information are correct, the ground control points are on a bad location (~10km slide as posted earlier).

It is a contradiction in the file itseld. Gdal tools uses GCPs in priority over header LL, LR, UL, UR corner data.

Workout: georeference a 'traditional' way if one faces a problem like this (and your data product is orthorectified). Always check the consistency of in-HDF GCPs and corner coordinates. In case of in-file contradiction georeference 'manually' using corner data. It worked for me and everything is 'in its correct place' now.

Changed 4 years ago by kuszinger

  • priority changed from high to low

Changed 3 years ago by neteler

  • priority changed from low to normal
  • version changed from unspecified to svn-trunk
  • type changed from enhancement to defect

I take liberty to change this to defect: The metadata (attached) define UTM/WGS84 project but the result of gdalwarp on AST3A1 HDF is UTM on GRS80: 

gdalinfo ast3a01_EOS_SWATH_prdat011.dat.tif
Files: ast3a01_EOS_SWATH_prdat011.dat.tif
Size is 5565, 5068        
Coordinate System is:    
PROJCS["UTM Zone 54, Northern Hemisphere", GEOGCS["Unknown datum based upon the GRS 1980 ellipsoid", 
DATUM["Not_specified_based_on_GRS_1980_spheroid", 
SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",141],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (257052.541733831982128,3920585.257126705255359)
Pixel Size = (0.000000003883098,-0.000000003883098)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  257052.542, 3920585.257) (138d19'29.56"E, 35d23'56.61"N)
Lower Left  (  257052.542, 3920585.257) (138d19'29.56"E, 35d23'56.61"N)
Upper Right (  257052.542, 3920585.257) (138d19'29.56"E, 35d23'56.61"N)
Lower Right (  257052.542, 3920585.257) (138d19'29.56"E, 35d23'56.61"N)
Center      (  257052.542, 3920585.257) (138d19'29.56"E, 35d23'56.61"N)
Band 1 Block=5565x1 Type=Byte, ColorInterp=Gray

GRS80 and WGS84 have subtle differences in the Semi-minor axis b (hence Inverse flattening).

Could this be fixed? thanks,

Markus

Changed 3 years ago by neteler

AST3A1 HDF metadata

Changed 3 years ago by neteler

Not sure if that's related, but all Image data have 0 as min and max: 

for i in *.tif ; do echo $i ; gdalinfo -mm $i | grep Computed ; done
ast3a01_DEMZS_Swath_DEMZSData.tif                                                         
    Computed Min/Max=-9999.000,3135.000
ast3a01_DEMZT_Swath_DEMZTData.tif      
    Computed Min/Max=-9999.000,3135.000
ast3a01_DEMZV_Swath_DEMFlagData.tif    
    Computed Min/Max=0.000,176.000
ast3a01_DEMZV_Swath_DEMZVData.tif 
    Computed Min/Max=-9999.000,3135.000
ast3a01_SWIR_Swath_ImageData4.tif      
    Computed Min/Max=0.000,0.000
ast3a01_SWIR_Swath_ImageData5.tif
    Computed Min/Max=0.000,0.000
ast3a01_SWIR_Swath_ImageData6.tif
    Computed Min/Max=0.000,0.000
ast3a01_SWIR_Swath_ImageData7.tif
    Computed Min/Max=0.000,0.000
ast3a01_SWIR_Swath_ImageData8.tif
    Computed Min/Max=0.000,0.000
ast3a01_SWIR_Swath_ImageData9.tif
    Computed Min/Max=0.000,0.000
ast3a01_TIR_Swath_ImageData10.tif
    Computed Min/Max=0.000,0.000
ast3a01_TIR_Swath_ImageData11.tif
    Computed Min/Max=0.000,0.000
ast3a01_TIR_Swath_ImageData12.tif
    Computed Min/Max=0.000,0.000
ast3a01_TIR_Swath_ImageData13.tif
    Computed Min/Max=0.000,0.000
ast3a01_TIR_Swath_ImageData14.tif
    Computed Min/Max=0.000,0.000
ast3a01_VNIR_Swath_ImageData1.tif
    Computed Min/Max=0.000,0.000
ast3a01_VNIR_Swath_ImageData2.tif
    Computed Min/Max=0.000,0.000
ast3a01_VNIR_Swath_ImageData3B.tif
    Computed Min/Max=0.000,0.000
ast3a01_VNIR_Swath_ImageData3N.tif
    Computed Min/Max=0.000,0.000

Only the included DEM data look ok. The file can be made available (personally).

Using svn trunk from a few days ago, on a 64 bit Linux box.

Markus

Changed 3 years ago by warmerdam

  • cc dnadeau added

Denis,

Do you have any interest in digging into this issue? If not, let me know, and I'll put it on Chaitanya's queue.

Changed 3 years ago by neteler

Here more detailed tests: 

# Look at original HDF file
gdalinfo -mm prdat011.dat                     
GDAL: GDALOpen(prdat011.dat, this=0x1b5ea10) succeeds as HDF4.          
Driver: HDF4/Hierarchical Data Format Release 4                         
Files: prdat011.dat                  
Size is 512, 512                     
Coordinate System is `'              
Metadata:                            
  HDFEOSVersion=HDFEOS_V2.12         
  IDOFASTERGDSDATAGRANULE=AST3A1 0801060139460912160052                 
  RECEIVINGCENTER=EDOS               
  PROCESSINGCENTER=ASTER-GDS         
...
  ELLIPSOIDANDDATUM4=WGS84, WGS84                
  DATUMPARA4=6378137.000000, 298.257224, 0.000000, 0.000000, 0.000000               
  RESMETHOD4=NN                                  
  MPMETHOD4=UTM                                  
  PROJECTIONPARAMETERS4=0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000                
  UTMZONECODE4=54
...


# wrong ellipsoid, but min/max correct:
[neteler@north aster3A1]$ gdalinfo -mm HDF4_EOS:EOS_SWATH:"prdat011.dat":VNIR_Swath:ImageData1 | tail                  
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
OGRCT: Target: +proj=utm +zone=54 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
OGRCT: Target: +proj=utm +zone=54 +ellps=GRS80 +units=m +no_defs
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
...
  PROCESSINGLEVELID=3                            
  MAPPROJECTIONNAME=Universal Transverse Mercator
Corner Coordinates:                              
Upper Left  (    0.0,    0.0)                    
Lower Left  (    0.0, 4902.0)                    
Upper Right ( 5112.0,    0.0)                    
Lower Right ( 5112.0, 4902.0)                    
Center      ( 2556.0, 2451.0)                    
Band 1 Block=5112x1 Type=Byte, ColorInterp=Gray  
    Computed Min/Max=0.000,255.000               


# translate to GeoTIFF: wrong ellipsoid applied
[neteler@north aster3A1]$ gdal_translate HDF4_EOS:EOS_SWATH:"prdat011.dat":VNIR_Swath:ImageData1 img1.tif
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs                    
OGRCT: Target: +proj=utm +zone=54 +ellps=GRS80 +units=m +no_defs      
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs                    
OGRCT: Target: +proj=utm +zone=54 +ellps=GRS80 +units=m +no_defs      


# Geotiff: wrong ellipsoid but still correct min/max:
[neteler@north aster3A1]$ gdalinfo -mm img1.tif
GDAL: GDALOpen(img1.tif, this=0x1bdc9f0) succeeds as GTiff.
Driver: GTiff/GeoTIFF
GDAL: GDALDefaultOverviews::OverviewScan()                 
Files: img1.tif
       img1.tif.aux.xml
Size is 5112, 4902 
Coordinate System is `'
GCP Projection =
PROJCS["UTM Zone 54, Northern Hemisphere",                 
    GEOGCS["Unknown datum based upon the GRS 1980 ellipsoid",
        DATUM["Not_specified_based_on_GRS_1980_spheroid",    
            SPHEROID["GRS 1980",6378137,298.2572221010002,   
                AUTHORITY["EPSG","7019"]]],                  
        PRIMEM["Greenwich",0],                               
        UNIT["degree",0.0174532925199433]],                  
    PROJECTION["Transverse_Mercator"],                       
    PARAMETER["latitude_of_origin",0],                       
    PARAMETER["central_meridian",141],                       
    PARAMETER["scale_factor",0.9996],                        
    PARAMETER["false_easting",500000],                       
    PARAMETER["false_northing",0],                           
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]                           
GCP[  0]: Id=1, Info=
          (2138893713408.5,0.5) -> (219962.042896839,3917467.00769924,0)
GCP[  1]: Id=2, Info=
          (4277787427314.5,0.5) -> (227450.101730825,3917467.00305552,0)
GCP[  2]: Id=3, Info=
          (6416681141220.5,0.5) -> (234938.162253144,3917466.9985378,0) 
GCP[  3]: Id=4, Info=
          (8555574855126.5,0.5) -> (242426.224417487,3917466.9941461,0) 
....
  MAPPROJECTIONNAME=Universal Transverse Mercator
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4902.0)
Upper Right ( 5112.0,    0.0)
Lower Right ( 5112.0, 4902.0)
Center      ( 2556.0, 2451.0)
Band 1 Block=5112x1 Type=Byte, ColorInterp=Gray
    Computed Min/Max=0.000,255.000
GDAL: GDALClose(img1.tif, this=0xe5d9f0)
GDAL: GDALDeregister_GTiff() called.


# warp map:
[neteler@north aster3A1]$ gdalwarp img1.tif img1b.tif
GDAL: GDALOpen(img1.tif, this=0x24c99f0) succeeds as GTiff.
WARP: Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE
GDAL: GDALClose(img1.tif, this=0x24c99f0)                  
Creating output file that is 5565P x 5068L.                
GDAL: GDALDriver::Create(GTiff,img1b.tif,5565,5068,1,Byte,(nil))
GDAL: GDALOpen(img1.tif, this=0x24ccf10) succeeds as GTiff.     
Processing input file img1.tif.                                 
0...10...20...30...40...50...60...70...80...90...100 - done.    
GDAL: GDALClose(img1.tif, this=0x24ccf10)                       
GDAL: GDALClose(img1b.tif, this=0x24ead60)                      
GDAL: GDALDeregister_GTiff() called.                            


# wrong ellipsoid and now wrong min/max
[neteler@north aster3A1]$ gdalinfo -mm img1b.tif
GDAL: GDALOpen(img1b.tif, this=0x1159ac0) succeeds as GTiff.
Driver: GTiff/GeoTIFF
GDAL: GDALDefaultOverviews::OverviewScan()                  
Files: img1b.tif
Size is 5565, 5068
Coordinate System is:
PROJCS["UTM Zone 54, Northern Hemisphere",                  
    GEOGCS["Unknown datum based upon the GRS 1980 ellipsoid",
        DATUM["Not_specified_based_on_GRS_1980_spheroid",    
            SPHEROID["GRS 1980",6378137,298.2572221010002,   
                AUTHORITY["EPSG","7019"]]],                  
        PRIMEM["Greenwich",0],                               
        UNIT["degree",0.0174532925199433]],                  
    PROJECTION["Transverse_Mercator"],
...
Origin = (257052.540050168550806,3920585.282846763730049)
Pixel Size = (0.000000003883098,-0.000000003883098)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
OGRCT: Source: +proj=utm +zone=54 +ellps=GRS80 +units=m +no_defs
OGRCT: Target: +proj=longlat +ellps=GRS80 +no_defs
Corner Coordinates:
Upper Left  (  257052.540, 3920585.283) (138d19'29.56"E, 35d23'56.61"N)
Lower Left  (  257052.540, 3920585.283) (138d19'29.56"E, 35d23'56.61"N)
Upper Right (  257052.540, 3920585.283) (138d19'29.56"E, 35d23'56.61"N)
Lower Right (  257052.540, 3920585.283) (138d19'29.56"E, 35d23'56.61"N)
Center      (  257052.540, 3920585.283) (138d19'29.56"E, 35d23'56.61"N)
Band 1 Block=5565x1 Type=Byte, ColorInterp=Gray
    Computed Min/Max=0.000,0.000
GDAL: GDALClose(img1b.tif, this=0x1159ac0)
GDAL: GDALDeregister_GTiff() called.

At this point the resulting map is all zero.

Note: See TracTickets for help on using tickets.