Opened 14 years ago

Closed 4 months ago

#772 closed defect (wontfix)

HDF - Georeferencing for ASTER HDF_EOS products

Reported by: warmerdam Owned by: warmerdam
Priority: normal Milestone: closed_because_of_github_migration
Component: GDAL_Raster Version: svn-trunk
Severity: major Keywords: hdf, aster
Cc: Markus Neteler, dnadeau, schut@…, Mateusz Łoskot

Description (last modified by Mateusz Łoskot)

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 (2)

ASTER_prdat011.hdr (20.0 KB) - added by kuszinger 10 years ago.
Aster product header acquired using gdalinfo
prdat011_dat.txt (26.0 KB) - added by Markus Neteler 9 years ago.
AST3A1 HDF metadata

Download all attachments as: .zip

Change History (14)

comment:1 Changed 14 years ago by warmerdam

Please ignore, just testing bugzilla.

comment:2 Changed 12 years ago by Markus Neteler

Cc: Markus Neteler added; neteler@… removed
Description: modified (diff)

comment:3 Changed 11 years ago by Mateusz Łoskot

Cc: Mateusz Łoskot added
Description: modified (diff)
Keywords: hdf added

comment:4 Changed 10 years ago by kuszinger

Keywords: aster added
Severity: normalmajor
Version: unspecified1.5.3

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 10 years ago by kuszinger

Attachment: ASTER_prdat011.hdr added

Aster product header acquired using gdalinfo

comment:5 Changed 10 years ago by kuszinger

Type: defectenhancement
Version: 1.5.3unspecified

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.

comment:6 Changed 10 years ago by kuszinger

Priority: highlow

comment:7 Changed 9 years ago by Markus Neteler

Priority: lownormal
Type: enhancementdefect
Version: unspecifiedsvn-trunk

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 9 years ago by Markus Neteler

Attachment: prdat011_dat.txt added

AST3A1 HDF metadata

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

comment:9 Changed 9 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.

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

comment:11 Changed 5 years ago by Jukka Rahkonen

This ticket needs refreshing. We would need a new "aster HDF EOS product that lacks any georeferencing" for testing because it can be that no "file lives in Frank's HDF collection" right now. And then we need at least one capable person to repeat the tests with current GDAL and report the findings here.

comment:12 Changed 4 months ago by Even Rouault

Milestone: closed_because_of_github_migration
Resolution: wontfix
Status: assignedclosed

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub?. If you believe this ticket is still valid, you may file it to https://github.com/OSGeo/gdal/issues if it is not already reported there.

Note: See TracTickets for help on using tickets.