Load hdf file with GDAL and Python, get NDVI

load_hdf.py

from osgeo import gdal
import os

layer_dict = {"reli":"reliability", "qual":"Quality", "ndvi":"NDVI", "evi":"EVI"}

def print_data(layer, data):
    print data
    print "data:", layer
    print "type:", data.dtype
    print "mean:", data.mean()
    print "max :", data.max()
    print "min :", data.min()
    return

def load_hdf(fname, layer="ndvi"):
    print "\nLoading %s" % os.path.split(fname)[1]
    # convenience lookup so we can use short names
    layer_key = layer_dict.get(layer)
    hdf = gdal.Open(fname)
    sdsdict = hdf.GetMetadata('SUBDATASETS')
    sdslist =[sdsdict[k] for k in sdsdict.keys() if '_NAME' in k]
    sds = []
    for n in sdslist:
        sds.append(gdal.Open(n))
    # make sure the layer we want is in the file
    if layer_key:
        full_layer = [i for i in sdslist if layer_key in i] # returns empty if not found
        idx = sdslist.index(full_layer[0])
        data = sds[idx].ReadAsArray()
        print_data(layer, data)
        return data
    else:
        print "Layer %s not found" % layer