#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from scipy.misc import bytescale

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import pprint

# References:
# https://earthdata.nasa.gov/sites/default/files/field/document/MODIS_True_Color.pdf
# http://www.idlcoyote.com/ip_tips/brightmodis.html

#----------------------------------------------------------------------------------------#
# inputs

file_name = 'MYD021KM.A2013189.1350.006.2013190155358.hdf'

file = SD(file_name, SDC.READ)

#----------------------------------------------------------------------------------------#

selected_sds = file.select('EV_250_Aggr1km_RefSB')

selected_sds_attributes = selected_sds.attributes()

for key, value in selected_sds_attributes.iteritems():
	#print key, value
	if key == 'reflectance_scales':
		reflectance_scales_250_Aggr1km_RefSB = np.asarray(value)		
	if key == 'reflectance_offsets':
		reflectance_offsets_250_Aggr1km_RefSB = np.asarray(value)	

sds_data_250_Aggr1km_RefSB = selected_sds.get()

print 'sds_data shape', sds_data_250_Aggr1km_RefSB.shape
print reflectance_scales_250_Aggr1km_RefSB.shape

#----------------------------------------------------------------------------------------#

selected_sds = file.select('EV_500_Aggr1km_RefSB')

selected_sds_attributes = selected_sds.attributes()

for key, value in selected_sds_attributes.iteritems():
	if key == 'reflectance_scales':
		reflectance_scales_500_Aggr1km_RefSB = np.asarray(value)	
	if key == 'reflectance_offsets':
		reflectance_offsets_500_Aggr1km_RefSB = np.asarray(value)	

sds_data_500_Aggr1km_RefSB = selected_sds.get()

print reflectance_scales_500_Aggr1km_RefSB.shape

#----------------------------------------------------------------------------------------#

data_shape = sds_data_250_Aggr1km_RefSB.shape

along_track = data_shape[1]
cross_trak = data_shape[2]

z = np.zeros((along_track, cross_trak,3))

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,0] = ( sds_data_250_Aggr1km_RefSB[0,i,j] - \
        reflectance_offsets_250_Aggr1km_RefSB[0] ) * \
        reflectance_scales_250_Aggr1km_RefSB[0] 

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,1] = ( sds_data_500_Aggr1km_RefSB[1,i,j] - \
        reflectance_offsets_500_Aggr1km_RefSB[1] ) * \
        reflectance_scales_500_Aggr1km_RefSB[1]  

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,2] = ( sds_data_500_Aggr1km_RefSB[0,i,j] - \
        reflectance_offsets_500_Aggr1km_RefSB[0] ) * \
        reflectance_scales_500_Aggr1km_RefSB[0] 

z[ z > 1 ] = 1.0
z[ z < 0 ] = 0.0

#----------------------------------------------------------------------------------------#
# Color enhancement

x = np.array([0,  30,  60, 120, 190, 255], dtype=np.uint8)
y = np.array([0, 110, 160, 210, 240, 255], dtype=np.uint8)

def scale_image(image, x, y):
    scaled = np.zeros((along_track, cross_trak), dtype=np.uint8)
    for i in range(len(x)-1):
        x1 = x[i]
        x2 = x[i+1]
        y1 = y[i]
        y2 = y[i+1]
        m = (y2 - y1) / float(x2 - x1)
        b = y2 - (m *x2)
        mask = ((image >= x1) & (image < x2))
        scaled = scaled + mask * np.asarray(m * image + b, dtype=np.uint8)

    mask = image >= x2
    scaled = scaled + (mask * 255)

    return scaled

z_color_enh = np.zeros((along_track, cross_trak,3), dtype=np.uint8)
z_color_enh[:,:,0] = scale_image(bytescale(z[:,:,0]), x, y)
z_color_enh[:,:,1] = scale_image(bytescale(z[:,:,1]), x, y)
z_color_enh[:,:,2] = scale_image(bytescale(z[:,:,2]), x, y)

#----------------------------------------------------------------------------------------#
# plot image using matplotlib

fig = plt.figure()

ax = fig.add_subplot(111)

img = plt.imshow(np.fliplr(z_color_enh), interpolation='nearest', origin='lower')

ticks_x_pos = [0,250,500,750,1000,1250]
ticks_x_label = ['0','250','500', '750', '1000','1250']

ticks_y_pos = [0,250,500,750,1000,1250,1500,1750,2000]
ticks_y_label = ['0','250','500', '750', '1000','1250','1500','1750','2000']

plt.xticks(ticks_x_pos, ticks_x_label, rotation=0, fontsize=8 )
plt.yticks(ticks_y_pos, ticks_y_label, rotation=0, fontsize=8 )

plt.title('MODIS RGB Image', fontsize=10)

plt.savefig("modis_granule_rgb_color_enh.png", bbox_inches='tight', dpi=200)

#plt.show()

plt.close()

About / Au sujet de

Ben
Created:
1 juin 2016 22:53:07

Updated:
1 juin 2016 22:53:07

License / Licence

MIT License

Abstract / Résumé