MODIS data / Sinusoidal projection

[Andreas Hilboll]

I'm working with MODIS data, which come in HDF-EOS2 format. I can read
these files using GDAL, and the projection information looks like this:

In [52]: ds.GetProjection()
Out[52]: '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]]'

In [53]: ds.GetGeoTransform()
Out[53]:
(5559752.598333, 926.625433055833, 0.0,
8895604.157333, 0.0, -926.6254330549995)

The meaning of the GDAL geotransform are as follows:

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

Can I represent this projection in Iris? iris.coord_systems doesn't
have a 'sinusoidal' member ...

Cheers,
Andreas.

[Andrew Dawson]

HI Andreas

This is not currently implemented. It needs implementing in cartopy first, and then an Iris coordinate system can follow. There is a ticket open for adding this to cartopy, for reference: https://github.com/SciTools/cartopy/issues/650

Andrew

[Luis García Carreras]

Hi Andreas,

I've had to read modis data with Iris, and my solution was to use the pyproj module to generate 2D lat/lon coordinates from Xgeo, Ygeo as you've defined them, and add these as auxiliary coordinates to the cube containing the modis data. If you have a cube of modis data 'modis':

import pyproj

   proj = pyproj.Proj(modis.attributes['PROJECTION'])
   (mx, my) = numpy.meshgrid(Xgeo, Ygeo)
   (lon, lat) = proj(mx, my, inverse=True)

   longitude = iris.coords.AuxCoord(lon, standard_name='longitude', units='degrees')
   latitude = iris.coords.AuxCoord(lat, standard_name='latitude', units='degrees')

   modis.add_aux_coord(longitude, data_dims=[1,2])
   modis.add_aux_coord(latitude, data_dims=[1,2])