Plot CALIOP CLay vertical profile product with python and matplotlib (ex: cloud phase)

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC # To read a HDF file

import numpy as np
import matplotlib.pyplot as plt
import sys

#----------------------------------------------------------------------------------------#
# bits_stripping Function

def bits_stripping(bit_start,bit_count,value):
    bitmask=pow(2,bit_start+bit_count)-1
    return np.right_shift(np.bitwise_and(value,bitmask),bit_start)

#----------------------------------------------------------------------------------------#
# Read HDF Files

caliop_file_name = 'CAL_LID_L2_01kmCLay-ValStage1-V3-01.2008-07-03T08-17-58ZD.hdf'
#caliop_file_name = 'CAL_LID_L2_05kmCLay-Prov-V3-01.2008-07-03T08-17-58ZD.hdf'
#caliop_file_name = 'CAL_LID_L2_333mCLay-ValStage1-V3-01.2008-07-03T08-17-58ZD.hdf'

caliop = SD(caliop_file_name, SDC.READ)

caliop_latitude = caliop.select('Latitude')
caliop_longitude = caliop.select('Longitude')
caliop_Profile_Time_UTC = caliop.select('Profile_UTC_Time')
caliop_Number_Layers_Found = caliop.select('Number_Layers_Found')
caliop_Layer_Top_Altitude = caliop.select('Layer_Top_Altitude')
caliop_Layer_Base_Altitude = caliop.select('Layer_Base_Altitude')
caliop_Feature_Classification_Flags = caliop.select('Feature_Classification_Flags')

caliop_latitude_data = caliop_latitude.get()
caliop_longitude_data = caliop_longitude.get()
caliop_Profile_Time_UTC_data = caliop_Profile_Time_UTC.get()
caliop_Number_Layers_Found_data = caliop_Number_Layers_Found.get()
caliop_Layer_Top_Altitude_data = caliop_Layer_Top_Altitude.get()
caliop_Layer_Base_Altitude_data = caliop_Layer_Base_Altitude.get()
caliop_Feature_Classification_Flags_data = caliop_Feature_Classification_Flags.get()

print "Read hdf done !"

#print caliop_Profile_Time_UTC_data[0,0]
#print caliop_Profile_Time_UTC_data.shape

shape = caliop_Profile_Time_UTC_data.shape
#print 'caliop_Profile_Time_UTC_data: ', shape

#print caliop_Profile_Time_UTC_data[0], caliop_Profile_Time_UTC_data[shape[0]-1]
#print caliop_Profile_Time_UTC_data[0] / 86400.0
#print caliop_Profile_Time_UTC_data[shape[0]-1] / 86400.0

#----------------------------------------------------------------------------------------#
# CALIOP granule start (cgs) time ?

tmp =  str( caliop_Profile_Time_UTC_data[0,0] )
#print tmp, float( tmp[5:] )

time_UTC = float( tmp[5:] )
cgs_fod = time_UTC
cgs_hour = int( time_UTC * 24 )
cgs_minute = int( time_UTC * 24 * 60 ) - cgs_hour * 60
cgs_second = int( time_UTC * 24 * 60 * 60 ) - cgs_hour * 60 * 60 - cgs_minute * 60

print 'CALIOP granule start time: ', cgs_hour, cgs_minute, cgs_second

#----------------------------------------------------------------------------------------#
# CALIOP granule end (cge) time ?

tmp = str( caliop_Profile_Time_UTC_data[shape[0]-1,0] )
#print tmp, float( tmp[5:] )

time_UTC = float( tmp[5:] )
cge_fod = time_UTC
cge_hour = int( time_UTC * 24 )
cge_minute = int( time_UTC * 24 * 60 ) - cge_hour * 60
cge_second = int( time_UTC * 24 * 60 * 60 ) - cge_hour * 60 * 60 - cge_minute * 60

print 'CALIOP granule end time : ', cge_hour, cge_minute, cge_second

#----------------------------------------------------------------------------------------#
# User Inputs 

start_hour = 8
start_minute = 35
start_second = 0

start_time = start_hour * 60 * 60 + start_minute * 60 + start_second
start_time_fod = start_time / 86400.0

end_hour = 8
end_minute = 40
end_second = 0

end_time = end_hour * 60 * 60 + end_minute * 60 + end_second
end_time_fod = end_time / 86400.0

#print start_time, start_time_fod
#print end_time, end_time_fod

#----------------------------------------------------------------------------------------#
# Check

if start_time_fod > cgs_fod and end_time_fod < cge_fod:
	print 'correct inputs !'
else:
	print 'incorrect inputs !'
	sys.exit() 

#----------------------------------------------------------------------------------------#
# Main loop (time filter)

fig = plt.figure(figsize=(12,6))

#fig.patch.set_facecolor('0.65')

fig = fig.add_subplot(111)

lat_min = 999
lat_max = -999

