1

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

doo
  • 31
  • 2
  • 6

2 Answers2

5

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)

enter image description here

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
Rutger Kassies
  • 61,630
  • 17
  • 112
  • 97
  • **I have tried your script above but an error occurs,** _>>> import gdal_ _>>> import matplotlib.pyplot as plt_ _>>>_ _>>> ds = gdal.Open("A2014037040500.L2_LAC.SeAHABS.hdf")_ _>>> data = ds.ReadAsArray()'_ **Traceback (most recent call last): File "", line 1, in data = ds.ReadAsArray() AttributeError: 'NoneType' object has no attribute 'ReadAsArray'** – doo Apr 14 '14 at 08:04
  • That means the variable `ds` is None, which in turn means that gdal wasn't able to open it. You could have provided the wrong location or something, or you gdal build doesn't have the necessary driver. Try running gdalinfo on it, that should give a more clear error message. – Rutger Kassies Apr 15 '14 at 06:26
  • Error occurs: Traceback (most recent call last): File "C:\path\to\file\1.py", line 5, in data = ds.ReadAsArray() File "C:\Python27\lib\site-packages\osgeo\gdal.py", line 729, in ReadAsArray return gdalnumeric.DatasetReadAsArray( self, xoff, yoff, xsize, ysize, buf_obj ) File "C:\Python27\lib\site-packages\osgeo\gdal_array.py", line 175, in DatasetReadAsArray datatype = ds.GetRasterBand(1).DataType AttributeError: 'NoneType' object has no attribute 'DataType' – doo Jul 07 '14 at 02:55
  • is it possible to get geolocation data like this? How could we plot the picture with the coordinate system stored in the data instead of the pixel numbers on the axis? – Shaun Oct 17 '18 at 14:36
1

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

Emiliano
  • 19
  • 1