Read/Open a modis aqua .hdf file and display/plot the output in gdal and matplotlib

Question

I have tried and search on how to solve this but still can't find a way on how to do read and plot this in gdal and matplotlib from a given Modis Aqua .hdf file. Any help is much appreciated. By the way am using Python 2.7.5 in Windows 7. The filename is A2014037040000.L2_LAC.SeAHABS.hdf.Among the Geophysical Datas of the hdf file I will only be using the chlor_a.

Update:

Here is the link of the sample file.

A2014037040500.L2_LAC.SeAHABS.hdf

Answer 1

The trick with HDF's is that most of the time you need a specific subdataset. If you use GDAL you need to open the HDF pointing directly to that subdataset:

import gdal
import matplotlib.pyplot as plt

ds = gdal.Open('HDF4_SDS:UNKNOWN:"MOD021KM.A2013048.0750.hdf":6')
data = ds.ReadAsArray()
ds = None

fig, ax = plt.subplots(figsize=(6,6))

ax.imshow(data[0,:,:], cmap=plt.cm.Greys, vmin=1000, vmax=6000)

You can also open the 'main' HDF file and inspect the subdatasets, and go from there:

# open the main HDF
ds = gdal.Open('MOD021KM.A2013048.0750.hdf')

# get the path for a specific subdataset
subds = [sd for sd, descr in ds.GetSubDatasets() if descr.endswith('EV_250_Aggr1km_RefSB (16-bit unsigned integer)')][0]

# open and read it like normal
dssub = gdal.Open(subds)
data = dssub.ReadAsArray()
dssub = None

ds = None

Answer 2

You should try setting datatype for the MODIS dataset. I guess it is 16 bits unsigned

ds= gdal.Open(hdfpath) data = ds.GetRasterBand(N).ReadAsArray().astype(numpy.uint16)

N is the band number for your data of interest. You can try open it with QGIS or ENVI to see the structure of the HDF file.

Remember that the bands starts at 1 and not as 0. First band is 1.

Hope it helps