for i in np.arange( shape[0] ):
	tmp =  str( caliop_Profile_Time_UTC_data[i,0] )
	time_UTC = float( tmp[5:] )
	if time_UTC > start_time_fod and time_UTC < end_time_fod:
		hour = int( time_UTC * 24 )
		minute = int( time_UTC * 24 * 60 ) - hour * 60
		second = int( time_UTC * 24 * 60 * 60 ) - hour * 60 * 60 - minute * 60		
		#print hour, minute, second, caliop_Number_Layers_Found_data[i,0], \
		#      caliop_latitude_data[i,0], caliop_longitude_data[i,0] 
		if caliop_latitude_data[i,0] < lat_min : 
			lat_min = caliop_latitude_data[i,0]
		if caliop_latitude_data[i,0] > lat_max : 
			lat_max = caliop_latitude_data[i,0]			
		if caliop_Number_Layers_Found_data[i,0] >= 1:
			for j in np.arange( caliop_Number_Layers_Found_data[i,0] ):
				feature_type = bits_stripping(0,3,
				               caliop_Feature_Classification_Flags_data[i,j])
				if feature_type == 2: 
					bar_xi = caliop_latitude_data[i-1,0] + \
					abs( caliop_latitude_data[i,0] - caliop_latitude_data[i-1,0] ) / 2.0
					bar_width = abs( caliop_latitude_data[i,0] - \
					                 caliop_latitude_data[i-1,0] )
					bar_base = caliop_Layer_Base_Altitude_data[i,j] 
					bar_height = caliop_Layer_Top_Altitude_data[i,j] - bar_base
					caliop_layer_phase = bits_stripping(5,2,
					caliop_Feature_Classification_Flags_data[i,j])
					bar_color = '0.55'
					bar_edgecolor = '0.55'
					if caliop_layer_phase == 1:
						bar_color = '#C9CDEC'
						bar_edgecolor = '#C9CDEC'
					elif caliop_layer_phase == 2:
						bar_color = '#0D0575'
						bar_edgecolor = '#0D0575'
					if caliop_layer_phase == 1 or caliop_layer_phase == 2:	
						plt.bar(left=bar_xi, height=bar_height, width=bar_width, 
						bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

print 'lat_min, lat_max: ', lat_min, lat_max

title = 'CALIOP Cloud Thermodynamic Phase'
output_file_name = 'caliop_' 
if '333mCLay' in caliop_file_name:
	title = title + '\n (333m Cloud Layer Product)'
	output_file_name = output_file_name + '333CLay_cloud_phase'
if '01kmCLay' in caliop_file_name:
	title = title + '\n (1km Cloud Layer Product)'
	output_file_name = output_file_name + '01CLay_cloud_phase'
if '05kmCLay' in caliop_file_name: 
	title = title + '\n (5km Cloud Layer Product)'
	output_file_name = output_file_name + '05CLay_cloud_phase'

fig.set_title(title)
fig.set_xlabel('Latitude')
fig.set_ylabel('Altitude (km)')

plt.xlim((lat_min, lat_max))
plt.ylim((0, 30))

plt.grid()

#fig.patch.set_facecolor('0.55')

plt.savefig(output_file_name, dpi=80 )
plt.show()

Plot MODIS granule RGB image (orthographic projection) using python and basemap

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from scipy.misc import bytescale
from mpl_toolkits.basemap import Basemap, cm 
from scipy import interpolate
from scipy.interpolate import griddata
from matplotlib import colors

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 = 'MYD021KM.A2013189.1350.006.2013190155358.hdf'
file_name_myd03 = 'MYD03.A2013189.1350.006.2013190152754.hdf'

file_myd021km = SD(file_name_myd021km, SDC.READ)

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

selected_sds = file_myd021km.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()

data_shape = sds_data_250_Aggr1km_RefSB.shape

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

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

selected_sds = file_myd021km.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()

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

file_myd03 = SD(file_name_myd03, SDC.READ)

selected_sds = file_myd03.select('Latitude')
myd03_lat = selected_sds.get()

selected_sds = file_myd03.select('Longitude')
myd03_long = selected_sds.get()

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

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

#----------------------------------------------------------------------------------------#
# Rough estimation of latitude and longitude at granule center (long_0, lat_0)

lat_min = myd03_lat[0,0]
lat_max = myd03_lat[along_track-1,cross_trak-1]

lat_0 = lat_min + (lat_max - lat_min) / 2.

long_min = min(myd03_long[0,0],myd03_long[along_track-1,cross_trak-1])
long_max = max(myd03_long[0,0],myd03_long[along_track-1,cross_trak-1])

lon_0 = long_min + (long_max - long_min) / 2.

#----------------------------------------------------------------------------------------#
# Orthographic Map Projection

fig = plt.figure()

ax = fig.add_subplot(111)

ax.patch.set_facecolor((0.75,0.75,0.75))

m1 = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution=None)

xpt0, ypt0 = m1(lon_0,lat_0) 

