hdf driver/gdalwarp: MODIS Sinusoidal grid seems to be shifted by 0.5 pixel

[Markus Neteler]

Frank,

comparing the 'gdalwarp' results to those of the MODIS Reprojection tool (MRT),
it appears that 'gdalwarp' results seem to be shifted by 0.5 pixel.

Procedures:
(1) MODIS Sinusoidal Grid -> MRT -> LatLong/WGS84: matches other coastline data

(2) MODIS Sinusoidal Grid -> gdalwarp -> Sinusoidal -> gdalwarp ->
LatLong/WGS84: shifted to other reference data

In details (test file: MYD10A1.A2004244.h18v04.004.2004256140521.hdf)
GRingPointLatitude:  39.8198,50.0070,49.9990,39.8144
GRingPointLongitude: 0.0001,-0.0087,15.5724,13.0379

####################################################
Case (1) MODIS Reprojection tool (MRT):
  output image corners (lat/lon):
    UL:  49.999999995507 0.000000000000 
    UR:  49.999999995507 15.561593960000 
    LL:  39.995426045507 0.000000000000 
    LR:  39.995426045507 15.561593960000 

The result seen in GDAL:
gdalinfo snow.Day_Tile_Snow_Cover.tif
Size is 3212, 2065
...
Origin = (0.002422,49.997578)
Pixel Size = (0.00484483,-0.00484483)
Metadata:
  TIFFTAG_SOFTWARE=MODIS Reprojection Tool  v3.2a July 2004
