Ticket #1239 (new defect)

Opened 2 years ago

Last modified 1 year ago

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

Reported by: neteler Assigned to: warmerdam
Priority: highest Milestone:
Component: GDAL_Raster Version: unspecified
Severity: blocker Keywords:
Cc: josegomez@gmx.net, neteler

Description (Last modified by mloskot)

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

Attachments

modis_LAI_shift_bug.png (178.0 kB) - added by neteler@itc.it on 11/28/06 12:54:34.
Screenshot of MODIS overlayed to VMAP0: shift after gdalwarp

Change History

11/15/06 13:12:20 changed by neteler@itc.it 

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

11/15/06 15:59:04 changed by 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.

11/28/06 12:53:44 changed by neteler@itc.it 

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

11/28/06 12:54:34 changed by neteler@itc.it

  • attachment modis_LAI_shift_bug.png added.

Screenshot of MODIS overlayed to VMAP0: shift after gdalwarp

02/27/07 08:57:52 changed by josegomez@gmx.net 

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.

04/03/07 16:04:24 changed by neteler

  • cc changed from josegomez@gmx.net to josegomez@gmx.net, neteler.

04/09/07 12:56:36 changed by mloskot

  • description changed.

(in reply to: ↑ description ; follow-up: ↓ 9 ) 04/16/07 05:32:53 changed by 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.

(in reply to: ↑ 8 ) 04/27/07 01:10:58 changed by 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