xpt1, ypt1 = m1(myd03_long[0,0],myd03_lat[0,0]) 
xpt2, ypt2 = m1(myd03_long[0,cross_trak-1],myd03_lat[0,cross_trak-1]) 
xpt3, ypt3 = m1(myd03_long[along_track-1,cross_trak-1], \
                myd03_lat[along_track-1,cross_trak-1])
xpt4, ypt4 = m1(myd03_long[along_track-1,0],myd03_lat[along_track-1,0])

llx = min(xpt1,xpt2,xpt3,xpt4) - xpt0  # lower left
lly = min(ypt1,ypt2,ypt3,ypt4) - ypt0

urx = max(xpt1,xpt2,xpt3,xpt4) - xpt0  # upper right
ury = max(ypt1,ypt2,ypt3,ypt4) - ypt0

m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution='l',\
    llcrnrx=llx,llcrnry=lly,urcrnrx=urx,urcrnry=ury)

x_igrid, y_igrid = m(myd03_long,myd03_lat)

x_igrid = x_igrid - xpt0
y_igrid = y_igrid - ypt0

z_igrid_01 = np.zeros((along_track, cross_trak))
z_igrid_02 = np.zeros((along_track, cross_trak))
z_igrid_03 = np.zeros((along_track, cross_trak))

for i in np.arange(2030):
    for j in np.arange(1354): 
        z_igrid_01[i,j] = z[i,j,0]
        z_igrid_02[i,j] = z[i,j,1]
        z_igrid_03[i,j] = z[i,j,2]

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()
z_igrid_01 = z_igrid_01.ravel()
z_igrid_02 = z_igrid_02.ravel()
z_igrid_03 = z_igrid_03.ravel()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T
xi, yi = np.mgrid[llx:urx:1000j, lly:ury:1000j]

z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='cubic')
z_02 = griddata(xy1_igrid, z_igrid_02, (xi, yi), method='cubic')
z_03 = griddata(xy1_igrid, z_igrid_03, (xi, yi), method='cubic')

rgb_projected = np.zeros((1000, 1000,3))
for i in np.arange(1000):
    for j in np.arange(1000): 
        rgb_projected[i,j,0] = z_01[i,j]
        rgb_projected[i,j,1] = z_02[i,j]
        rgb_projected[i,j,2] = z_03[i,j]

#rgb_projected[ z > 1 ] = 1.0
#rgb_projected[ z < 0 ] = 0.0
whereAreNaNs = np.isnan(rgb_projected);
rgb_projected[whereAreNaNs] = 0.75;

img = m.imshow(np.rot90(np.fliplr(rgb_projected)), origin='lower')

m.drawcoastlines()

m.drawparallels(np.arange(-90.,120.,5.), color='k', labels=[True,False,False,False])
m.drawmeridians(np.arange(0.,420.,5.), color='k', labels=[False,False,False,True])

ax.set_xlabel("", fontsize=10)
ax.set_ylabel("", fontsize=10)       

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

plt.savefig("modis_granule_rgb_orth_projection.png", dpi=100)
plt.show()

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

Plot gridded data on a global map using python and basemap

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm

import matplotlib.pyplot as plt

import matplotlib as mpl

import numpy as np

data = np.loadtxt("gridded_data.txt")

m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,\
            llcrnrlon=-180,urcrnrlon=180,resolution='c')

m.drawcoastlines()

m.drawlsmask(land_color='white',ocean_color='white',lakes=True)

m.drawparallels([-90,-60,-30,0,30,60,90],labels=[1,0,0,0],fontsize=10)

m.drawmeridians([-180,-120,-60,0,60,120,180],labels=[0,0,0,1],fontsize=10)

m.drawmapboundary(fill_color='aqua')

img = m.imshow(data, interpolation='none', cmap='jet' )

m.colorbar(img)

plt.title('Plot gridded data on a map')

plt.savefig('gridded_data_global_map.png', pad_inches=0.5, bbox_inches='tight')

plt.show()

Find corresponding MODIS granule files for a month of data using python

					
#!/usr/bin/env python

import os
import glob
import numpy as np
import calendar

#----------------------------------------------------------------------------------------#
# user config.

year = 2008
month = 07

path_to_myd021km = '/Volumes/Data/MYD021KM/'
path_to_myd03 = '/Volumes/Data/Myd03/'
path_to_myd35 = '/Volumes/DATA_03/modis/Myd35/'
path_to_myd06_c5 = '/Volumes/DATA_03/modis/Myd06_c5/' 
path_to_myd06_c6 = '/Volumes/DATA_03/modis/Myd06_c6/'

result_file_name = 'modis_granule_corresponding_files_'+ str(year) + \
                   '_' + "%02d" % month +'.txt'

#----------------------------------------------------------------------------------------#
# create file

f = open(result_file_name, 'w')

day_of_year_0 = 0
for i in range(month-1):
	day_of_year_0 = day_of_year_0 + calendar.monthrange(year,i+1)[1]

