#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

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

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

file_name = 'MYD021KM.A2008001.2040.006.2012066152535.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] 

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

norme = 0.4 # factor to increase the brightness ]0,1]

rgb = np.zeros((along_track, cross_trak,3))
for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        rgb[i,j,0] = z[i,j,0] / norme 
        rgb[i,j,1] = z[i,j,1] / norme 
        rgb[i,j,2] = z[i,j,2] / norme

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

#print 'MAX rgb', np.max( rgb[:,:,0] )

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

fig = plt.figure()

ax = fig.add_subplot(111)

img = plt.imshow(np.fliplr(rgb), 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.png", bbox_inches='tight', dpi=200)

#plt.show()

plt.close()

About / Au sujet de

Ben
Created:
6 mai 2016 21:47:08

Updated:
6 mai 2016 21:47:08

License / Licence

MIT License

Abstract / Résumé