test.py

def prep_hdf(f):
    from osgeo import gdal

    print "Loading %s" % os.path.split(f.local)[1]
    hdf = gdal.Open(f.local)
    sdsdict = hdf.GetMetadata('SUBDATASETS')
    sdslist = [sdsdict[k] for k in sdsdict.keys() if '_NAME' in k]
    f.sds = [gdal.Open(n) for n in sdslist]
    return f

def extract_data(aws, d, f, layer):
    f.dataset = d.layer2name[layer]
    f.layer = f.sds[layer-1].ReadAsArray(q) #layer indices are "normal" numbers - don't start at 0
    f = process_into_npz(aws, d, f)

    aws.done = 1

    return aws, f