#!/usr/bin/env python

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

import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------------------------#
# 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_1km = SD('CAL_LID_L2_01kmCLay-ValStage1-V3-01.2008-07-03T08-17-58ZD.hdf', SDC.READ)

caliop_latitude_1km = caliop_1km.select('Latitude')
caliop_longitude_1km = caliop_1km.select('Longitude')
caliop_Profile_Time_UTC_1km = caliop_1km.select('Profile_UTC_Time')
caliop_Number_Layers_Found_1km = caliop_1km.select('Number_Layers_Found')
caliop_Layer_Top_Altitude_1km = caliop_1km.select('Layer_Top_Altitude')
caliop_Layer_Base_Altitude_1km = caliop_1km.select('Layer_Base_Altitude')
caliop_Feature_Classification_Flags_1km = \
                                 caliop_1km.select('Feature_Classification_Flags')

caliop_latitude_1km_data = caliop_latitude_1km.get()
caliop_longitude_1km_data = caliop_longitude_1km.get()
caliop_Profile_Time_UTC_1km_data = caliop_Profile_Time_UTC_1km.get()
caliop_Number_Layers_Found_1km_data = caliop_Number_Layers_Found_1km.get()
caliop_Layer_Top_Altitude_1km_data = caliop_Layer_Top_Altitude_1km.get()
caliop_Layer_Base_Altitude_1km_data = caliop_Layer_Base_Altitude_1km.get()
caliop_Feature_Classification_Flags_1km_data = \
                                      caliop_Feature_Classification_Flags_1km.get()

print "Read caliop 1km done"

caliop_5km = SD('CAL_LID_L2_05kmCLay-Prov-V3-01.2008-07-03T08-17-58ZD.hdf', SDC.READ)

caliop_latitude_5km = caliop_5km.select('Latitude')
caliop_longitude_5km = caliop_5km.select('Longitude')
caliop_Profile_Time_UTC_5km = caliop_5km.select('Profile_UTC_Time')
caliop_Number_Layers_Found_5km = caliop_5km.select('Number_Layers_Found')
caliop_Layer_Top_Altitude_5km = caliop_5km.select('Layer_Top_Altitude')
caliop_Layer_Base_Altitude_5km = caliop_5km.select('Layer_Base_Altitude')
caliop_Feature_Classification_Flags_5km = \
                                 caliop_5km.select('Feature_Classification_Flags')

caliop_latitude_5km_data = caliop_latitude_5km.get()
caliop_longitude_5km_data = caliop_longitude_5km.get()
caliop_Profile_Time_UTC_5km_data = caliop_Profile_Time_UTC_5km.get()
caliop_Number_Layers_Found_5km_data = caliop_Number_Layers_Found_5km.get()
caliop_Layer_Top_Altitude_5km_data = caliop_Layer_Top_Altitude_5km.get()
caliop_Layer_Base_Altitude_5km_data = caliop_Layer_Base_Altitude_5km.get()
caliop_Feature_Classification_Flags_5km_data = \
                                      caliop_Feature_Classification_Flags_5km.get()

print "Read caliop 5km done"

#print caliop_latitude_1km_data.shape
#print caliop_latitude_5km_data.shape

caliop_01km_shape = caliop_latitude_1km_data.shape
caliop_05km_shape = caliop_latitude_5km_data.shape

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

start_hour = 8
start_minute = 31
start_second = 24

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

end_hour = 8
end_minute = 44
end_second = 53

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

#----------------------------------------------------------------------------------------#
# Main Loop

fig = plt.figure(figsize=(12,6))
#fig.patch.set_facecolor('0.65')
fig = fig.add_subplot(111)

lat_min = 999
lat_max = -999

cloudy_layer_base_altitude_1km = np.zeros(10)
cloudy_layer_top_altitude_1km = np.zeros(10)
cloudy_layer_phase_1km = np.zeros(10)

cloudy_layer_base_altitude_5km = np.zeros(10)
cloudy_layer_top_altitude_5km = np.zeros(10)
cloudy_layer_phase_5km = np.zeros(10)

merged_layer_base_altitude = np.zeros(10)
merged_layer_top_altitude = np.zeros(10)
merged_layer_phase_1km = np.zeros(10)

tmp_merged_layer_base_altitude = np.zeros(10)
tmp_merged_layer_btop_altitude = np.zeros(10)
tmp_merged_layer_phase_1km = np.zeros(10)

