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. _______________________________________________ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user