[FWTools] Extract to binary with meta, manipulate, then project

[Frank Warmerdam]

Simon Hook wrote:
>> This *should* exist for 1B ASTER products, so you can extract the lat and
>> lon as distinct files if you want.  But, the gdalwarp program does not
>> currently support this kind of geolocation array, so you can't ultimately
>> use them to geocode your image unless you have something else to do the
>> actual work.
> How do I extract these lat/lon arrays, I am lost as to what utility to 
> use, and how to use it. Note the embedded lat/lon arrays are a subset, 
> they are not a 1 to 1 match. This is also true with MODIS.


Simon,

I have an ASTER product here, and when I do a gdalinfo on it I get:


Driver: HDF4Image/HDF4 Dataset
Size is 2493, 2368
Coordinate System is `'
GCP Projection = PROJCS["UTM Zone 30, Northern Hemisphere",GEOGCS["WGS 84",DATUM
["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TO
WGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["E
PSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["L
at",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Me
rcator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARA
METER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_
northing",0],UNIT["Meter",1]]
GCP[  0]: Id=, Info=
           (0.5,0.5) -> (678758.000023268,4361582.99987027,0)
GCP[  1]: Id=, Info=
           (295.5,0.5) -> (688088.000033167,4361582.99987018,0)
... lots of gcps, metadata, etc...
   NUMBEROFBADPIXELS3N=0, 0, 0
   NUMBEROFBADPIXELS3B=0, 0, 0
Geolocation:
   SRS=PROJCS["UTM Zone 30,...] 
X_DATASET=HDF4_EOS:EOS_SWATH_GEOL:"ASTER_DEM20031027113839.hdf":DefaultSwath:Longitude
   X_BAND=1
 
Y_DATASET=HDF4_EOS:EOS_SWATH_GEOL:"ASTER_DEM20031027113839.hdf":DefaultSwath:Latitude
   Y_BAND=1
   PIXEL_OFFSET=0
   PIXEL_STEP=295
   LINE_OFFSET=0
   LINE_STEP=311
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 2368.0)
Upper Right ( 2493.0,    0.0)
Lower Right ( 2493.0, 2368.0)
Center      ( 1246.5, 1184.0)
Band 1 Block=2493x1 Type=Float32, ColorInterp=Gray

The GCPs shown are derived from the geolocation arrays (and further
subsampled).  But the "Geolocation:" metadata shows how to get the
actual geolocation arrays.  So based on that, I can do something like:

gdal_translate 
HDF4_EOS:EOS_SWATH_GEOL:"ASTER_DEM20031027113839.hdf":DefaultSwath:Longitude
longitude.tif


gdal_translate 
HDF4_EOS:EOS_SWATH_GEOL:"ASTER_DEM20031027113839.hdf":DefaultSwath:Latitude
latitude.tif

to grab out the longitude and latitude arrays.

The PIXEL_OFFSET/STEP and LINE_OFFSET/STEP values provide the relationship
between the geolocation array and the underlying imagery as described in
the RFC I quoted earlier.

Now this particular aster product was only a single main dataset.  I think
most 1B products actually have a bunch of subdatasets, so you would need
to do the gdalinfo on one of those to actually get the geolocation metadata.

I realize this is messy and complicated, but I would claim it is partly
because the ASTER products and HDF4 in general is messy and complicated.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org