Opened 19 years ago
Closed 5 years 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 )
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)
Change History (14)
comment:2 by , 16 years ago
| Cc: | added; removed |
|---|---|
| Description: | modified (diff) |
comment:3 by , 15 years ago
| Cc: | added |
|---|---|
| Description: | modified (diff) |
| Keywords: | hdf added |
comment:4 by , 15 years ago
| Keywords: | aster added |
|---|---|
| Severity: | normal → major |
| Version: | unspecified → 1.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
by , 15 years ago
| Attachment: | ASTER_prdat011.hdr added |
|---|
Aster product header acquired using gdalinfo
comment:5 by , 15 years ago
| Type: | defect → enhancement |
|---|---|
| Version: | 1.5.3 → unspecified |
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 by , 15 years ago
| Priority: | high → low |
|---|
comment:7 by , 14 years ago
| Priority: | low → normal |
|---|---|
| Type: | enhancement → defect |
| Version: | unspecified → svn-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
comment:8 by , 14 years ago
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 by , 14 years ago
| Cc: | 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 by , 14 years ago
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 by , 9 years ago
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 by , 5 years ago
| Milestone: | → closed_because_of_github_migration |
|---|---|
| Resolution: | → wontfix |
| Status: | assigned → closed |
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.