for day in np.arange(calendar.monthrange(year,month)[1]):

	day_of_year = day_of_year_0 + day + 1

	#------------------------------------------------------------------------------------# 
	# create paths

	d = str(year) + '/' + str(year) + '_' + "%02d" % month + '_' + "%02d" % (day+1) + '/'

	path_to_myd021km_day = path_to_myd021km + d
	path_to_myd03_day = path_to_myd03 + d
	path_to_myd35_day = path_to_myd35 + d
	path_to_myd06_c5_day = path_to_myd06_c5 + d
	path_to_myd06_c6_day = path_to_myd06_c6 + d

	#------------------------------------------------------------------------------------# 
	# get all file names

	myd021km_day_list = glob.glob(path_to_myd021km_day + '/*hdf' )
	myd03_day_list = glob.glob(path_to_myd03_day + '/*hdf' )
	myd35_day_list = glob.glob(path_to_myd35_day + '/*hdf' )
	myd06_c5_day_list = glob.glob(path_to_myd06_c5_day + '/*hdf' )
	myd06_c6_day_list = glob.glob(path_to_myd06_c6_day + '/*hdf' )

	#------------------------------------------------------------------------------------# 
	# create file

	for modis_granule_name in myd03_day_list:

		offset = len(path_to_myd03_day)
		myd03_name = modis_granule_name[offset:]
		hour = modis_granule_name[offset+15:offset+17]
		minute = modis_granule_name[offset+17:offset+19]
		#print hour, minute

		result = str(year) + ' ' + str(month) + ' ' + str(day+1) + \
		         ' ' + hour + ' ' + minute

		result = result + ' ' + myd03_name

		#--------------------------------------------------------------------------------#
		# MYD021KM
		
		file_name_prefix = 'MYD021KM.A'+ str(year) + "%03d" % day_of_year + \
		                   '.' + hour + minute
		found_file = 0
		for myd021km_granule_name in myd021km_day_list:
			offset = len(path_to_myd021km_day)
			if file_name_prefix == myd021km_granule_name[offset:offset+22]:
				found_file = 1
				myd021km_name = myd021km_granule_name[offset:]
		if found_file == 0: myd021km_name = 'NoFile'

		result = result + ' ' + myd021km_name

		#--------------------------------------------------------------------------------#
		# MYD35

		file_name_prefix = 'MYD35_L2.A'+ str(year) + "%03d" % day_of_year + \
		                   '.' + hour + minute
		found_file = 0
		for myd35_granule_name in myd35_day_list:
			offset = len(path_to_myd35_day)
			if file_name_prefix == myd35_granule_name[offset:offset+22]:
				found_file = 1
				myd35_name = myd35_granule_name[offset:]
		if found_file == 0: myd35_name = 'NoFile'

		result = result + ' ' + myd35_name

		#--------------------------------------------------------------------------------#
		# MYD06_c5

		file_name_prefix = 'MYD06_L2.A'+ str(year) + "%03d" % day_of_year + \
		                   '.' + hour + minute
		found_file = 0		
		for myd06_c5_granule_name in myd06_c5_day_list:
			offset = len(path_to_myd06_c5_day)
			if file_name_prefix == myd06_c5_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c5_name = myd06_c5_granule_name[offset:]
		if found_file == 0: myd06_c5_name = 'NoFile'		

		result = result + ' ' + myd06_c5_name

		#--------------------------------------------------------------------------------#
		# MYD06_c6

		file_name_prefix = 'MYD06_L2.A'+ str(year) + "%03d" % day_of_year + \
		                   '.' + hour + minute
		found_file = 0		
		for myd06_c6_granule_name in myd06_c6_day_list:
			offset = len(path_to_myd06_c6_day)
			if file_name_prefix == myd06_c6_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c6_name = myd06_c6_granule_name[offset:]
		if found_file == 0: myd06_c6_name = 'NoFile'	

		result = result + ' ' + myd06_c6_name

		#--------------------------------------------------------------------------------#
		# save result

		f.write(result + '\n')

f.close()

Download a month of MODIS L2 hdf files from lads using python

					
#!/usr/bin/env python

from ftplib import FTP

import numpy as np
import os
import calendar

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

year = 2012
month = 2

#----------------------------------------------------------------------------------------#
# Create directory 

os.system( 'mkdir ' + str(year) )

for day in np.arange(calendar.monthrange(year,month)[1]):
	d = './' + str(year) + '/' + str(year) + '_' + "%02d" % month + \
		 '_' + "%02d" % (day+1) + '/'
	os.system('mkdir ' + d)

#----------------------------------------------------------------------------------------#
# download files from NASA ladsftp 

download_myd06_c6 = False
download_myd03 = False
download_myd021km = False
download_myd06_c5 = True

print 'Opening FTP connection'

ftp = FTP('ladsftp.nascom.nasa.gov')
ftp.login('anonymous','')

day_of_year = 0 
for i in range(month-1):
    day_of_year = day_of_year + calendar.monthrange(year,i+1)[1]

