Reprojecting MODIS Sea Ice Extent and IST data using GDAL

Question

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:

1. Which EPSG-code should I use in the first place?
2. 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:

1. I could not find any way to reproject pre-2011 MODIS IST data
2. I was not able to use the recommended CRS EPSG:6932/EPSG:6933 for post-2011 MODIS IST data

Answer 1

I tried the files with the GDAL utilities and see the same problem. So I then went into the HDF file with hdfview to see if the data values looked OK and where the projection information was being stored in the metadata.

SR-ORG:6932 is definitely not EPSG:6932 as it is a local area map projection for North America. The EPSG definition is at http://epsg.io/6932

gdalsrsinfo appears to be obtaining the wrong information, the file StructMetadata.0 global attribute contains the following for the southern hemisphere projection:

GROUP=GRID_2
    GridName="MOD_Grid_Seaice_4km_South"
    XDim=4501
    YDim=4501
    UpperLeftPointMtrs=(-9026314.402000,9026314.402000)
    LowerRightMtrs=(9026314.402000,-9026314.402000)
    Projection=GCTP_LAMAZ
    ProjParams=(6371228,0,0,0,0,-90000000,0,0,0,0,0,0,0)
    SphereCode=0
    GridOrigin=HDFE_GD_UL
    GROUP=Dimension
    END_GROUP=Dimension
    GROUP=DataField
        OBJECT=DataField_1
            DataFieldName="Sea_Ice_by_Reflectance_SP"
            DataType=DFNT_UINT8
            DimList=("XDim","YDim")
            CompressionType=HDFE_COMP_DEFLATE
            DeflateLevel=9
        END_OBJECT=DataField_1
        OBJECT=DataField_2
            DataFieldName="Ice_Surface_Temperature_SP"
            DataType=DFNT_UINT16
            DimList=("XDim","YDim")
            CompressionType=HDFE_COMP_DEFLATE
            DeflateLevel=9
        END_OBJECT=DataField_2
    END_GROUP=DataField
    GROUP=MergedFields
    END_GROUP=MergedFields
END_GROUP=GRID_2

A search for GCTP_LAMAZ found a document NASA GSFC (SDST_096_RevC_Final_092804.doc) that explained the ProjParms codes. From that, the following should have worked but didn't:

gdal_translate -of GTiff \
  -a_srs "+proj=laea +lat_0=90s +ellps=sphere +a=6371228.0 +b=6371228.0" \
  -a_ullr -9026314.402000 9026314.402000 9026314.402000 -9026314.402000 \
  -projwin -3000000 3000000 3000000 -3000000 \
  'HDF4_EOS:EOS_GRID:"MOD29E1D.A2009345.005.2009346083234.hdf":MOD_Grid_Seaice_4km_South:Ice_Surface_Temperature_SP' \
  MOD29E1D.A2009345.005.2009346083234.tif

The resulting GeoTIFF is in the correct position, but the values are garbled. This makes me think GDAL is also not properly handling the compression ('HDFE_COMP_DEFLATE'). Other than writing a Python script to read the data array from the HDF and write the GeoTIFF using the GDAl library, I would investigate extracting the array using one of the HDF-EOS tools at http://hdfeos.org/software/tool.php (HDF-EOS2 Dumper).

Update to the above:

The values in the data appears to be OK and understandable when the files are read using the PyHDF library (http://hdfeos.org/software/pyhdf.php). Putting this though the following Python script generates a correctly positioned GeoTIFF file with ice surface temperature values in degrees Kelvin:

#!/usr/bin/python

import os, sys
import numpy as N

# Use pyhdf library - see http://hdfeos.org/software/pyhdf.php
from pyhdf.SD import SD, SDC

import osgeo.gdal as gdal
import osgeo.osr as osr

# Target dataset
dataset_name = 'Ice_Surface_Temperature_SP'

# Open HDF file
hdffn = 'MOD29E1D.A2014345.005.2014351005302.hdf'
hdf = SD(hdffn, SDC.READ)

# List available SDS datasets.
# print hdf.datasets()

# Get dataset attributes from StructMetadata text attribute
metadata = hdf.attributes()
struct_metadata = metadata['StructMetadata.0'].split('\n')
lineno = 0
for txtline in struct_metadata:
    if txtline.find('=') > -1:
        valtxt = txtline.split('=')[1]
    if txtline.find('XDim=') > -1:
        xdim = int(valtxt)
    elif txtline.find('YDim=') > -1:
        ydim = int(valtxt)
    elif txtline.find('UpperLeftPointMtrs=') > -1:
        ulx,uly = map(float, (valtxt.strip('()')).split(',') )
    elif txtline.find('LowerRightMtrs=') > -1:
        lrx,lry = map(float, (valtxt.strip('()')).split(',') )
    elif txtline.find('Projection=') > -1:
        projection = valtxt
    elif txtline.find('ProjParams=') > -1:
        projparams = map(float, (valtxt.strip('()')).split(','))
    elif txtline.find('SphereCode=') > -1:
        sphere = int(valtxt)
    elif txtline.find('GridOrigin=') > -1:
        origin = valtxt
    elif txtline.find( ("DataFieldName=\"%s\"" % dataset_name) ) > -1:
        # This condition breaks us out of reviewing the metadata with the latest relevant values
        break
    lineno = lineno + 1
# Extract values from projparams
sphere_radius = projparams[0]
proj_latitude_origin = projparams[5] / 1000000.0

# Set the grid resolution
resx = (lrx-ulx) / xdim
resy = (uly-lry) / ydim

# Set proj4 string for projection
proj4str = ("+proj=laea +lat_0=%f +ellps=sphere +a=%f +b=%f" % \
    (proj_latitude_origin,sphere_radius,sphere_radius))
srs = osr.SpatialReference()
srs.ImportFromProj4(proj4str)

# Read Ice Surface Temperatures (IST) dataset:
#   array data, scale and offset for values, data key
ist = hdf.select(dataset_name)
ist_array = N.array(ist[:,:],dtype=N.float32)
ist_scale = ist.scale_factor
ist_offset = ist.add_offset
print ("IST scale  = %8.3f" % ist_scale)
print ("IST offset = % 8.3f" % ist_offset)
# The key attribute contains information on the range of values to be expected
# In this case, anything above 50 is considered to be a valid temperature
print ist.Key
# Apply scale and offset to get IST in degrees Kelvin, and values at 50 below
# representing certain flagged issues with data at that location
ist_array = (ist_scale * ist_array) + ist_offset
# Create an index to valid IST values
idx = N.nonzero( ist_array > 50.0 )
# Show some basic IST array information
print "IST array size:", ist_array.shape
print ("IST   minimum = %6.2f   maximum = %6.3f" % \
    (N.min(ist_array[idx]), N.max(ist_array[idx])))

# Write array to GeoTIFF file
# See https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#create-raster-from-array
gtiffdrv = gdal.GetDriverByName('Gtiff')
tiffn = ("%s.%s.tif" % (hdffn[:-4],dataset_name))
outds = gtiffdrv.Create(tiffn, xdim, ydim, 1, gdal.GDT_Float32)
outds.SetProjection(srs.ExportToWkt())
outds.SetGeoTransform((ulx, resx, 0, uly, 0, -resy))
outband = outds.GetRasterBand(1)
outband.WriteArray(ist_array)
outband.FlushCache()
outds = None

Loading the GeoTIFF file in QGIS, the coastline mask values (25.0 in the data) match up nicely to coastline vector datasets like GSHHS and Antarctic Digital Database.