[Gdal-dev] RFC 4 - Geolocation Arrays

[Frank Warmerdam]

Andrey Kiselev wrote:
> On Wed, Aug 02, 2006 at 10:39:16AM -0400, Frank Warmerdam wrote:
>> I have prepared RFC4 in an effort to address adding geolocation array
>> support to GDAL in a reasonably comprehensive manner.  Currently the
>> RFC is still considered to be "under development", but I thought it
>> wouldn't hurt to put it out there for feedback from the community.
> 
> Hi Frank,
> 
> Just two notes from my side.
> 
> 1. The third coordinate should be added, i.e. Z_DATASET/Z_BAND metadata
> records. Some products contain the third value (MODIS). Well, it is
> possible to have more geolocation arrays in the HDF file, so it is worth
> to add X_NAME/Y_NAME/Z_NAME field with a name of parameter
> ("Latitude"/"Longitude"/"Height"). With these names we need only one
> step to add geolocation/data maps and completely mimic the HDF-EOS
> structure :-)

Andrey,

Excellent point.  I'll update accordingly.

> 2. It should be noted that both OFFSET and STEP parameters could be
> negative in which case we have geolocation array larger than
> corresponding data array. The following is a quote from the HDF-EOS
> manual:
> 
> 	"In many cases, the size of a geolocation dimension will be
> 	different from the size of the corresponding data dimension. To
> 	take care of such occurrences, there are two pieces of
> 	information that must be supplied when defining a dimension map:
> 	the offset and the increment.  The offset tells how far along a
> 	data dimension you must travel to find the first point to have a
> 	corresponding entry along the geolocation dimension. The
> 	increment tells how many points to travel along the data
> 	dimension before the next point is found for which there is a
> 	corresponding entry along the geolocation dimension. Figure 5-3
> 	depicts a dimension map.
> 
> 	The "data skipping" method described above works quite well if
> 	there are fewer regularly spaced geolocation points than data
> 	points along a particular pair of mapped dimensions of a Swath.
> 	It is conceivable, however, that the reverse is true ­ that
> 	there are more regularly spaced geolocation points than data
> 	points. In that event, both the offset and increment should be
> 	expressed as negative values to indicate the reversed
> 	relationship. The result is shown in Figure 5-4. Note that in
> 	the reversed relationship, the offset and increment are applied
> 	to the geolocation dimension rather than the data dimension."
> 
> I haven't never seen such arrays "in wild", but everything is possible
> in the HDF world. Also negative offset is a common thing in ASTER
> products.

I am not clear why negative values need to express something special.
I have been working with the assumption that the offset could be negative
if the geolocation information extends beyond the data it is geolocating,
and that the offset and step could be floating point values.  So for
instance, if we have denser geolocation information than the data itself
the step would be between 0.0 and 1.0.

It seems to me that the "negative option" in HDF EOS can be translated
back into "forward facing" offset/step values without too much work if
we encounter it.  I'd rather not incorporate this unnecessary complexity
(of a negative option) into the gdal data model if possible.

> Translating geolocation arrays is an interesting task. The one approach
> could be a "-geolocation <name>" option of gdal_translate. When this option
> supplied, the all other options (-of,-ot,-scale,dst_dataset etc.) tells
> how the input geolocation array should be saved. The question is: do we
> need to support translating HDF-EOS datasets into the valid another
> HDF-EOS with all geolocations properly transferred? If not, that simple
> approach should work.

I have no requirement to produce HDF-EOS datasets at this point.  But it
seems that the CreateCopy() method could - in theory - be able to setup
a swath with geolocation information based on what it gets from GDAL
if needed.

> PS.
> I can take responsibility to add geolocation support to HDF4 and L1B
> drivers.

I've addressed this for the HDF4's drivers EOS_SWATH handler already
(and now this change is committed).  I introduced a new "subtype"
called EOS_SWATH_GEOL which provides "raw" access to swath fields without
all the metadata, gcps and so forth.  So the "GEOLOCATION" metadata for
a dataset can look like this:

   Metadata (GEOLOCATION):
     SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 
84",6378137,298.257223563,A
UTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM
["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHOR
ITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]
]   X_DATASET=HDF4_EOS:EOS_SWATH_GEOL:"MOD02QKM.A2005185.1820.005.2005186014821.hd
f":MODIS_SWATH_Type_L1B:Longitude
   X_BAND=1
   Y_DATASET=HDF4_EOS:EOS_SWATH_GEOL:"MOD02QKM.A2005185.1820.005.2005186014821.hd
f":MODIS_SWATH_Type_L1B:Latitude
   Y_BAND=1
   PIXEL_OFFSET=0
   PIXEL_STEP=4
   LINE_OFFSET=1
   LINE_STEP=4

I'd appreciate if you could convince yourself that the geolocation
implemention there now is correct.

I would also love to have it in the L1B driver.   I was thinking of taking
a crack at it, but I'd be most comfortable if you could do it.

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 OSGF, http://osgeo.org