#----------------------------------------------------------------------------------------#
# download MYD06 C6 File

if download_myd06_c6:

	for day in np.arange(calendar.monthrange(year,month)[1]):

		directory = 'allData/6/MYD06_L2/'+str(year)+'/'+"%03d"%(day_of_year+day+1)
		ftp.cwd(directory)
		print directory

		yyy = []
		ftp.retrlines('NLST', yyy.append)
		print 'len(yyy)', len(yyy)

		if len(yyy) != 288:
			print 'Warning missing data !'
		for j in np.arange(len(yyy)):
			file_name = yyy[j]
			print file_name	
			
			ftp.retrbinary('RETR ' + file_name, open(file_name, 'wb').write)
			command = 'mv ' + file_name + ' ' + './' + str(year) + '/' + str(year) + \
			'_' + "%02d" % month + '_' + "%02d" % (day+1) + '/' + file_name 
			os.system(command)

#----------------------------------------------------------------------------------------#
# download MYD03 File 

if download_myd03:

	for day in np.arange(calendar.monthrange(year,month)[1]):

		directory = '/allData/6/MYD03/'+str(year)+'/'+"%03d"%(day_of_year+day+1)
				
		ftp.cwd(directory)

		yyy = []
		ftp.retrlines('NLST', yyy.append)

		print 'len(yyy)', len(yyy)

		if len(yyy) != 288:
			print 'Warning missing data !'
		for j in np.arange(len(yyy)):
			file_name = yyy[j]
			print file_name	
			
			ftp.retrbinary('RETR ' + file_name, open(file_name, 'wb').write)
			command = 'mv ' + file_name + ' ' + './' + str(year) + '/' + str(year) + \
			'_' + "%02d" % month + '_' + "%02d" % (day+1) + '/' + file_name 
			os.system(command)

#----------------------------------------------------------------------------------------#
# download MYD021km File 

if download_myd021km:

	for day in np.arange(calendar.monthrange(year,month)[1]):

		directory = '/allData/6/MYD021KM/'+str(year)+'/'+"%03d"%(day_of_year+day+1)
		ftp.cwd(directory)

		yyy = []
		ftp.retrlines('NLST', yyy.append)

		print 'len(yyy)', len(yyy)

		if len(yyy) != 288:
			print 'Warning missing data !'
		for j in np.arange(len(yyy)):
			file_name = yyy[j]
			print file_name	
			
			ftp.retrbinary('RETR ' + file_name, open(file_name, 'wb').write)
			command = 'mv ' + file_name + ' ' + './' + str(year) + '/' + str(year) + \
			'_' + "%02d" % month + '_' + "%02d" % (day+1) + '/' + file_name 
			os.system(command)

#----------------------------------------------------------------------------------------#
# download MYD06 C51 File

if download_myd06_c5:

	for day in np.arange(calendar.monthrange(year,month)[1]):

		directory = '/allData/51/MYD06_L2/'+str(year)+'/'+"%03d"%(day_of_year+day+1)
		ftp.cwd(directory)

		yyy = []
		ftp.retrlines('NLST', yyy.append)

		print 'len(yyy)', len(yyy)

		if len(yyy) != 288:
			print 'Warning missing data !'
		for j in np.arange(len(yyy)):
			file_name = yyy[j]
			print file_name	
			
			ftp.retrbinary('RETR ' + file_name, open(file_name, 'wb').write)
			command = 'mv ' + file_name + ' ' + './' + str(year) + '/' + str(year) + \
			'_' + "%02d" % month + '_' + "%02d" % (day+1) + '/' + file_name 
			os.system(command)

#----------------------------------------------------------------------------------------#
# Closing FTP connection

print 'Closing FTP connection'
ftp.close()

Plot MODIS granule RGB image using python (with color enhancement)

					
#!/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()

Download MODIS L2 granule files from NASA lads using python and ftp

					
#!/usr/bin/env python

from ftplib import FTP

import numpy as np
import os
import calendar
import pprint
import sys

InpArgLIst = sys.argv

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

download_myd06_c6 = True
download_myd03 = True
download_myd021km = True
download_myd06_c5 = True

year = 2005 
month = 4 
day = 1 
hour = 13 
minute = 15 

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

if len(InpArgLIst) == 6 :
	year = int( InpArgLIst[1] )
	month = int( InpArgLIst[2] )
	day = int( InpArgLIst[3] )
	hour = int( InpArgLIst[4] ) 
	minute = int( InpArgLIst[5] ) 

print 'Year: ', year
print 'month: ', month
print 'day: ', day
print 'hour: ', hour
print 'minute: ', minute

print 'Opening FTP connection'

ftp = FTP('ladsftp.nascom.nasa.gov')
ftp.login('anonymous','')

day_of_year = 0 
for i in range(month-1):
    day_of_year = day_of_year + calendar.monthrange(year,i+1)[1]

day_of_year = day_of_year + day

day_of_year = "%03d" % day_of_year

