[Gdal-dev] warping MODIS data with gdal

Wed Nov 22 15:20:51 EST 2006

Hi Steve;

Usually a shift in (mostly) the latitudinal direction is a sign of a
problem in converting between datums. My guess would be (and it is only
a guess) is that GDAL is having a problem with converting from a sphere
to an ellipsoid. Map projection routines (like GCTP) convert from
projection coord.s to lat, lon. But the lat is the geocentric lat in
this case. This latitude when converted to WGS84 (or what ever) must be
change to geodetic latitude. Although, 19 km seems large for just a
problem with the datum convertion, but I'd have to do some calculations
to see if that makes sense.

BTW, I just saw the LDCM Science Team has been picked. Please contact
them and tell them that the users don't want a "standard file format"
and a "standard map projection" or we'll end up like what the EOS
Science Team determined was good for us (HDF, Sinusoidal).....

Chuck

=====================================
Charles Wivell
MN DNR Forestry/Resource Assessment
413 SE 13th Street
Grand Rapids, MN 55744
218-327-4449 x 233
=====================================

>>> Stephen Hagen <steve.hagen at unh.edu> 11/22/2006 1:38 PM >>>
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

_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org 
http://lists.maptools.org/mailman/listinfo/gdal-dev