This is an additional exploration of the problem described at Reprojecting MODIS Ice Sea Temperature products with gdalwarp which was never answered in a satisfying way (I apologise for the length of the post).
I want to reproject Antarctic MODIS Sea Ice Extent and IST Daily L3 Global 4km EASE-Grid Day data to Antarctic Polar Stereographic using GDAL. I am interested in both AQUA (product MYD29E1D) and TERRA (product MOD29E1D) and I am interested in the full time series available.
Processing of version 6 is currently ongoing and should be completed by the end of 2015 according to https://nsidc.org/data/MYD29E1D and http://nsidc.org/data/MOD29E1D. Therefore I try to use version 5 that is made available via FTP from ftp://n5eil01u.ecs.nsidc.org/SAN/MOSA/MYD29E1D.005/ and ftp://n5eil01u.ecs.nsidc.org/SAN/MOST/MOD29E1D.005/.
These MODIS products use the rather esoteric EASE-Grid projection which exists in several versions described here: http://nsidc.org/data/ease/versions.html. Apparently there is an orignal version of EASE-Grid that was released in 1992 and a version 2.0 released in 2011. The website sited above states that six EPSG codes are available to handle EASE-Grid. For each of the two version there exist three EPSG codes relating to "North", "South" and "Global".
Original North: 3408, South: 3409, Global 3410
V2.0 North: 6931, South: 6932, Global: 6933
Which adds to the complexity is that there are two deprecated EASE-GRID EPSG codes
Deprecated North: 3973, South: 3974, Global: 3975
The question that I had now were:
- Which EPSG-code should I use in the first place?
- If there are two versions of EASE-GRID, does that mean I have to use a different EPSG-code/CRS for data created pre-2011?
GDALSRSINFO applied to the fourth band (south pole IST) results in the following definition, regardless if it is applied to a 2009 or 2014 dataset (2009 and 2014 being chosen because these data are well before and after the transition to EASE-Grid 2.0):
gdalsrsinfo 'HDF4_EOS:EOS_GRID:"MOD29E1D.A2009345.005.2009346083234.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]]
As I could not find any conclusive documentation, I tried to convert a dataset from 2009 and a dataset from 2014 using all available "Global" and "South" codes (3409, 3410, 6932, 6933, 3974, 3975). The two datasets I used are ftp://n5eil01u.ecs.nsidc.org/SAN/MOST/MOD29E1D.005/2009.12.11/ and ftp://n5eil01u.ecs.nsidc.org/SAN/MOST/MOD29E1D.005/2014.12.11/.
Regarding the software versions I use:
gdalinfo --version GDAL 1.11.2, released 2015/02/10 proj Rel. 4.9.0, 27 October 2013
I first extract the relevant band (MODIS product) and then convert it to Antarctic Polar Stereographic e.g.
gdal_translate -unscale -srcwin 500 500 3500 3500 -a_srs epsg:3409 HDF4_EOS:EOS_GRID:"MOD29E1D.A2009345.005.2009346083234.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP MOD29E1D.A2009345.005.2009346083234_tmp.tif
(The -srcwin is necessary to get rid of reprojection errors resulting from values outside the projection area)
followed by e.g.
gdalwarp -overwrite -s_srs EPSG:3409 -t_srs EPSG:3031 MOD29E1D.A2009345.005.2009346083234_tmp.tif MOD29E1D.A2009345.005.2009346083234_3409.tif
(Side remark: if I use the two commands without specifying the source CRS, the results are unusable)
The first problem that I encountered is that 6932 and 6933 are not EPSG-codes supported by the GDAL/PROJ version that I use. The errors that I get are:
ERROR 6: EPSG PCS/GCS code 6932 not found in EPSG support files. Is this a valid EPSG coordinate system? ERROR 6: EPSG PCS/GCS code 6933 not found in EPSG support files. Is this a valid EPSG coordinate system?
I then checked if the proj I use knows about EPSG:6932/6933 by using (here using the newest stable version of proj - just to make sure).
cs2cs -v +init=epsg:6932 Using from definition: init=epsg:6932 Rel. 4.9.1, 04 March 2015 <cs2cs>: projection initialization failure cause: no options found in 'init' file program abnormally terminated cs2cs -v +init=epsg:6933 Using from definition: init=epsg:6933 Rel. 4.9.1, 04 March 2015 <cs2cs>: projection initialization failure cause: no options found in 'init' file program abnormally terminated
Proj4 apparently does not know about these two EPSG codes. I hoped for a proj or WKT definition that I could use, but http://epsg.io/6932 does not give either of those. The coordinate reference systems 6932 and 6933 are therefore not usable.
I therefore applied transformations to the four remaining coordinate reference systems (3409, 3410, 3974, 3975) to the two example files. Of he eight resulting files, only one created an output with a reasonable fit: the "old" EPSG:3409 or the deprecated EPSG:3974 when applied to the 2014 dataset.
So to summarise:
- I could not find any way to reproject pre-2011 MODIS IST data
- I was not able to use the recommended CRS EPSG:6932/EPSG:6933 for post-2011 MODIS IST data