MODIS HDF LAI product: ERROR 1: latitude or longitude exceeded limits in gdalwarp

[Markus Neteler]

Frank,

I was trying to warp a MODIS 15A2 HDF file today (Leaf area Index etc), but got this error:

gdalwarp -of HFA -s_srs '+proj=latlong +ellps=sphere' -t_srs EPSG:4326 'HDF4_EOS:EOS_GRID:"MOD15A2.A2003153.h18v04.004.2003171141042.hdf":MOD_Grid_MOD15A2:Lai_1km' MOD_Grid_MOD15A2_Lai_1km_LL_WGS84.img
ERROR 1: latitude or longitude exceeded limits
ERROR 1: GDALSuggestedWarpOutput() failed because the passed
transformer failed.

If I don't use +ellps=sphere it generates a map, but the coordinates are incorrect:

gdalinfo MOD_Grid_MOD15A2_Lai_1km_LL_WGS84.img
Driver: HFA/Erdas Imagine Images (.img)
Size is 1200, 1200
Coordinate System is:
GEOGCS["WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972],
        TOWGS84[0,0,0,0,0,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (0.000000,5559752.598333)
Pixel Size = (926.62543306,-926.62543306)
Corner Coordinates:
Upper Left  (       0.000, 5559752.598) (  0d 0'0.01"E,5559752d35'21474836534.00"N)
Lower Left  (       0.000, 4447802.079) (  0d 0'0.01"E,4447802d 4'17179869227.20"N)
Upper Right ( 1111950.520, 5559752.598) (1111950d31'4294967306.80"E,5559752d35'21474836534.00"N)
Lower Right ( 1111950.520, 4447802.079) (1111950d31'4294967306.80"E,4447802d 4'17179869227.20"N)
Center      (  555975.260, 5003777.338) (555975d15'35.40"E,5003777d20'17179869202.60"N)

The indicated HDF file is 5.2MB - maybe you could take a look?

thanks in advance

Markus

[neteler@…]

Frank,

would it be possible to verify this problem?
The file is available at:
http://mpa.itc.it/markus/tmp/MOD15A2.A2003153.h18v04.004.2003171141042.hdf

Best regards
Markus

[warmerdam]

Markus, 

I get a different error, that 441 of 441 points failed
to reproject.  When I do a report on the file it indicates the
final is in a sinusoidal projection so your -s_srs setting
seems inappropriate. 

warmerda@amd64[90]% gdalinfo -nomd 'HDF4_EOS:EOS_GRID:"MOD15A2.A2003153.h18v04.004.2003171141042.hdf":MOD_Grid_MOD15A2:Lai_1km' 
Driver: HDF4Image/HDF4 Dataset
Size is 1200, 1200
Coordinate System is:
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 = (926.62543306,-926.62543306)
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.339) (  7d 4'15.84"E, 45d 0'0.00"N)
Band 1 Block=1200x1 Type=Byte, ColorInterp=Gray

so far I'm not seeing anything I'd call a bug. 

[neteler@…]

Frank,

doing a re-warp to LatLong WGS84, the result is heavily shifted.

gdalwarp 'HDF4_EOS:EOS_GRID:"MOD15A2.A2003153.h18v04.004.2003171141042.hdf":MOD_Grid_MOD15A2:Lai_1km' -t_srs EPSG:4326 /tmp/modis_lai.tif

gdalinfo modis_lai.tif
Driver: GTiff/GeoTIFF
Size is 1606, 1032
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (0.000000,50.189228)
Pixel Size = (0.00968967,-0.00968967)
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.5616056,  50.1892278) ( 15d33'41.78"E, 50d11'21.22"N)
Lower Right (  15.5616056,  40.1894911) ( 15d33'41.78"E, 40d11'22.17"N)
Center      (   7.7808028,  45.1893595) (  7d46'50.89"E, 45d11'21.69"N)
Band 1 Block=1606x5 Type=Byte, ColorInterp=Gray

overlayed to:

ogrinfo -so ~/data/vmap0/europe_countries.shp europe_countries
INFO: Open of `/ssi0/ssi/neteler/data/vmap0/europe_countries.shp'
      using driver `ESRI Shapefile' successful.
Layer name: europe_countries
Geometry: Polygon
Feature Count: 3356
Extent: (-31.265751, 30.138889) - (44.834999, 80.834053)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["wgs84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
cat: Real (11.0)
vmap_na2: String (80.0)
country: String (80.0)

results in the attached map.

Only the dirty tricks as indicated above and today on the
list work:
http://lists.maptools.org/pipermail/gdal-dev/2006-November/010939.html
My script
 http://mpa.itc.it/markus/useful/modis_hdf2erdas_ll_wgs84.sh
works nicely and the result is no longer shifted (but the approach
is pretty strange).

I still suspect a problem.

Best regards,
Markus

[neteler@…]

Screenshot of MODIS overlayed to VMAP0: shift after gdalwarp

[josegomez@…]

The "dirty trick" mentioned on the mailing list <http://lists.maptools.org/pipermail/gdal-dev/2006-November/010939.html>, and in Markus Neteler's script <http://mpa.itc.it/markus/useful/modis_hdf2erdas_ll_wgs84.sh> does not work with the current version of GDAL (1.4.0, FWTools 1.2.2 on Win32). There is a shift in latitude of around 20-21km when converting from the MODIS dataset into WGS84/LatLong as per Markus' suggestion. This procedure works fine for GDAL v 1.3.2 in our system.

[1gray]

I suppose, it's not the issue with MODIS L3 files, nor with GDAL. It looks to me now that it's `v.proj' in GRASS which does the wrong thing. I'll try to post further explanations to the grass-user mailing list.

[1gray]

Replying to 1gray:

I suppose, it's not the issue with MODIS L3 files, nor with GDAL. It looks to me now that it's `v.proj' in GRASS which does the wrong thing. I'll try to post further explanations to the grass-user mailing list.

I think I've got the point. The projection description in MOD15 is wrong; the projection used is Sinusoidal on sphere, yet they feed the GCTP with WGS84 coordinates.

After reading [1] (as suggested by Maciej Sieczka) and becoming aware of the `+wktext' option, I suggest the following sample of `gdalwarp' invocation:

$ gdalwarp -of HFA \

-s_srs '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext' \ -t_srs TARGET-SRS \ 'HDF4_EOS:EOS_GRID:"fpar8.hdf":MOD_Grid_MOD15A2:Fpar_1km' \ fpar8.hfa

Note the key `+nadgrids=@null +wktext' options. In TARGET-SRS, `+ellps=sphere' becomes unnecessary (unless, of course, one wants the dataset to be warped into a projection on a sphere.)

[1] ​http://proj.maptools.org/faq.html#sphere_as_wgs84

[Kyle Shannon]

Tested in trunk against:

​http://e4ftl01.cr.usgs.gov//MOTA/MCD15A2.005/2002.07.04/MCD15A2.A2002185.h04v11.005.2007172152154.hdf

kyle@gcs-workstation:~/Downloads$ gdalwarp -of HFA -s_srs '+proj=latlong +ellps=sphere' -t_srs EPSG:4326 HDF4_EOS:EOS_GRID:"MCD15A2.A2002185.h04v11.005.2007172152154.hdf":MOD_Grid_MOD15A2:Lai_1km test.img
Creating output file that is 1200P x 1200L.
Processing input file HDF4_EOS:EOS_GRID:MCD15A2.A2002185.h04v11.005.2007172152154.hdf:MOD_Grid_MOD15A2:Lai_1km.
Using internal nodata values (e.g. 255) for image HDF4_EOS:EOS_GRID:MCD15A2.A2002185.h04v11.005.2007172152154.hdf:MOD_Grid_MOD15A2:Lai_1km.
Copying nodata values from source HDF4_EOS:EOS_GRID:MCD15A2.A2002185.h04v11.005.2007172152154.hdf:MOD_Grid_MOD15A2:Lai_1km to destination test.img.
0...10...20...30...40...50...60...70...80...90...100 - done.

Closing as stale and working in trunk. Reopen if needed.