for i in np.arange( caliop_01km_shape[0] ):
	tmp =  str( caliop_Profile_Time_UTC_1km_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		
		if caliop_latitude_1km_data[i,0] < lat_min : 
			lat_min = caliop_latitude_1km_data[i,0]
		if caliop_latitude_1km_data[i,0] > lat_max : 
			lat_max = caliop_latitude_1km_data[i,0]	
		bar_xi = caliop_latitude_1km_data[i-1,0] + \
		abs( caliop_latitude_1km_data[i,0] - caliop_latitude_1km_data[i-1,0] ) / 2.0
		bar_width = abs( caliop_latitude_1km_data[i,0] - 
		            caliop_latitude_1km_data[i-1,0] )
		#---------- 1km cloud vertical profile ----------#
		number_cloudy_layers_found_1km = 0
		if caliop_Number_Layers_Found_1km_data[i,0] >= 1 :
			for j in np.arange( caliop_Number_Layers_Found_1km_data[i,0] ):
				feature_type = bits_stripping(0,3,
				               caliop_Feature_Classification_Flags_1km_data[i,j])
				if feature_type == 2: 
					cloudy_layer_base_altitude_1km[number_cloudy_layers_found_1km] = \
					caliop_Layer_Base_Altitude_1km_data[i,j] 
					cloudy_layer_top_altitude_1km[number_cloudy_layers_found_1km] = \
					caliop_Layer_Top_Altitude_1km_data[i,j] 		
					cloudy_layer_phase_1km[number_cloudy_layers_found_1km] = \
					bits_stripping(5,2,caliop_Feature_Classification_Flags_1km_data[i,j])
					number_cloudy_layers_found_1km = number_cloudy_layers_found_1km + 1		
		#---------- 5km cloud vertical profile ----------#
		k = int(i/5.0)
		number_cloudy_layers_found_5km = 0
		if caliop_Number_Layers_Found_5km_data[k,0] >= 1 :
			for j in np.arange( caliop_Number_Layers_Found_5km_data[k,0] ):
				feature_type = bits_stripping(0,3,
				caliop_Feature_Classification_Flags_5km_data[k,j])
				if feature_type == 2: 				
					cloudy_layer_base_altitude_5km[number_cloudy_layers_found_5km] = \
					caliop_Layer_Base_Altitude_5km_data[k,j] 
					cloudy_layer_top_altitude_5km[number_cloudy_layers_found_5km] = \
					caliop_Layer_Top_Altitude_5km_data[k,j] 
					cloudy_layer_phase_5km[number_cloudy_layers_found_5km] = \
					bits_stripping(5,2,caliop_Feature_Classification_Flags_5km_data[k,j])		
					number_cloudy_layers_found_5km = number_cloudy_layers_found_5km + 1			
		#---------- Merged cloud vertical profile ----------#
		merged_number_cloudy_layers_found = 0
		for j in np.arange( number_cloudy_layers_found_1km ):
			merged_layer_base_altitude[merged_number_cloudy_layers_found] = \
			cloudy_layer_base_altitude_1km[j]
			merged_layer_top_altitude[merged_number_cloudy_layers_found] = \
			cloudy_layer_top_altitude_1km[j]
			merged_layer_phase_1km[merged_number_cloudy_layers_found] = \
			cloudy_layer_phase_1km[j]			
			merged_number_cloudy_layers_found = merged_number_cloudy_layers_found + 1
		for j in np.arange( number_cloudy_layers_found_5km ):
			merged_layer_base_altitude[merged_number_cloudy_layers_found] = \
			cloudy_layer_base_altitude_5km[j]
			merged_layer_top_altitude[merged_number_cloudy_layers_found] = \
			cloudy_layer_top_altitude_5km[j]
			merged_layer_phase_1km[merged_number_cloudy_layers_found] = \
			cloudy_layer_phase_5km[j]			
			merged_number_cloudy_layers_found = merged_number_cloudy_layers_found + 1
		#---------- Plot ----------#					
		if merged_number_cloudy_layers_found >= 1:
			for j in np.arange( merged_number_cloudy_layers_found ):
				bar_base = merged_layer_base_altitude[j]
				bar_height = merged_layer_top_altitude[j] - bar_base
				caliop_layer_phase = merged_layer_phase_1km[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

title = 'CALIOP Cloud Thermodynamic Phase \n Merged 01 & 05 CLay'
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("caliop_01_and_05CLay_cloud_phase.png", dpi=80 )

plt.show()

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2016-10-18

Edited / Modifié: 2016-10-18

License / Licence: MIT License

Abstract / Résumé