3

I have problems projecting MODIS Ice Sea Temperature (IST) products to EPSG:3031 (Antartic Polar Stereographic) using gdalwarp.

The product I'm working with is MOD29E1D, documented here, p33:

The global daily level-3 sea ice product containing sea ice extent and ice surface temperature (IST) is in EASE-Grid polar projection at approximately 4 km spatial resolution.

For reproducibility, the file I'm trying to process is located here and named MOD29E1D.A2010357.005.2010358085416.hdf.

I'm interested in the fourth dataset, i.e. the ice sea temperature for the South Pole. Here is the gdalinfo output for this dataset:

$ gdalinfo -stats 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2010357.005.2010358085416.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP'
Driver: HDF4Image/HDF4 Dataset
Files: MOD29E1D.A2010357.005.2010358085416.hdf
       MOD29E1D.A2010357.005.2010358085416.hdf.aux.xml
Size is 4501, 4501
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the Clarke 1866 ellipsoid",
        DATUM["Not specified (based on Clarke 1866 spheroid)",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",-5156620156.177409],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (-9026314.402000000700355,9026314.402000000700355)
Pixel Size = (4010.804000000000087,-4010.804000000000087)
  SHORTNAME=MOD29E1D
  SOUTHBOUNDINGCOORDINATE=-90.0
  SPSOPARAMETERS=none
  units=degree_Kelvin
  valid_range=22320, 31320
  VERSIONID=5
  WESTBOUNDINGCOORDINATE=-180.0
Corner Coordinates:
Upper Left  (-9026314.402, 9026314.402) (-2147483648d-2147483648'  nan"E,-2147483648d-2147483648'  nan"N)
Lower Left  (-9026314.402,-9026314.402) (-2147483648d-2147483648'  nan"E,-2147483648d-2147483648'  nan"N)
Upper Right ( 9026314.402, 9026314.402) (-2147483648d-2147483648'  nan"E,-2147483648d-2147483648'  nan"N)
Lower Right ( 9026314.402,-9026314.402) (-2147483648d-2147483648'  nan"E,-2147483648d-2147483648'  nan"N)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,Invalid angle)
Band 1 Block=4501x222 Type=UInt16, ColorInterp=Gray
  Description = Estimated sea ice surface temperature 4 km South Pole grid
  Min=0.000 Max=29668.000 
  Minimum=0.000, Maximum=29668.000, Mean=2238.933, StdDev=4362.846
  NoData Value=7
  Unit Type: degree_Kelvin
  Offset: 0,   Scale:0.01
  Metadata:
    STATISTICS_MAXIMUM=29668
    STATISTICS_MEAN=2238.9331848594
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=4362.8462148872

The output of gdalsrsinfo gives us the proj.4 string for this dataset:

$ gdalsrsinfo 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2010357.005.2010358085416.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP'
PROJ.4 : '+proj=laea +lat_0=-5156620156.177409 +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the Clarke 1866 ellipsoid",
        DATUM["Not specified (based on Clarke 1866 spheroid)",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",-5156620156.177409],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

From there, I can use gdalwarp to reproject the data layer:

$ gdalwarp  -s_srs '+proj=laea +lat_0=-5156620156.177409 +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs ' -t_srs EPSG:3031 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2010357.005.2010358085416.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP' foo.tif

However, when compared with another dataset, the resulting layer presents a displacement. In the following image, the transparent, green layer is a reference layer in the Ross Sea region, while the dark purple is the land mask of the MODIS data:

enter image description here

I assume there is a problem in the definition of the projection, any pointer would be welcome.

EDIT

The configuration on my machine:

$ proj
Rel. 4.8.0, 6 March 2012
$ gdalinfo --version
GDAL 1.10.0, released 2013/04/24

EDIT2

Even more suprises: the exact same command is failing on one of our clusters:

$ gdalwarp  -s_srs '+proj=laea +lat_0=-5156620156.177409 +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs ' -t_srs EPSG:3031 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2010357.005.2010358085416.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP' foo.tif
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: tolerance condition error
ERROR 1: Reprojection failed, err = -20, further errors will be supressed on the transform object.
Creating output file that is -2147483648P x -2147483648L.
ERROR 1: Attempt to create -2147483648x-2147483648 dataset is illegal,sizes must be larger than zero.

The configuration on this machine is:

$ proj
Rel. 4.9.0, 22 June 2013
$ gdalinfo --version
GDAL 1.11.0, released 2014/04/16
Improve this question
1
  • 2
    "latitude_of_center" looks wrong, and look at the bbox in longlat from GDAL - you'll need to sort out what it's meant to be
    – mdsumner
    Commented Sep 3, 2014 at 2:02

1 Answer 1

Reset to default
1

So, I have partially solved my question:

Turns out there are two versions of the EASE-Grid projection. For the MODIS Ice Surface Temperature data produced after the 1st of January 2011, the 2nd version is used, while the first version is used (I think) for any data acquired before.

My solution only covers the second version of the EASE-Grid:

# Clip data, apply scale and affset, and affect EPSG:3974 as data SRS
gdal_translate -unscale -srcwin 500 500 3500 3500 -a_srs epsg:3974 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2010357.005.2010358085416.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP' ist_3974.tif

# Project data to EPSG:3031
gdalwarp -s_srs EPSG:3974 -t_srs EPSG:3031 ist_3974.tif ist_3031_kelvin.tif

For bonus points, you can convert to Kelvin and sort out NoData values with gdal_calc.py:

# Convert temperatures from Kelvin to Celcius
gdal_calc.py -A ist_3031_kelvin.tif --outfile=ist_3031_celcius.tif --calc="A-273.15" --type=Float32

# Sort NoData values (here flagged to -99999)
gdal_calc.py -A ist_3031_celcius.tif --outfile=ist_3031_final.tif --calc="(A>-30.15)*A" --type=Float32 --NoDataValue=-99999

Unfortunately, we could not find any gdal_translate command to fix data collected before the 01/01/2011. Download some sample data here if you want to give it a shot!

For reference here is what the data looks like when trying to reproject from the first version of EASE-Grid with the correct EPSG code, which is 3409 according to the EASE-Grid documentation:

enter image description here

Improve this answer
1
  • 1
    The two EASE grids differ on the ellipsoid/sphere used. The old one uses a sphere (either based on Clarke 1866 as gdal reports, or International 1924 as the page you link suggests), the new uses WGS84. That explains the offset you see. Apart from that, GDAL support for the old datasets seems to be broken, as explained in the linked new question.
    – AndreJ
    Commented Jul 23, 2015 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.