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

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2016-10-18

Edited / Modifié: 2016-10-18

License / Licence: MIT License

Abstract / Résumé