Corner Coordinates:
Upper Left  (   0.0024224,  49.9975776) (  0d 0'8.72"E, 49d59'51.28"N)
Lower Left  (   0.0024224,  39.9930036) (  0d 0'8.72"E, 39d59'34.81"N)
Upper Right (  15.5640164,  49.9975776) ( 15d33'50.46"E, 49d59'51.28"N)
Lower Right (  15.5640164,  39.9930036) ( 15d33'50.46"E, 39d59'34.81"N)
Center      (   7.7832194,  44.9952906) (  7d46'59.59"E, 44d59'43.05"N)
Band 1 Block=3212x1 Type=Byte, ColorInterp=Gray

-> This looks ok. So the MRT produces correct results.


####################################################
Case (2) rewarping with 'gdalwarp':

- using gdalwarp on the HDF file, then again gdalwarp to produce LatLong:

gdalinfo snow_cover_ll2.tif
Driver: GTiff/GeoTIFF
Size is 3211, 2064
...
Origin = (0.000000,50.189228)
Pixel Size = (0.00484483,-0.00484483)
Corner Coordinates:
Upper Left  (   0.0000000,  50.1892278) (  0d 0'0.00"E, 50d11'21.22"N)
Lower Left  (   0.0000000,  40.1894911) (  0d 0'0.00"E, 40d11'22.17"N)
Upper Right (  15.5567608,  50.1892278) ( 15d33'24.34"E, 50d11'21.22"N)
Lower Right (  15.5567608,  40.1894911) ( 15d33'24.34"E, 40d11'22.17"N)
Center      (   7.7783804,  45.1893595) (  7d46'42.17"E, 45d11'21.69"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

-> Size: we have lost 1 row, 1 column
-> result of gdalwarp is shifted

####################################################

Verification of the shift, looking at North value of upper left corner:
49.999999995507-49.9975776 =  0.002422396
0.002422396 * 2 = 0.004844792

-> suggests that gdalwarp operates on MODIS Sinusoida Grid with
 0.5 pixel shift. Probably there is a mistake in confusing cell corner
 coordinates with cell center coordinates.

Summary: The results of case (2),gdalwarp, should be identical to
case (1), MRT. The MRT results are correct.

Best regards

 Markus

[warmerdam]

Andrey,

Could you take a look at this?  Do you think the geotransform setting
code in the EOS_GRID case needs to be adjusted by half a pixel? 

[neteler@…]

The MODIS scene here which I used for the experiments is here:
http://mpa.itc.it/markus/tmp/frank/
(also as tar.gz)

Thanks

 Markus

[dron]

Markus,

It seems you are right. HDF-EOS documentation is not clear on this topic. I have
made required changes.

Andrey

[neteler@…]

Andrey,

thanks for this! I have tried again, but still a major north-south
shift appears. I am not sure if it is on my side:

gdalwarp -t_srs EPSG:4326
'HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA'
MOD_Grid_Snow_500m_Day_Tile_Snow_Cover
GDAL: Auto register /usr/local/lib/gdalplugins/gdal_GRASS.so using
GDALRegister_GRASS.
HDF4Image: List of dimensions in grid Snow_Spatial_QA: YDim,XDim
HDF4Image: Grid projection: projection code: 16, zone code -1, sphere code -1
HDF4Image: Grid geolocation: top left X 0.000000, top left Y 5559752.598333, low
right X 1111950.519667, low right Y 4447802.078667, cols 2400, rows 2400
GDAL:
GDALOpen(HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA)
succeeds as HDF4Image.

GDAL: GDALOpen(MOD_Grid_Snow_500m_Day_Tile_Snow_Cover) succeeds as GTiff.

GDAL: GDALWarpKernel()::GWKNearestNoMasksByte()
Src=0,0,2400x2400 Dst=0,0,3211x2064
:0...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(MOD_Grid_Snow_500m_Day_Tile_Snow_Cover)

GDAL:
GDALClose(HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Snow_Spatial_QA)
GDAL: GDALDeregister_GTiff() called.


Observations: sphere code -1, isn't Sinusoidal on a sphere?
 http://nsidc.org/data/modis/faq.html#7

Maybe this introduces the shift when warping.

I have attached modis_aqua_snow_qgis_vmap2.jpg to illustrate the shift
(see the islands such as Elba or corsica).

Best regards

 Markus

[neteler@…]

Still shift to be observed in MODIS Grid

[dron]

I made a mistake in my previous fix. Please, try the new version. Spheroid code
-1 is normal, that means that spheroid parameters should be taken from the HDF
file, not from the predefined table.

Regards,
Andrey

[neteler@…]

Andrey,

it's still shifted... I'll run
make distclean
and try again.

Best regards,
Markus

[dron]

Markus,

Now I'm agree that there is a problem with vertical shift. The problem is not in
missed half pixel deltas, the shift is much bigger. I don't know where we have a
bug, I will dig into this issue more deeply.

Andrey

[dron]

Markus,

Are you sure that this scene has the correct georeferencing? Do you have the
same problem on other MODIS scenes? I just don't know how to solve this. I just
fetch the parameters from HDF and pass them as is.

Andrey

[neteler@…]

Andrey,

well, I was able to correctly geocode it with MRT.
Maybe they do tricks for SIN (they did once for the ISIN)?
Basically their code uses gctp (when you download MRT after
registration, you get all the source code as well).

I have put it here:
 http://mpa.itc.it/markus/tmp/MRTLinux.zip
 69MB

Maybe their reading code contains something specific?

Sorry for all these problems!

 Markus

[dron]

Hmm, I have spent even more time digging into this issue. Now I'm not sure that
we have a bug in GDAL. It is likely mismatch in ellipsoids parameters. What I
did step by step:

$ gdalwarp -t_srs "+proj=latlong +ellps=sphere"
HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover
gdal_out.tif

$ resample -p mrt.prm

where the mrt.prm is

INPUT_FILENAME = MYD10A1.A2004244.h18v04.004.2004256140521.hdf
OUTPUT_FILENAME = mrt.tif
RESAMPLING_TYPE = BILINEAR
OUTPUT_PROJECTION_TYPE = GEO
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.18 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 )

and I got the coincided results. Image sizes are different, but when I opened
them in OpenEV all coast lines become the same. Could you try this out? Note,
that you should choose the same ellipsoid types both in gdalwarp and in MRT
projection parameters.

Best regards,
Andrey

PS. Please, use the latest CVS snapshot.

[neteler@…]

Andrey,

I have just updated from CVS, it currently fails to compile:

../.libs/libgdal.so: undefined reference to `HDF4Dataset::GetDataTypeSize(long)
const'
../.libs/libgdal.so: undefined reference to `HDF4Dataset::GetDataType(long) const'
collect2: ld returned 1 exit status
make[1]: *** [ogrinfo] Fehler 1
make[1]: Leaving directory `/home/neteler/soft/gdal/ogr'
make: *** [ogr-apps] Fehler 2

even after 'make distclean'.

Do I err on my side?

Thanks,
Markus

[dron]

Ooops... please, try again.

Andrey

[neteler@…]

Andrey,

compiles, works, but still something odd.
When using gdalwarp for reprojecting from Sinusoidal
to LatLong/WGS84 I have to use +ellps=sphere (!).
Then it almost (but not perfectly) matches my
WGS84 vector data.

I assume that the hdf driver is using the wrong sphere
definition for MODIS Sinusoidal which not the PROJ4 sphere
but a "private" MODIS sphere:

GRS80 Authalic Sphere is r=6370997m,
but the MODIS sphere has instead a radius of 6371007.181.
(see http://edcdaac.usgs.gov/landdaac/tools/mrtswath/info/ReleaseNotes.pdf).

I don't know how to verify this but would like to ask you
to do it...

Once MODIS SIN is assigned to the correct sphere I am sure that
it finally will work.

Best regards

 Markus

[neteler@…]

Andrey,

compiles, works, but still something odd.
When using gdalwarp for reprojecting from Sinusoidal
to LatLong/WGS84 I have to use +ellps=sphere (!).
Then it almost (but not perfectly) matches my
WGS84 vector data.

I assume that the hdf driver is using the wrong sphere
definition for MODIS Sinusoidal which not the PROJ4 sphere
but a "private" MODIS sphere:

GRS80 Authalic Sphere is r=6370997m,
but the MODIS sphere has instead a radius of 6371007.181.
(see http://edcdaac.usgs.gov/landdaac/tools/mrtswath/info/ReleaseNotes.pdf).

I don't know how to verify this but would like to ask you
to do it...

Once MODIS SIN is assigned to the correct sphere I am sure that
it finally will work.

Best regards

 Markus

[dron]

Markus,

HDF driver uses the _right_ sphere.

$ gdalinfo
HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover

...

PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the custom spheroid",
        DATUM["Not specified (based on custom spheroid)",
            SPHEROID["Custom spheroid",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["Meter",1]]
...

Note that spheroid is the same as you are mentioned (6371007.181). This is
fetched from the HDF file in agreement with the GCTP reference. I just can't
find nothing criminal in the georeferencing code. Have you tried the
reprojection scheme I described previously? I think that the source of problem
in your case is wrong spheroid assumptions for different layers.

Regards,
Andrey

[neteler@…]

> From Andrey Kiselev  2005-04-19 10:23  [reply] -------
> Markus,
> HDF driver uses the _right_ sphere.

Sorry for the confusion. Yes, the sphere is right. However,
rewarped maps result shifted. No idea, why...

Below the procedure (sorry, long, but hopefully useful):

# from MODIS HDF to normal sphere:
gdalwarp -t_srs "+proj=latlong +ellps=sphere" \
 
HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover
\
  gdal_normal_sphere.tif

gdalinfo gdal_normal_sphere.tif | head -25
Driver: GTiff/GeoTIFF
Size is 3211, 2064
Coordinate System is:
GEOGCS["Normal Sphere (r=6370997)",
    DATUM["unknown",
        SPHEROID["unnamed",6370997,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (0.000000,50.000000)
Pixel Size = (0.00484489,-0.00484489)
Corner Coordinates:
Upper Left  (   0.0000000,  50.0000000) (  0d 0'0.01"E, 50d 0'0.00"N)
Lower Left  (   0.0000000,  40.0001409) (  0d 0'0.01"E, 40d 0'0.51"N)
Upper Right (  15.5569513,  50.0000000) ( 15d33'25.02"E, 50d 0'0.00"N)
Lower Right (  15.5569513,  40.0001409) ( 15d33'25.02"E, 40d 0'0.51"N)
Center      (   7.7784756,  45.0000705) (  7d46'42.51"E, 45d 0'0.25"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

# From normal sphere to WGS84
gdalwarp -t_srs EPSG:4326 gdal_normal_sphere.tif gdal_wgs84_LL.tif

gdalinfo gdal_wgs84_LL.tif
Driver: GTiff/GeoTIFF
Size is 3211, 2064
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (0.000000,50.189228)
Pixel Size = (0.00484485,-0.00484485)
Corner Coordinates:
Upper Left  (   0.0000000,  50.1892281) (  0d 0'0.01"E, 50d11'21.22"N)
Lower Left  (   0.0000000,  40.1894659) (  0d 0'0.01"E, 40d11'22.08"N)
Upper Right (  15.5568006,  50.1892281) ( 15d33'24.48"E, 50d11'21.22"N)
Lower Right (  15.5568006,  40.1894659) ( 15d33'24.48"E, 40d11'22.08"N)
Center      (   7.7784003,  45.1893470) (  7d46'42.24"E, 45d11'21.65"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

#Comparing to Vector map:
ogrinfo -summary polbnda_europe_wgs84.shp polbnda_europe_wgs84 | head -18
INFO: Open of `polbnda_europe_wgs84.shp'
using driver `ESRI Shapefile' successful.

Layer name: polbnda_europe_wgs84
Geometry: Polygon
Feature Count: 4378
Extent: (-31.265751, 17.224695) - (46.725971, 80.834053)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

... should match, right? Unfortunately it doesn't, see attachment
modis_snow_vmap0_wgs84_gdal.png

No idea what's wrong.

########
This is what MRT does:
  input image corners (lat/lon):
    UL:  50.00 0.00 
    UR:  50.00 15.56 
    LL:  40.00 0.00 
    LR:  40.00 13.05 

  input image spatial subset corners (lat/lon):
    UL:  50.00 0.00 
    UR:  50.00 15.56 
    LL:  40.00 0.00 
    LR:  40.00 15.56 

  SINUSOIDAL PROJECTION PARAMETERS:
   Radius of Sphere:     6371007.181000 meters
   Longitude of Center:     0.000000 degrees
   False Easting:      0.000000 meters 
   False Northing:     0.000000 meters 

  Output image info
  -----------------
  output image corners (lat/lon):
    UL:  49.999999995507 0.000000000000 
    UR:  49.999999995507 15.562266608930 
    LL:  39.994423054211 0.000000000000 
    LR:  39.994423054211 15.562266608930 
  output image corners (X-Y projection units):
    UL:  0.000000000000 49.999999995507 
    UR:  15.562266608930 49.999999995507 
    LL:  0.000000000000 39.994423054211 
    LR:  15.562266608930 39.994423054211 
    
     band               type lines smpls pixsiz      min    max   fill
  1) Day_Tile_Snow_Cover  UINT8  1698  2641 0.0059      0    255    255

The MRT result is shown in the second attachment, which looks good:
modis_snow_vmap0_wgs84_MRT.png

File details:
gdalinfo modis_snow_vmap0_wgs84_MRT.Day_Tile_Snow_Cover.tif
Driver: GTiff/GeoTIFF
Size is 2641, 1698
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (0.002946,49.997054)
Pixel Size = (0.00589257,-0.00589257)
Metadata:
  TIFFTAG_SOFTWARE=MODIS Reprojection Tool  v3.2a July 2004
Corner Coordinates:
Upper Left  (   0.0029463,  49.9970537) (  0d 0'10.61"E, 49d59'49.39"N)
Lower Left  (   0.0029463,  39.9914768) (  0d 0'10.61"E, 39d59'29.32"N)
Upper Right (  15.5652129,  49.9970537) ( 15d33'54.77"E, 49d59'49.39"N)
Lower Right (  15.5652129,  39.9914768) ( 15d33'54.77"E, 39d59'29.32"N)
Center      (   7.7840796,  44.9942652) (  7d47'2.69"E, 44d59'39.35"N)
Band 1 Block=2641x1 Type=Byte, ColorInterp=Gray

Find also attached the gdalinfo output of GDALwarped file and MRTwarped
file.

Overall, I don't want to stay on your nerves but just help to find the
problem.

Best regards

 Markus
 

[neteler@…]

MODIS gdalwarp result vs VMAP0 (both WGS84)

[neteler@…]

MODIS MRT program result vs VMAP0 (both WGS84)

[neteler@…]

gdalinfo output of gdalwarp result

[neteler@…]

gdalinfo output of MRT result

[dron]

Markus,

What about file with parameters for transforation with MRT?

> Overall, I don't want to stay on your nerves but just help to find the problem.

No problm, I want to solve this poblm too.

I see the shift when I use different ellipsoids and long/lat "projection". As
soon as I have transformed both images to the same ellipsoids the shift goes away.

Andrey

[neteler@…]

MODIS MRT parameter file

[neteler@…]

MODIS MRT log file

[neteler@…]

> I see the shift when I use different ellipsoids and long/lat "projection". As
> soon as I have transformed both images to the same ellipsoids the shift goes away.

Yes, but according to the gdalinfo output the ellipsoids are identical
Do you see a mistake in my procedure as posted earlier today?

Thanks

 Markus

[dron]

I do not see any mistake, also I heve read through the RT sources and found that
they are calculate bounding coordinates identically to me (anyway, the half
pixel error cannot be a reason of so huge shift).

Andrey

[dron]

Well, one more experiment. gdalinfo reports about source MODIS dataset the
following:

Size is 2400, 2400
PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the custom spheroid",
        DATUM["Not specified (based on custom spheroid)",
            SPHEROID["Custom spheroid",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["Meter",1]]
Origin = (0.000000,5559752.598333)
Pixel Size = (463.31271653,-463.31271653)
....
Corner Coordinates:
Upper Left  (       0.000, 5559752.598) (  0d 0'0.01"E, 50d 0'0.00"N)
Lower Left  (       0.000, 4447802.079) (  0d 0'0.01"E, 40d 0'0.00"N)
Upper Right ( 1111950.520, 5559752.598) ( 15d33'26.06"E, 50d 0'0.00"N)
Lower Right ( 1111950.520, 4447802.079) ( 13d 3'14.66"E, 40d 0'0.00"N)
Center      (  555975.260, 5003777.338) (  7d 4'15.84"E, 45d 0'0.00"N)
Band 1 Block=2400x1 Type=Byte, ColorInterp=Gray

I do

$ gdalwarp -t_srs EPSG:4326
HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover
gdal_wgs84.tif

And now I have

Size is 3211, 2064
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (0.000000,50.189228)
Pixel Size = (0.00484483,-0.00484483)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (   0.0000000,  50.1892278) (  0d 0'0.01"E, 50d11'21.22"N)
Lower Left  (   0.0000000,  40.1894911) (  0d 0'0.01"E, 40d11'22.17"N)
Upper Right (  15.5567608,  50.1892278) ( 15d33'24.34"E, 50d11'21.22"N)
Lower Right (  15.5567608,  40.1894911) ( 15d33'24.34"E, 40d11'22.17"N)
Center      (   7.7783804,  45.1893595) (  7d46'42.17"E, 45d11'21.69"N)
Band 1 Block=3211x2 Type=Byte, ColorInterp=Gray

Note, that the warped image is 11' shifted to north! Thatis exactly where we
have a problem. I don't know why this happens, probably this is an error in the
datum shifting code which arises in cases of transformation to unprojected
scene. I will do more researches in this direction

Andrey

[dron]

Loooks like this bug should be reassigned against PROJ.4 (or I just
misunderstand something). Let's try the datum shifting to geographic coordinates:

$ cat coords.lst
771180 4594849
$ cs2cs +proj=sinu +ellps=sphere +to +proj=latlong +ellps=sphere coords.lst 
9d14'5.33"E     41d19'20.961"N 0.000
$ cs2cs +proj=sinu +ellps=sphere +to +proj=latlong +ellps=WGS84 coords.lst 
9d14'5.33"E     41d30'48.042"N 2210.028

Well, there is 11' shift we are walking around. I'm not too knowledgeable in
PROJ.4 code, but I shall try to find the error.

Andrey

[warmerdam]

Andrey mentioned in IRC that he expected the lat/long location to be the same
on both ellipsoids (sphere and WGS84).  

Internally PROJ.4 is converting the location on the sphere into a geocentric
location, and then convert it back to lat/long on WGS84.  Because their shapes
are quite different the angles are also noticably different (as well as the big
shift in Z). 

[dron]

Markus,

I have discussed this issue with Frank and what we are decided for now:

1. PROJ.4 is correct in its transformation, it properly shifts coordinates
between different ellipsoids.
2. MRT is not completely correct. It should shift the image too.
3. If you want use both PROJ.4 and MRT results simultaneously you should not
warp the image between datums based on different ellipsoids, but just rewrite
the projection definition, i.e.

$ gdalwarp -t_srs "+proj=latlong +ellps=sphere"
HDF4_EOS:EOS_GRID:"MYD10A1.A2004244.h18v04.004.2004256140521.hdf":MOD_Grid_Snow_500m:Day_Tile_Snow_Cover
gdal_sphere.tif
$ gdal_translate -a_srs EPSG:4326 gdal_sphere.tif gdal_wgs84.tif

gdal_wgs84.tif now should be compatible with the MRT output.

Regards,
Andrey