Genrate grid information file from MODIS HDFEOS data

Question

Is there a way to generate grid information (lat-lon) from the MODIS MCD19A2 files in python?. The file is downloaded from Link to the data file .In MATLAB it can be done using the following block code

import matlab.io.hdf4.*
import matlab.io.hdfeos.*

% Open the HDF-EOS2 Grid file.
FILE_NAME='MCD19A2.A2010010.h25v06.006.2018047103710.hdf';
file_id = gd.open(FILE_NAME, 'rdonly');

% Read data from a data field.
GRID_NAME='grid1km';
grid_id = gd.attach(file_id, GRID_NAME);

DATAFIELD_NAME='Optical_Depth_055';
[data, lat, lon] = gd.readField(grid_id, DATAFIELD_NAME, [], [], []);

In short, I am looking for a pyhdf/python equivalent for gd.readField of MATLAB package

Answer 1

HDF-EOS Tools and Information Center Help was so nice to provide a script to deal with grid definition. This can be found here. In case the link is not working, here is the code:

"""
Copyright (C) 2014-2019 The HDF Group
Copyright (C) 2014 John Evans

This example code illustrates how to access and visualize an LP DAAC MCD19A2
v6 HDF-EOS2 Sinusoidal Grid file in Python.

If you have any questions, suggestions, or comments on this example, please use
the HDF-EOS Forum (http://hdfeos.org/forums).  If you would like to see an
example of any other NASA HDF/HDF-EOS data product that is not listed in the
HDF-EOS Comprehensive Examples page (http://hdfeos.org/zoo), feel free to
contact us at [email protected] or post it at the HDF-EOS Forum
(http://hdfeos.org/forums).

Usage:  save this script and run

    $python MCD19A2.A2010010.h25v06.006.2018047103710.hdf.py


Tested under: Python 3.7.3 :: Anaconda custom (64-bit)
Last updated: 2019-09-20
"""
import os
import re
import pyproj

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

from pyhdf.SD import SD, SDC
from mpl_toolkits.basemap import Basemap

FILE_NAME = 'MCD19A2.A2010010.h25v06.006.2018047103710.hdf'
DATAFIELD_NAME = 'Optical_Depth_055'
hdf = SD(FILE_NAME, SDC.READ)

# Read dataset.
data3D = hdf.select(DATAFIELD_NAME)
data = data3D[1,:,:].astype(np.double)

# Read attributes.
attrs = data3D.attributes(full=1)
lna=attrs["long_name"]
long_name = lna[0]
vra=attrs["valid_range"]
valid_range = vra[0]
fva=attrs["_FillValue"]
_FillValue = fva[0]
sfa=attrs["scale_factor"]
scale_factor = sfa[0]        
ua=attrs["unit"]
units = ua[0]
aoa=attrs["add_offset"]
add_offset = aoa[0]

# Apply the attributes to the data.
invalid = np.logical_or(data < valid_range[0], data > valid_range[1])
invalid = np.logical_or(invalid, data == _FillValue)
data[invalid] = np.nan
data = (data - add_offset) * scale_factor
data = np.ma.masked_array(data, np.isnan(data))

# Construct the grid.  The needed information is in a global attribute
# called 'StructMetadata.0'.  Use regular expressions to tease out the
# extents of the grid.
fattrs = hdf.attributes(full=1)
ga = fattrs["StructMetadata.0"]
gridmeta = ga[0]
ul_regex = re.compile(r'''UpperLeftPointMtrs=\(
                          (?P<upper_left_x>[+-]?\d+\.\d+)
                          ,
                          (?P<upper_left_y>[+-]?\d+\.\d+)
                          \)''', re.VERBOSE)

match = ul_regex.search(gridmeta)
x0 = np.float(match.group('upper_left_x'))
y0 = np.float(match.group('upper_left_y'))

lr_regex = re.compile(r'''LowerRightMtrs=\(
                          (?P<lower_right_x>[+-]?\d+\.\d+)
                          ,
                          (?P<lower_right_y>[+-]?\d+\.\d+)
                          \)''', re.VERBOSE)
match = lr_regex.search(gridmeta)
x1 = np.float(match.group('lower_right_x'))
y1 = np.float(match.group('lower_right_y'))

nx, ny = data.shape
x = np.linspace(x0, x1, nx)
y = np.linspace(y0, y1, ny)
xv, yv = np.meshgrid(x, y)

sinu = pyproj.Proj("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext")
wgs84 = pyproj.Proj("+init=EPSG:4326") 
lon, lat= pyproj.transform(sinu, wgs84, xv, yv)


# There is a wrap-around issue to deal with, as some of the grid extends
# eastward over the international dateline.  Adjust the longitude to avoid
# a smearing effect.
lon[lon < 0] += 360

m = Basemap(projection='cyl', resolution='l',
            llcrnrlat=np.min(lat), urcrnrlat = np.max(lat),
            llcrnrlon=np.min(lon), urcrnrlon = np.max(lon))                
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(np.floor(np.min(lat)), np.ceil(np.max(lat)), 5),
                labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(np.floor(np.min(lon)), np.ceil(np.max(lon)), 5),
                labels=[0, 0, 0, 1])

# Subset data if you don't see any plot due to limited memory.
# m.pcolormesh(lon[::2,::2], lat[::2,::2], data[::2,::2], latlon=True)
m.pcolormesh(lon, lat, data, latlon=True)


cb = m.colorbar()
cb.set_label(units)

basename = os.path.basename(FILE_NAME)
plt.title('{0}\n{1}'.format(basename, long_name))
fig = plt.gcf()
pngfile = "{0}.py.png".format(basename)
fig.savefig(pngfile)

Answer 2

Thanks for sharing the python file, but there is a problem:

the latitude and longitude read by python and matlab are not always the same. For example, the dimension of h23v04 read by python starts at 50°N, while the value read by matlab starts at 49.9958°N