[Gdal-dev] warping MODIS data with gdal

[Stephen Hagen]

Hello,


I'm trying to reproject raw MODIS data (MOD09A1 Level 3 reflectance data
in HDF-EOS format) from its original SINUSOIDAL projection to UTM using
gdal in FWTools. I have encountered several problems. And I need help. I
put the hdf file that I am working with here: 

ftp://eos.sr.unh.edu/pub/outgoing/Frolking/Hagen/
MOD09A1.A2001289.h08v05.004.2003162200605.hdf.gz


1) gdalinfo does not recognize the projection information in the MOD09A1
hdf file:


$>gdalinfo MOD09A1.A2001289.h08v05.004.2003162200605.hdf


Driver: HDF4/Hierarchical Data Format Release 4
Size is 512, 512
Coordinate System is `'
Metadata:.....
Subdatasets:
1_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":0
SUBDATASET_1_DESC=[2400x2400] sur_refl_b01 (16-bit integer)

2_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":1
SUBDATASET_2_DESC=[2400x2400] sur_refl_b02 (16-bit integer)

3_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":2
SUBDATASET_3_DESC=[2400x2400] sur_refl_b03 (16-bit integer)

4_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":3
SUBDATASET_4_DESC=[2400x2400] sur_refl_b04 (16-bit integer)

5_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":4
SUBDATASET_5_DESC=[2400x2400] sur_refl_b05 (16-bit integer)

6_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":5
SUBDATASET_6_DESC=[2400x2400] sur_refl_b06 (16-bit integer)

7_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":6
SUBDATASET_7_DESC=[2400x2400] sur_refl_b07 (16-bit integer)

8_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":7
SUBDATASET_8_DESC=[2400x2400] sur_refl_qc_500m (32-bit unsigned integer)

9_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":8
SUBDATASET_9_DESC=[2400x2400] sur_refl_szen (16-bit integer)

10_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":9
SUBDATASET_10_DESC=[2400x2400] sur_refl_vzen (16-bit integer)

11_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":10
SUBDATASET_11_DESC=[2400x2400] sur_refl_raz (16-bit integer)

12_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":11
12_DESC=[2400x2400] sur_refl_state_500m (16-bit unsigned integer)

13_NAME=HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.hdf":12
13_DESC=[2400x2400] sur_refl_day_of_year (16-bit unsigned integer)

Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 512.0)
Upper Right ( 512.0, 0.0)
Lower Right ( 512.0, 512.0)
Center ( 256.0, 256.0)


2. gdalinfo does not recognize the sub data set projection information
either:

$>gdalinfo HDF4_SDS:UNKNOWN:"MOD09A1.A2001289.h08v05.004.2003162200605.
hdf":0

Driver: HDF4Image/HDF4 Dataset
Size is 2400, 2400
Coordinate System is `'.....
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 2400.0)
Upper Right ( 2400.0, 0.0)
Lower Right ( 2400.0, 2400.0)
Center ( 1200.0, 1200.0)

Band 1 Block=2400x1 Type=Int16, ColorInterp=Gray
Offset: -0, Scale:0.0001


3. I know the projection information from the hdf-eos header. It is
sinusoidal with a radius of 6371007.181 m. I also know the bounding box
coordinates from the header. So, I'll use gdal_translate to "force" in
projection information.


(As a side note, I found this information about the MODIS sphere in
ellipsoid.csv:


7048,GRS 1980 Authalic Sphere,6371007,9001,,6371007,0,Authalic sphere
derived from GRS 1980 ellipsoid (code 7019). (An authalic sphere is one
with a surface area equal to the surface area of the ellipsoid). 1/f is
infinite.,EPSG,EPSG,2001-06-25,,0

But I don't know how to use this information.)


$>gdal_translate -of GTiff -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0
a=6371007.181 +b=6371007.181 +units=m no_defs' -a_ullr -11119505.196667
4447802.078667 -10007554.677000 3335851.559000 -sds
MOD09A1.A2001289.h08v05.004.2003162200605.hdf
MOD09A1.A2001289.h08v05.004.2003162200605.tif


Using gdalinfo again, this appears to have worked:

$>gdalinfo MOD09A1.A2001289.h08v05.004.2003162200605.tif1

Driver: GTiff/GeoTIFF
Size is 2400, 2400
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6371007.181,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Sinusoidal"],
PARAMETER["longitude_of_center",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-11119505.196667,4447802.078667)
Pixel Size = (463.31271653,-463.31271653).....

Corner Coordinates:
Upper Left (-11119505.197, 4447802.079) (130d32'26.62"W, 40d 0'0.00"N)
Lower Left (-11119505.197, 3335851.559) (115d28'12.19"W, 30d 0'0.00"N)
Upper Right (-10007554.677, 4447802.079) (117d29'11.96"W, 40d 0'0.00"N)
Lower Right (-10007554.677, 3335851.559) (103d55'22.97"W, 30d 0'0.00"N)
Center (-10563529.937, 3891826.819) (115d58'24.91"W, 35d 0'0.00"N)

Band 1 Block=2400x1 Type=Int16, ColorInterp=Gray
Offset: -0, Scale:0.0001


4) However, when I warp this file into UTM or geographic, the resulting
map is shifted by about 19 km in the latitudinal direction.

$>gdalwarp -t_srs "EPSG:26912" -tr 463.3127 463.3127
MOD09A1.A2001289.h08v05.004.2003162200605.tif1
MOD09A1.A2001289.h08v05.004.2003162200605_utm.tif1

I've tried subsetting the data first to avoid having an image cross
several UTM zones. The results are no better.


I know I can use HEG or MRT. Those work fine and accurately. But I love 
GDAL and want to use it in every step of my giant processing stream.
MODIS data (and the file formats and projection that comes along) are
commonly used. Shouldn't gdal be able to do this? Have I missed
something?

So, my ultimate question:

WHY CAN'T GDAL ACCURATELY WARP A MODIS FILE FROM ITS ORIGINAL SINUSOIDAL
PROJECTION INTO ANOTHER PROJECTION?


OR


WHAT HAVE I DONE WRONG?


Any insights would be very much appreciated.


Best regards,

Steve