2

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
  • 1
    Small question, why aren't you using HEGTools for this, instead of GDAL? – Mikkel Lydholm Rasmussen Jul 23 '15 at 9:49
  • 1
    HEGTools seems to work with the older data, but GDAL only on the newer ones. It would be nice to have this work in GDAL completely. – AndreJ Jul 23 '15 at 15:19
  • Thanks for the comments. I want to use GDAL because I work in an HPC environment where all software tools have to be build from source. I therefore try to avoid specialized other tools. Besides, the NSIDC website states that PROJ and GDAL are supported which was in the first place one of the reasons why I tried this. AndreJ gave an additional reason. – Markus M. Jul 23 '15 at 20:25
2

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.

Output GeoTIFF file as displayed in QGIS with coastlines in light blue.

  • 1
    You get the same garbage without adding CRS information. So I think this is not a reprojection issue, but of data storage or compression. Strangely the newer files work almost out of the box. – AndreJ Jul 23 '15 at 18:40
  • Thank you Polarnix for your reply. I corrected my question regarding the mix-up with SR-ORG:6932 that I made. I also tried the proj params that you provided to the 2014 data which worked, but showed a significant shift. I appreciate very much that you basically confirm my findings (summarised under 1. and 2. above). I upvoted your post, but did not accept it as the solution as my problem as it is still unsolved. As AndreJ pointed out, the problem with pre-2011 data might also have to do with other problems. I will contact NSIDC and ask for a clarification regarding this problem. – Markus M. Jul 23 '15 at 21:04
  • The EASE grid webpage notes that the new data would be projected on the WGS84 ellipsoid, so you can try that instead of the sphere Polamix used. That might help with the shift problem. – AndreJ Jul 24 '15 at 3:55
  • I contacted NSIDC regarding this issue and received a response. I was informed that the data sets that I am using are in original EASE-Grid, even current data, so the EASE 2.0 info is not applicable. The used EPSG codes would be North:3408, South:3409, Global: 3410. My conclusion therefore is that the problems results from GDAL not being able to read the pre-2011 HDF files properly. I will try the approach suggested by Polarnix. I really appreciate the effort that is spent to help me solve this issue. – Markus M. Jul 26 '15 at 20:46
  • Another update: 1. although EASE-Grid 2.0 is not relevant for the v5 products MOD29E1D and MYD29E1D, it can be used with the coordinate reference system definitions provided by the EPSG registy e.g. epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::6932 2. I submitted a ticket to the GDAL devs about the problem with older MOD29E1D and MYD29E1D HDF files and will update this post when we learn something new. – Markus M. Jul 30 '15 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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