ind = str(year) + day_of_year + '.' +  "%02d" % hour + "%02d" % minute
print ind

#----------------------------------------------------------------------------------------#
# get MYD06 C6 File

if download_myd06_c6:

	directory = 'allData/6/MYD06_L2/'+str(year)+'/'+day_of_year
	ftp.cwd(directory)
	print directory

	yyy = []
	ftp.retrlines('NLST', yyy.append)
	#pprint.pprint(yyy)

	matching = [s for s in yyy if ind in s]
	print matching
	print matching[0]

	e = 1
	n = 0 
	while e == 1 and n < 10 :
		try:
			ftp.retrbinary('RETR ' + matching[0], open(matching[0], 'wb').write)
			e = 0 
		except:
			n = n + 1
			print("Error while downloading")

#----------------------------------------------------------------------------------------#
# get MYD03 File 

if download_myd03:

	directory = '/allData/6/MYD03/'+str(year)+'/'+day_of_year
	ftp.cwd(directory)

	yyy = []
	ftp.retrlines('NLST', yyy.append)

	matching = [s for s in yyy if ind in s]
	print matching
	print matching[0]

	e = 1
	n = 0 
	while e == 1 and n < 10 :
		try:
			ftp.retrbinary('RETR ' + matching[0], open(matching[0], 'wb').write)
			e = 0 
		except:
			n = n + 1
			print("Error while downloading")

#----------------------------------------------------------------------------------------#
# get MYD021km File 

if download_myd021km:

	directory = '/allData/6/MYD021KM/'+str(year)+'/'+day_of_year
	ftp.cwd(directory)

	yyy = []
	ftp.retrlines('NLST', yyy.append)

	matching = [s for s in yyy if ind in s]
	print matching
	print matching[0]

	e = 1
	n = 0 
	while e == 1 and n < 10 :
		try:
			ftp.retrbinary('RETR ' + matching[0], open(matching[0], 'wb').write)
			e = 0 
		except:
			n = n + 1
			print("Error while downloading")

#----------------------------------------------------------------------------------------#
# get MYD06 C51 File

if download_myd06_c5:

	directory = '/allData/51/MYD06_L2/'+str(year)+'/'+day_of_year
	ftp.cwd(directory)

	yyy = []
	ftp.retrlines('NLST', yyy.append)

	matching = [s for s in yyy if ind in s]
	print matching
	print matching[0]

	e = 1
	n = 0 
	while e == 1 and n < 10 :
		try:
			ftp.retrbinary('RETR ' + matching[0], open(matching[0], 'wb').write)
			e = 0 
		except:
			n = n + 1
			print("Error while downloading")

#----------------------------------------------------------------------------------------#
# Closing FTP connection

print 'Closing FTP connection'
ftp.close()

Plot Lorentz-Mie scattering calculation results using python and matplotlib

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------------------------#
# read Mie Parameter, refractive index, qe, albedo, gfact

f = open('mie_code_output.txt', 'r')

list = f.readline()
x = list.split()[0]
x = float(x.strip())
x = round(x,2) 
nr = list.split()[1]
nr = float(nr.strip())
nr = round(nr,2) 
ni = list.split()[2]
ni = float(ni.strip())
ni = round(ni,2) 

list = f.readline()
list.strip()
qe = list.split()[0]
qe = float(qe.strip())
qe = round(qe,2) 
albedo = list.split()[1]
albedo = float(albedo.strip())
albedo = round(albedo,2)
gfact = list.split()[2]
gfact = float(gfact.strip())
gfact = round(gfact,2) 

f.closed

#----------------------------------------------------------------------------------------#
# read scattering sngle, p11, p12, p33, p34

c1,c2,c3,c4,c5  = np.loadtxt("mie_code_output.txt", skiprows=2, unpack=True)

#----------------------------------------------------------------------------------------#
# print p11

fig = plt.figure()

ax = fig.add_subplot(111)

plt.plot(c1,c2,'k-')
plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{11}$")

plt.text(0.70,0.85,'x = '+str(x), fontsize=12, transform = ax.transAxes)
plt.text(0.70,0.80,'n=('+str(nr)+','+str(ni)+')', fontsize=12, transform = ax.transAxes)     

plt.text(0.70,0.70,'qe = '+str(qe), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.65,'w0 = '+str(albedo), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.60,'g = '+str(gfact), fontsize=12, transform = ax.transAxes)     

plt.savefig('p11.png')

