[GRASS-user] HDF-EOS swath to grid (GeoTiff) conversion

[Ivan Shmakov]

>>>>> Nikos Alexandris <nikos.alexandris at felis.uni-freiburg.de> writes:

[...]

 >>> you can use gdal_translate to add (per scanline) GCPs to the metadata
 >>> before running gdalwarp.  gdal_translate -gcp pixel line easting
 >>> northing [elevation]

 >> Indeed, it works, at least for the AIRS data.  (Although I've
 >> had to upgrade GDAL 1.3.2, as supplied with Debian 4.0 r1, to
 >> 1.4.3.)

 >> For those who are interested, the script I've used looks roughly
 >> as follows:

 > You think your script would be applicable to Surface Reflectance
 > products (L2) ?

	Do you mean MODIS L2 data?  Well, I'll assume it for now...

 > Sorry for bothering again, but since I've limited experience with
 > scripts... I ask before I I try to modify it for my need.

	I'll try to comment on some of the points where the script is
	tied closely to the AIRS L2 data format.

 >> --cut--
 >> #!/bin/bash
 >> file=AIRS.2007.07.09.079.L2.RetStd.v4.0.9.102.D07190164053.hdf
 >> file_gcp=airs-2007-07-09-gran-79-l2-retstd-tsurfstd-gcp+0.5.tiff
 >> file_out=airs-2007-07-09-gran-79-l2-retstd-tsurfstd-gcp+0.5-laea-10km.tiff

 >> ## FIXME: specify ``no projection'' instead?
 >> source_srs='+proj=latlong'

	First of all, this seems to specify the projection system for
	the GCPs.  Thus, it's essential.  No FIXME: needed.

 >> target_srs='+proj=laea +lat_0=55.0 +lon_0=90'

 >> ## FIXME: more clean error handling
 >> set -x -e

 >> long_lat_to_gcp () {
 >>     gawk '$1 != -9999 && $2 != -9999 {
 >>               print "-gcp", (NR - 1) % 30 + .5, int ((NR - 1) / 30) + .5, $1, $2;

	Here, the `-gcp MINOR MAJOR LONGITUDE LATITUDE' lines are made.
	It's assumed that the data is stored as an 30x45 (30 pixels
	along the minor dimension) array, i. e.:

a [1, 1] ... a [1, 30] a [2, 1] ... a [2, 30] ... a [45, 30]

	These are to be tailored for the data to be processed (e. g.,
	1354 x (10 * scans) for 1 km MODIS L2 swath.)

 >>           }'
 >> }

 >> gdal_translate \
 >>     $(paste <(hdfdump --text "$file" {Long,Lat}itude) | long_lat_to_gcp) \

	Extraction of the geolocation arrays is done with `hdfdump' [1],
	yet another ``HDF to plain-something'' conversion tool.  It
	could probably be done with `hdp' or GRASS (through the ``x, y''
	location) as well.

	Now, the `-gcp ...' lines are put into the command line for
	`gdal_translate'.  Command line arguments arrays are subject to
	OS- (and configuration-)dependent limits, and it's assumed that
	there's enough room for that many options.

	This assumption could easily become false for MODIS data on most
	OS, with some 1000000 points (or more.)  Either some points are
	to be discarded, or the GCP setting is to be applied in multiple
	passes.  I haven't tried doing either.

 >>     "HDF4_EOS:EOS_SWATH:\"$file\":L2_Standard_atmospheric&surface_product:TSurfStd" \

	This line have to be tailored according to the way GDAL ``sees''
	the file.  One could use `gdalinfo' on the file in question to
	find the pattern one needs.

 >>     "$file_gcp"
 >> gdalwarp \
 >>     -tps -s_srs "$source_srs" -t_srs "$target_srs" \
 >>     -tr 1e4 1e4 \

	Since MODIS data has finer resolution than AIRS, the destination
	resolution should be changed as well.  E. g., `-tr 1e3 1e3'
	could probably be used for 1 km MODIS L2 data.

 >>     "$file_gcp" "$file_out"
 >> --cut--

 >> Is there any chance that GDAL will obtain ``GCPs'' from the
 >> HDF-EOS file directly?

	Of course, this would eliminate all the aforementioned issues,
	and the script itself.

 >> [...]

 >>>> The one more problem with the MODIS L2 data is that it consists of
 >>>> overlapping scans -- MODIS could ``see'' some parts of the Earth in
 >>>> consequent scans.

	This is one more question I've had to check.  May be the
	conversion has to be done scan-wise -- each scan is to be
	selected and georeferenced (with `gdal_translate' and, e. g.,
	the above script), then reprojected with `gdalwarp'.  The
	resulting reprojected strips are to be somehow merged after.

[...]

[1] http://theory.asu.ru/~ivan/devel/hdfdump/

PS.  I feel that the gdal-dev mailing list is a better place to discuss
	this kind of problems.  Could we move there?