Unable to read HDF4 dataset using pyhdf (pyhdf.err

Question

I am trying to read a HDF4 file (https://www.dropbox.com/s/5d40ukfsu0yupwl/MOD13A2.A2016001.h23v05.006.2016029070140.hdf?dl=0).

import os
import numpy as np
from pyhdf.SD import SD, SDC

# Open file.
FILE_NAME = 'MOD13A2.A2016001.h23v05.006.2016029070140.hdf'
hdf = SD(FILE_NAME, SDC.READ)

# List available SDS datasets.
print (hdf.datasets())

# Read dataset.
DATAFIELD_NAME="1_km_16_days_NDVI"
data2D = hdf.select(DATAFIELD_NAME)
data = data2D[:,:] 

When I executed this script, I am getting following error: Traceback (most recent call last): File "Test.py", line 15, in data2D = hdf.select(DATAFIELD_NAME) File "C:\Python35\lib\site-packages\pyhdf\SD.py", line 1599, in select raise HDF4Error("select: non-existent dataset") pyhdf.error.HDF4Error: select: non-existent dataset

I have used similar python code for reading other HDF4 files and it works well. But I am unable to understand problem in this case.

Answer 1

If you look carefully at the output of print(hdf.datasets()) you would notice the following line:

 ...
 '1 km 16 days NDVI': (('YDim:MODIS_Grid_16DAY_1km_VI',
                        'XDim:MODIS_Grid_16DAY_1km_VI'),
                       (1200, 1200),
                       22,
                       0),
 ...

The key is the name of the dataset. Note that the name uses spaces to separate words, not the underscores as in your example.

If you replace the DATAFIELD_NAME with DATAFIELD_NAME='1 km 16 days NDVI'. That is replace underscores with spaces in the dataset name, your code would work without the error.

This is enough to access the data within the tile, but geolocation requires more work. At the moment you already can visualize the data with

from matplotlib import pyplot as plt
plt.imshow(data)
plt.colorbar()
plt.show()

What remains to be done is

• scaling to scientific units
• geolocation

As an aside, do use the pprint module to output large dictionaries with complex structure, i.e. replace print(hdf.datasets()) with

import pprint
...
pprint.pprint(hdf.datasets())