import sys
print(sys.version)

%pylab inline
import os, os.path
from osgeo import gdal

fn = "MOD021KM.A2004199.2140.005.2010140125955.hdf"
geofn = "MOD03.A2004199.2140.005.2010140193808.hdf"

dat = gdal.Open(fn)
geodat = gdal.Open(geofn)
print(type(dat))

metadata = dat.GetMetadata_Dict()
print(metadata['RANGEBEGINNINGDATE'], metadata['RANGEBEGINNINGTIME'])

dat.GetMetadata_Dict()['ANCILLARYINPUTTYPE'], dat.GetMetadata_Dict()['ANCILLARYINPUTPOINTER']

subdatasets = dat.GetSubDatasets()
geosubdatasets = geodat.GetSubDatasets()

subdatasets[4][1], subdatasets[7][1], geosubdatasets[8][1], geosubdatasets[9][1]

highres = gdal.Open(subdatasets[4][0], gdal.GA_ReadOnly)
midres = gdal.Open(subdatasets[7][0], gdal.GA_ReadOnly)
latsds = gdal.Open(geosubdatasets[8][0], gdal.GA_ReadOnly)
lonsds = gdal.Open(geosubdatasets[9][0], gdal.GA_ReadOnly)

print(highres.RasterCount, midres.RasterCount, latsds.RasterCount, lonsds.RasterCount)

red = highres.GetRasterBand(1)
green = midres.GetRasterBand(2)
blue = midres.GetRasterBand(1)

for color in [red, green, blue]:
    print(color.XSize, color.YSize)

imgdata = numpy.zeros([blue.YSize, blue.XSize, 3], np.uint8)
imgdata.shape

for idx, color in enumerate([red, green, blue]):
    min, max = color.ComputeRasterMinMax()
    data = np.ma.masked_greater_equal(color.ReadAsArray(), max)
    imgdata[:, :, idx] = np.multiply(255 / (max - min), data - min).astype(np.uint8)

fig = plt.figure(figsize=(18, 27))
plt.imshow(imgadata)

fig = plt.figure(figsize=(18, 27))
imgplt = plt.imshow(imgdata[0:300, 200:1200, :])

lat = latsds.GetRasterBand(1).ReadAsArray()
lon = lonsds.GetRasterBand(1).ReadAsArray()
print(np.min(lat), np.max(lat))

fig = plt.figure(figsize=(18, 12))
plt.pcolormesh(lon, lat, imgdata[..., 2], cmap='bone')