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.