#----------------------------------------------------------------------------------------#
# print p12

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(c1,c3,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{12}$")    

plt.savefig('p12.png')

#----------------------------------------------------------------------------------------#
# print p33

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(c1,c4,'k-')
plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{33}$")    

plt.savefig('p33.png')

#----------------------------------------------------------------------------------------#
# print p34

fig = plt.figure()

ax = fig.add_subplot(111)

plt.plot(c1,c5,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{34}$")    

plt.savefig('p34.png')

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

Lorentz-Mie scattering calculation using fortran

					
program mie_code

! reference: Hong Du: Mie-Scattering Calculation, Applied Optics, 
!                     Vol. 43, Issue 9, pp. 1951-1956 (2004)

implicit none

integer :: i,j,k

integer, parameter :: n = 100  ! number of coefficients 

double precision :: x = 5.0  ! Mie Parameter (2*pi*r/wl)

double complex :: ri = (1.33d0,-0.0d0)  ! refractive index 

double complex, dimension(n) :: ma , mb  

double precision, dimension(:,:), allocatable :: pi , tau  

double precision, dimension(0:n) :: phi , chi 
double complex, dimension(0:n) :: zeta  

double complex, dimension(0:n) :: rn  

double complex, dimension(:), allocatable :: sa1 , sa2 

double precision, dimension(:), allocatable :: p11,p12,p33,p34  

double precision :: qe

double precision :: albedo , gfact

integer :: nb_angles

double precision :: teta
double precision :: teta_resol = 1.0 ! scattering angle resolution

double precision :: s , t

double precision :: temp1 , temp2 , temp3
double complex :: ctemp1 , ctemp2 ,ctemp3

double precision :: coef1 , coef2 , coef3 , coef4

double precision :: norm

double precision :: cdr = 0.01745329252

double complex :: ci = (0.0d0,1.0d0) 

!----------------------------------------------------------------------------------------!

if (n < int(x)) write(6,*) 'Warning: number of coeffs n is not enough !'

nb_angles = int(180/teta_resol) + 1

allocate( pi(nb_angles,n) )
allocate( tau(nb_angles,n) )

allocate( sa1(nb_angles) )
allocate( sa2(nb_angles) )

allocate( p11(nb_angles) )
allocate( p12(nb_angles) )
allocate( p33(nb_angles) )
allocate( p34(nb_angles) )

do i = 1 , n
	ma(i)  = (0.0d0,0.0d0)
	mb(i)  = (0.0d0,0.0d0)
end do

do i = 1 , nb_angles
do j = 1 , n
	pi(i,j)  = 0.0d0
	tau(i,j) = 0.0d0
end do
end do

do i = 1 ,nb_angles
	sa1(i) = 0.0d0
	sa2(i) = 0.0d0
end do

!----------------------------------------------------------------------------------------!
! Step 1: Calculations of the Mie angular functions.

teta = 0.0d0

do i = 1 , nb_angles

!---- pi function

temp1 = dcos( teta * cdr )

temp2 = 0.0d0
temp3 = 1.0d0
pi(i,1) = temp3
do j = 2 , n
	s = temp1 * temp3
	t = s - temp2
	pi(i,j) = s + t + t / ( j - 1 )
	temp2 = temp3
	temp3 = pi(i,j)
end do

!----- tau function

tau(i,1) = temp1
do j = 2 , n
	t = temp1 * pi(i,j) - pi(i,j-1)
	tau(i,j) = j * t - pi(i,j-1)
end do

teta = teta + teta_resol
end do

!----------------------------------------------------------------------------------------!
! Step 2: Calculations of the Riccati-Bessel functions by upward recurrence.

!----- phi function

temp1 = dcos( x )
phi(0) = dsin( x )
temp2 = phi(0) 
do i = 1 , n
	phi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
	temp1 = temp2
	temp2 = phi(i)
end do

!----- chi function

chi(0) = dcos(x)
temp1  = chi(0) 
chi(1) = dcos(x) / x + dsin(x)
temp2  = chi(1) 
do i = 2 , n
	chi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
	temp1 = temp2
	temp2 = chi(i)
end do

!----- zeta fcuntion

do i = 0 , n
	zeta(i) = phi(i) + ci * chi(i) 
end do

!----------------------------------------------------------------------------------------!
! Step 3: Calculations of the complex functions rn(mx) by upward recurrence.

temp1 = dble( ri ) * x
temp2 = aimag( ri ) * x

coef1 = dtan( temp1 ) - dexp( - 2.0 * temp2 ) * dtan( temp1 )
coef2 = 1.0d0 + dexp( - 2.0 * temp2 )
coef3 = - 1.0d0 + dexp( - 2.0 * temp2 )
coef4 = dtan( temp1 ) + dexp( - 2.0 * temp2 ) * dtan( temp1 )

temp1 = coef1 * coef3 + coef2 * coef4
temp2 = coef2 * coef3 - coef1 * coef4
temp3 = coef3 * coef3 + coef4 * coef4

rn(0) = ( temp1 / temp3 ) + ci * ( temp2 / temp3 )

do i = 1 , n
	ctemp1 = ( 2 * i - 1 ) / ( ri * x ) - rn(i-1)
	rn(i) = 1.0d0 / ctemp1
end do

!----------------------------------------------------------------------------------------!
! Step 4: Calculations of the Mie coefficients an et bn.

do i = 1 , n
	ctemp1 = rn(i) / ri + i * ( 1.0d0 - 1.0d0 / ( ri * ri ) ) / x
	ma(i) = ( ctemp1 * phi(i) - phi(i-1) ) / ( ctemp1 * zeta(i) - zeta(i-1) )
	mb(i) = ( rn(i) * ri  * phi(i) - phi(i-1) ) / ( rn(i) * ri  * zeta(i) - zeta(i-1) )
end do

!----------------------------------------------------------------------------------------!
! Step 5: scattering amplitudes

do i = 1 , nb_angles
	ctemp2 = (0.0d0,0.0d0)
	ctemp3 = (0.0d0,0.0d0)
	do j = 1 , n
		temp1 = DBLE( 2 * j + 1 ) / DBLE( j * ( j + 1 ) )
		ctemp2 = ctemp2 + temp1 * ( ma(j) * pi(i,j) + mb(j) * tau(i,j) )
		ctemp3 = ctemp3 + temp1 * ( mb(j) * pi(i,j) + ma(j) * tau(i,j) )
	end do
	sa1(i) = ctemp2
	sa2(i) = ctemp3
end do

!----------------------------------------------------------------------------------------!
! Step 6: phase function.

!----- norm (need to be improved <- integration too simple)

norm = 0.0
do i = 1 , n
	temp1 = abs( ma(i) ) * abs( ma(i) ) + abs( mb(i) ) * abs( mb(i) )
	norm = norm + ( 2.0d0 * i + 1.0d0 ) * temp1
end do

do i = 1 , nb_angles
	p11(i) = ( sa1(i) * conjg( sa1(i) ) + sa2(i) * conjg( sa2(i) ) ) / norm
	p12(i) = ( sa2(i) * conjg( sa2(i) ) - sa1(i) * conjg( sa1(i) ) ) / norm 
	p33(i) = ( sa2(i) * conjg( sa1(i) ) + sa1(i) * conjg( sa2(i) ) ) / norm 
	p34(i) = ci * ( sa1(i) * conjg( sa2(i) ) - sa2(i) * conjg( sa1(i) ) ) / norm 
end do

!----------------------------------------------------------------------------------------!
! extinction efficiency.

temp1 = 0.0d0
do i = 1 , n
	temp1 = temp1 + ( 2 * i + 1 ) * dble( ma(i) + mb(i) ) 
end do

qe = ( 2.0d0 / ( x * x ) ) * temp1

!----------------------------------------------------------------------------------------!
! asymmetry factor.

temp1 = 0.0d0
do i = 1 , n - 1
	ctemp1 = ma(i) * conjg( ma(i+1) ) + mb(i) * conjg( mb(i+1) ) 
	ctemp2 = ma(i) * conjg( mb(i) )
	temp2 = dble( i * ( i + 2.0d0 ) ) / dble( i + 1.0d0 )
	temp3 = dble( 2.0d0 * i + 1.0d0 ) / dble( i * ( i + 1.0d0 ) )
	temp1 = temp1 + ( temp2 * dble( ctemp1 ) + temp3 * dble( ctemp2 ) )
end do

gfact = ( 2.0d0 / norm ) * temp1

!----------------------------------------------------------------------------------------!
! single scattering albedo.

albedo = ( 2.0d0 / ( x * x ) ) * norm / qe

!----------------------------------------------------------------------------------------!
! save data in file

teta = 0.0d0
open(1,file='mie_code_output.txt')
	write(1,*) x, real(real(ri)), real(AIMAG(ri)) !, wl
	write(1,*) qe,albedo,gfact
	do i = 1 , nb_angles
		write(1,*) teta,p11(i),p12(i),p33(i),p34(i)
		teta = teta + teta_resol
	end do
close(1)

end program mie_code

How to plot a MODIS granule RGB image using python 3 ?

					
#!/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.items():
	#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.items():
	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))

z[:,:,0] = ( sds_data_250_Aggr1km_RefSB[0,:,:] - reflectance_offsets_250_Aggr1km_RefSB[0] ) * reflectance_scales_250_Aggr1km_RefSB[0] 
z[:,:,1] = ( sds_data_500_Aggr1km_RefSB[1,:,:] - reflectance_offsets_500_Aggr1km_RefSB[1] ) * reflectance_scales_500_Aggr1km_RefSB[1] 
z[:,:,2] = ( sds_data_500_Aggr1km_RefSB[0,:,:] - 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))

rgb = z / norme

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

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

fig = plt.figure()

ax = fig.add_subplot(111)

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

l = [int(i) for i in np.linspace(0,cross_trak,6)]
plt.xticks(l, [i for i in reversed(l)], rotation=0, fontsize=7 )

l = [int(i) for i in np.linspace(0,along_track,9)]
plt.yticks(l, l, rotation=0, fontsize=7 )

plt.xticks(fontsize=7)
plt.yticks(fontsize=7)

plt.title('How to plot a MODIS RGB image \n using python 3 ?', fontsize=8)

plt.savefig("modis_granule_rgb.png", bbox_inches='tight', dpi=100)

#plt.show()

plt.close()