[GRASS-user] MODIS HDF conversion in TIF tool.

[Nikos Alexandris]

Well,

after searching and reading I found a script written by Markus Neteler. I
tried to change the parameters in it and get directly from sinusoidal to
HGRS87 (below the script only with the parameters changed). Initially it
looked like it works... but when comparing with some vectors data (projected
in HGRS87) a big shift appears.

After bothering Markus again about "how-to reproject MODIS surface
reflectance data into the HGRS87 (EPSG:2100)" and following his
recommendation to make it in 2 steps (1. sinusoidal to LatLong_WGS84 and 2.
(re-)reprojecting into HGRS87) it worked.

I don't understand why can't it be done in one step.

The script:

#!/bin/sh

#Warp MODIS HDF to GeoTiff - HGRS87
#Based on "modis_hdf2erdas_ll_wgs84.sh" script of Markus Neteler, 2005: warp
MODIS HDF to Erdas/Img - LatLong/WGS84

PROG="$0"

if [ $# -lt 1 -o "$1" = "-h" -o "$1" = "--help" -o "$1" = "-help" ] ; then
 echo "Script to warp all layers of MODIS HDF to GeoTiff - HGRS87"
 echo "Usage:"
 echo "      $PROG M[OD|MYD]blabla.hdf"
 exit 1
fi

FILE=$1
HDFlayerNAMES="`gdalinfo $FILE | grep SUBDATASET_ | grep _NAME | cut -d'='
-f2`"

#aarg! NDVI/EVI MOD13 contains spaces:
HDFlayerNAMES_NOSPACES="`gdalinfo $FILE | grep SUBDATASET_ | grep _NAME |
cut -d'=' -f2 | tr -s ' ' '|'`"

#MODIS does not use the "Normal Sphere (r=6370997)"!!!
# MODIS sphere, radius of 6371007.181
# http://edcdaac.usgs.gov/landdaac/tools/mrtswath/info/ReleaseNotes.pdf
#
#GGRS87 (or HGRS87 or EGSA87)
#<2100> +proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0
+ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs <>


# cs2cs -le
SPHERE="+ellps=GRS80"

#loop over all layers:
for i in ${HDFlayerNAMES_NOSPACES} ; do
        #transform back:
        INPUT="`echo ${i} | tr -s '|' ' '`"
	NEWNAME="`echo ${i} | cut -d':' -f4,5 | tr ':' '_' | tr -s '|' '_'`"
	# apply GCPs
	gdalwarp -t_srs "+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000
+y_0=0 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs $SPHERE"
"${INPUT}" $NEWNAME.tmp.tif
	# fix the LatLong reference cut to fit into world map (ulx uly lrx lry):
????
	# -projwin -180 90 180 -90 ???
	# not desired for NDVI etc ???
	gdal_translate -of HFA -a_srs "EPSG:2100" \
            $NEWNAME.tmp.tif $NEWNAME.HGRS87.tif
	rm -f $NEWNAME.tmp.tif
done
-- 
View this message in context: http://www.nabble.com/MODIS-HDF-conversion-in-TIF-tool.-tf4799564.html#a13870118
Sent from the Grass - Users mailing list archive at Nabble.com.