Download a month of CALIOP L2 hdf files from ICARE using python

					
#!/usr/bin/env python

import numpy as np
import os
import calendar

from ftplib import FTP

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

year = 2008
month = 07

icare_username = 'john.doe'
icare_password = '********'

#----------------------------------------------------------------------------------------#
# 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)

#----------------------------------------------------------------------------------------#
# Opening ICARE FTP connection

print 'Opening FTP connection'

ftp = FTP('ftp.icare.univ-lille1.fr')
ftp.login(icare_username,icare_password)

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

#----------------------------------------------------------------------------------------#
# download files 

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

	directory = '/DATA/LIENS/CALIOP/01kmCLay.v4.10/'+str(year)+'/' + \
				str(year)+'_'+ "%02d" % month + '_' + "%02d" % (day+1) + '/'

	ftp.cwd(directory)

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

	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)

print 'Closing FTP connection'
print ftp.close()

Plot MODIS L3 chlorophyll concentration with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC
from mpl_toolkits.basemap import Basemap, cm
from matplotlib.colors import LogNorm

from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import math

#---------- read data ----------#

File_Name = 'A20080012008031.L3m_MO_CHL_chlor_a_9km'
File = SD(File_Name, SDC.READ)
l3chlor = File.select('l3m_data')
l3chlor_data = l3chlor.get()

print l3chlor_data.shape
print l3chlor_data.max()
print l3chlor_data.min()

FillValue = -32767.0

l3chlor_data[ (l3chlor_data < 0.01) & (l3chlor_data != -32767.0) ] = 0.011
l3chlor_data[ l3chlor_data > 10 ] = 10.0
l3chlor_data[ l3chlor_data == FillValue ] = 0.00001

#---------- Plot Data Global Map ----------#

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(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
m.drawmapboundary(fill_color='aqua')

colors = [(0.33,0.33,0.33)] + [(plt.cm.jet(i)) for i in xrange(1,256)]
new_map = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)

img = m.imshow(np.flipud(l3chlor_data), norm=LogNorm(vmin=0.01, vmax=10.0), \
			   cmap=new_map, interpolation='bilinear')
cb = m.colorbar(img,"bottom", size="5%", pad='4%')

plt.title('MODIS/Aqua Chlorophyll Concentration \n January 2008')
plt.savefig('ExampleChlorophyllConcentration.png',dpi=200)
plt.show()

Read and extract hours and minutes from a POLDER granule using python

					
#!/usr/bin/env python

from numpy import *
import numpy as np
import sys
import subprocess
import glob
import os

print 'Number of arguments:', len(sys.argv), 'arguments.'
print 'Argument List:', str(sys.argv)

print sys.argv[0], sys.argv[1]
polder_file_name = sys.argv[1]

command = 'parasolascii -g ' + polder_file_name + ' phase > ' + polder_file_name + '.txt'
os.system(command)

hour = []
i=0
command = ['parasolascii', polder_file_name, 'utc_hour']
process = subprocess.Popen(command, stdout=subprocess.PIPE)
while True:
	text = process.stdout.readline()
	if text != '':
		hour.append( int(text) )
		#print hour[i],i
		i=i+1
	else:
		break

minute = []
i=0
command = ['parasolascii', polder_file_name, 'utc_minute']
process = subprocess.Popen(command, stdout=subprocess.PIPE)
while True:
	text = process.stdout.readline()
	if text != '':
		minute.append( int(text) )
		#print minute[i],i
		i=i+1
	else:
		break
print len(hour),len(minute)
time = []
for i in range(len(hour)):
	time.append( hour[i] * 60.0 + minute[i]	)
start_time = min(time)
end_time = max(time)
if (end_time-start_time)>60.0:		
	print 'hello', alen(hour)
	for j in np.arange(alen(hour)):
		if(hour[j]<12.0): 
			time[j] = 24.0*60.0+hour[i]*60.0+minute[j]	
			start_time = min(time)
			end_time = max(time) - 24*60
start_hour = int( start_time / 60 )
start_minute = int( start_time - start_hour * 60.0 )
end_hour = int( end_time / 60 )
end_minute = int( end_time - end_hour * 60.0 )	
print start_hour, start_minute, end_hour, end_minute
if start_hour < 10:
	str_start_hour = '0' + str( start_hour )	
else:
	str_start_hour = str( start_hour )
if start_minute < 10:
	str_start_minute = '0' + str( start_minute )
else:
	str_start_minute = str( start_minute )		
if end_hour < 10:
	str_end_hour = '0' + str( end_hour )	
else:
	str_end_hour = str( end_hour )		
if end_minute < 10:
	str_end_minute = '0' + str( end_minute )
else:
	str_end_minute = str( end_minute )		
command = 'mv ' + polder_file_name + '.txt' + ' ' + polder_file_name + '_' + str_start_hour + str_start_minute + '_' + str_end_hour + str_end_minute +  '.txt'
os.system(command)

Diffraction de Fraunhofer pour une ouverture rectangulaire en fortran 90

					
program rectangular_aperture

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: wl ! WaveLength (micrometer) !
real(pr) :: wk ! WaveNumber (micrometer) !
real(pr) :: lx , ly ! aperture sides !

real(pr) :: D ! Distance from the screen to the aperture (micrometer) !
real(pr) :: x_min , x_max ! screen x_min , x_max !
real(pr) :: y_min , y_max ! screen y_min , y_max !
real(pr) :: resol_x , resol_y ! x and y resolution !

integer :: nb_x_point , nb_y_point

real(pr) :: tmp_01 , tmp_02 , tmp_03 , tmp_04
real(pr) :: x , y 
real(pr) :: r
real(pr) :: Intensity

integer :: output ! if 1 -> 3D !
    			  ! if 2 -> 2D !

!----------------------------------------------------------------------------------------!
! User Inputs

wl = 0.55d0 
lx = 10.0d0
ly = 10.0d0
D = 1000000.0d0
x_min = - 300000.0d0
x_max =	  300000.0d0
y_min = - 300000.0d0
y_max =   300000.0d0
resol_x = 1000.0d0
resol_y = 1000.0d0
output = 1

!----------------------------------------------------------------------------------------!
! Check Resolution 

wk = ( 2.0d0 * pi ) / wl

if( resol_x > ( x_max - x_min ) / 2.0d0 )then
write(6,*) "the resolution along x is too large"
end if

if( resol_y > ( y_max - y_min ) / 2.0d0 )then
write(6,*) "the resolution along y is too large"
end if

!----------------------------------------------------------------------------------------!
! Main code

open(1,file='output_data')

!----------------------------------------------------------------------------------------!
! 3d

if( output == 1 )then
nb_x_point = int( ( x_max - x_min ) / resol_x ) + 1
nb_y_point = int( ( y_max - y_min ) / resol_y ) + 1
x = x_min
do i = 1 , nb_x_point
y = y_min
do j = 1 , nb_y_point
r = dsqrt( x * x + y * y + D * D )
tmp_01 = wk * ( x / r ) * ( lx / 2.0d0 )
tmp_02 = fsc( tmp_01 ) 
tmp_03 = wk * ( y / r ) * ( ly / 2.0d0 )
tmp_04 = fsc( tmp_03 ) 
Intensity = ( tmp_02 * tmp_02 ) * ( tmp_04 * tmp_04 ) 
if( Intensity <= 0.01d0 )then ! To show the existence of the weak secondary maxima !
write(1,*) x , y , Intensity  
else 
write(1,*) x , y , 0.01d0
end if
y = y + resol_y
end do
x = x + resol_x
end do
end if

!----------------------------------------------------------------------------------------!
! 2d

if( output == 2 )then
nb_x_point = int( ( x_max - x_min ) / resol_x ) + 1
nb_y_point = int( ( y_max - y_min ) / resol_y ) + 1
x = x_min
y = 0.0d0
do i = 1 , nb_x_point 
r = dsqrt( x * x + y * y + D * D )
tmp_01 = wk * ( x / r ) * ( lx / 2.0d0 )
tmp_02 = fsc( tmp_01 ) 
tmp_03 = wk * ( y / r ) * ( ly / 2.0d0 )
tmp_04 = fsc( tmp_03 ) 
Intensity = ( tmp_02 * tmp_02 ) * ( tmp_04 * tmp_04 ) 
write(1,*) x , y , Intensity 
x = x + resol_x
end do
end if

close(1)

contains

real function fsc(x) 

real(pr) :: x

if( abs(x) > 0.001d0 ) then
	fsc = sin( x ) / x
else
	fsc = 1.0d0 - x * x / 6.0d0
end if

return

end function fsc

end program rectangular_aperture

Plot POLDER global gridded L2 cloud pressure with python

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from datetime import datetime
from pylab import *
from pyhdf.SD import SD, SDC # To read a HDF file
import numpy as np
import matplotlib.pyplot as plt
import os
import glob
import subprocess

resol = 0.5
Data_matrix = np.zeros( (180/resol,360/resol) )
mean_obs_matrix = np.zeros( (180/resol,360/resol) )
print Data_matrix.shape
Data_matrix[:,:] = - 1.0

year = '2008'
month = '01'
polder_path = '/Volumes/Data/Polder/RB2/'

PolderPressureMethod = 'po2' # 'po2' (oxygen) or 'pray' (Rayleigh)

if PolderPressureMethod == 'pray':
	PolerLabel = 'CTP (Rayleigh)'
	PolderFileLabel = 'Rayleigh'
else:
	PolerLabel = 'CTP (Oxygen)'
	PolderFileLabel = 'Oxygen'

for i in range(31):
	if i+1 < 10:
		day = '0' + str(i+1)
	else:
		day = str(i+1)
	polder_full_path = polder_path + year + '/' + year + '_' + month + '_' + day + '/' 		
	polder_day_list = glob.glob(polder_full_path + '/*KD' )
	for polder_granule_name in polder_day_list:
		lat = []
		long = []
		pressure = []
		command = ['parasolascii', '-g', polder_granule_name, PolderPressureMethod]
		process = subprocess.Popen(command, stdout=subprocess.PIPE)
		print polder_granule_name
		while True:
			text = process.stdout.readline()
			if text != '':
				polder_lat, polder_long, polder_pressure = text.split()
				lat.append( float(polder_lat) )
				long.append( float(polder_long) )
				pressure.append( float(polder_pressure) )            
			else:
				break
		#print len(lat)
		for d in range(len(lat)):
			index_i = (lat[d] + 90)/resol 
			index_j = (long[d] + 180)/resol 
			if pressure[d] > 0.0 and pressure[d] < 1100:
				mean_obs_matrix[index_i,index_j] = mean_obs_matrix[index_i,index_j] + 1.0
				Data_matrix[index_i,index_j] = Data_matrix[index_i,index_j] + pressure[d]
	print "Day",i ,"done"

MatrixShape = Data_matrix.shape
for i in range(MatrixShape[0]):
	for j in range(MatrixShape[1]):
		if mean_obs_matrix[i,j] != 0:
			Data_matrix[i,j] = Data_matrix[i,j] / mean_obs_matrix[i,j]

np.savetxt('PolderMeanCTP_'+PolderFileLabel+'.txt', Data_matrix )

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(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
m.drawmapboundary(fill_color='aqua')
plt.title("POLDER monthly mean" + PolerLabel + "\n January 2008")
#colors = [(1,1,1)] + [(cm.gist_rainbow(i)) for i in reversed(xrange(1,256))]
colors = [(1,1,1)] + [(cm.jet(i)) for i in xrange(1,256)]
new_map = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)
#img = m.imshow(A, interpolation='nearest' )
img = m.imshow(Data_matrix, cmap=new_map, \
	interpolation='bilinear', vmin=-1, vmax=1100.0 )
bounds=[0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100]
cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, \
	ticks=[0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100])
plt.savefig("PlotPolderMonthlyProductsFromL2_CTP_"+PolderFileLabel+".png", dpi=100 )
plt.show()

Plot POLDER global gridded L2 cloud cover with python

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from datetime import datetime
from pylab import *
from pyhdf.SD import SD, SDC # To read a HDF file
import numpy as np
import matplotlib.pyplot as plt
import os
import glob
import subprocess

resol = 0.5
coverage_matrix = np.zeros( (180/resol,360/resol) )
mean_obs_matrix = np.zeros( (180/resol,360/resol) )
print coverage_matrix.shape
coverage_matrix[:,:] = - 1.0

year = '2008'
month = '01'
polder_path = '/Volumes/Data/Polder/RB2/'

for i in range(31):
	if i+1 < 10:
		day = '0' + str(i+1)
	else:
		day = str(i+1)
	polder_full_path = polder_path + year + '/' + year + '_' + month + '_' + day + '/' 		
	polder_day_list = glob.glob(polder_full_path + '/*KD' )
	for polder_granule_name in polder_day_list:
		lat = []
		long = []
		cc = []
		command = ['parasolascii', '-g', polder_granule_name, 'cc']
		process = subprocess.Popen(command, stdout=subprocess.PIPE)
		#print polder_granule_name
		while True:
			text = process.stdout.readline()
			if text != '':
				polder_lat, polder_long, polder_cc = text.split()
				lat.append( float(polder_lat) )
				long.append( float(polder_long) )
				cc.append( float(polder_cc) )            
			else:
				break
		#print len(lat)
		for d in range(len(lat)):
			index_i = (lat[d] + 90)/resol 
			index_j = (long[d] + 180)/resol 
			mean_obs_matrix[index_i,index_j] = mean_obs_matrix[index_i,index_j] + 1.0
			coverage_matrix[index_i,index_j] = coverage_matrix[index_i,index_j] + 100.0 * cc[d]
	print "Day",i ,"done"

coverage_matrix[:,:] = coverage_matrix[:,:] / mean_obs_matrix[:,:]
	
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(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
m.drawmapboundary(fill_color='aqua')
plt.title("POLDER monthly cloud cover \n (defined as the number of cloudy pixels \n divided by the total number of pixels)")
#colors = [(1,1,1)] + [(cm.gist_rainbow(i)) for i in reversed(xrange(1,256))]
colors = [(0,0,0)] + [(cm.jet(i)) for i in xrange(1,256)]
new_map = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)
#img = m.imshow(A, interpolation='nearest' )
img = m.imshow(coverage_matrix, cmap=new_map, interpolation='bilinear', vmin=-1, vmax=100 )
bounds=[0,10,20,30,40,50,60,70,80,90,100]
cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, ticks=[0,10,20,30,40,50,60,70,80,90,100])
plt.savefig("PlotPolderMonthlyProductsFromL2_CC.png", dpi=100 )
plt.show()

Plot CloudSat radar reflectivity using python and matplotlib

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC # To read a HDF file (SD data)
from pyhdf.HDF import *
from pyhdf.VS import *

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

from mpl_toolkits.axes_grid1 import make_axes_locatable
import pylab

#---------- 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 (VD data) ----------#

f = HDF('2008001002019_08922_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf') 
vs = f.vstart() 
Latitude = vs.attach('Latitude')
Longitude = vs.attach('Longitude')

a = Latitude[:]

Latitude.detach() # "close" the vdata
Longitude.detach() # "close" the vdata
vs.end() # terminate the vdata interface
f.close() 

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

cpr_2b_geoprof = SD('2008001002019_08922_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf', SDC.READ)

cpr_2b_geoprof_height = cpr_2b_geoprof.select('Height')
cpr_2b_geoprof_cloud_mask = cpr_2b_geoprof.select('CPR_Cloud_mask')
cpr_2b_geoprof_radar_reflectivity = cpr_2b_geoprof.select('Radar_Reflectivity')

cpr_2b_geoprof_height_data = cpr_2b_geoprof_height.get()
cpr_2b_geoprof_cloud_mask_data = cpr_2b_geoprof_cloud_mask.get()
cpr_2b_geoprof_radar_reflectivity_data = cpr_2b_geoprof_radar_reflectivity.get()

fillvalue = 15360
missing = -8888

img = np.zeros((800,125))
img.fill(-30)

for i in np.arange(800):
for j in np.arange(125):
	k = int( 125 * (cpr_2b_geoprof_height_data[i,j] + 5000) / 35000 )
	if cpr_2b_geoprof_radar_reflectivity_data[i,j] > -3000 and \
	cpr_2b_geoprof_radar_reflectivity_data[i,j] < 2100: 		
		img[i,k] = cpr_2b_geoprof_radar_reflectivity_data[i,j] / 100

ax = plt.subplot(111)
im = ax.imshow(img.T, interpolation='bilinear', cmap=cm.jet, origin='lower', extent=[0,200,-10,60])

plt.title( 'CLOUDSAT Radar Reflectivity' )

plt.ylabel('Height (km)',fontsize=14)
pylab.yticks([-10,0,10,20,30,40,50,60], [-5,0,5,10,15,20,25,30])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.10)
plt.colorbar(im, cax=cax)

plt.savefig("cloudsat_radar_reflectivity.png")
plt.show()

Implémenter l'algorithme de Kruskal en fortran 90

					
program kruskal_example

implicit none

integer, parameter :: pr = selected_real_kind(15,3)
integer, parameter :: n = 7 ! Number of Vertice
integer, parameter :: nr = int( ( n * n - n ) / 2 )

integer :: i,i1,j,j1,k

real(pr) :: tmp_01, tmp_02, tmp_03
real(pr), dimension(n,n) :: EdgesArray 
real(pr), dimension(nr,3) :: ReducedEdgesArray 
real(pr), dimension(n-1,3) :: MST

integer :: NbEdges
integer, dimension(n,n) :: LV

logical :: FormACycle, FoundVertice_01, FoundVertice_02
logical :: FindNewEdge

data EdgesArray / 0.0, 7.0, 0.0, 5.0, 0.0, 0.0, 0.0, &
				& 7.0, 0.0, 8.0, 9.0, 7.0, 0.0, 0.0, &
				& 0.0, 8.0, 0.0, 0.0, 5.0, 0.0, 0.0, &
				& 5.0, 9.0, 0.0, 0.0, 15.0, 6.0, 0.0, &
				& 0.0, 7.0, 5.0, 15.0, 0.0, 8.0, 9.0, &
				& 0.0, 0.0, 0.0, 6.0, 8.0, 0.0, 11.0, &
				& 0.0, 0.0, 0.0, 0.0, 9.0, 11.0, 0.0 /

!do i=1,n
!write(6,*) (EdgesArray(i,j),j=1,n)
!end do

!----------------------------------------------------------------------------------------!
! (Step 1) Data preparartion

ReducedEdgesArray = 0

k = 0
do i=2,n
do j=1,i-1
if( EdgesArray(i,j) > 0 )then
k=k+1
ReducedEdgesArray(k,1) = i
ReducedEdgesArray(k,2) = j
ReducedEdgesArray(k,3) = EdgesArray(i,j)
end if
end do
end do

!write(6,*) "Number of edges: ", k

!----------------------------------------------------------------------------------------!
! (Step2) Sort the edges in order of increasing weight 

do i=1,k
do j=1,k
if( ReducedEdgesArray(i,3) < ReducedEdgesArray(j,3) )then
tmp_01 = ReducedEdgesArray(i,1)
tmp_02 = ReducedEdgesArray(i,2)
tmp_03 = ReducedEdgesArray(i,3)
ReducedEdgesArray(i,1) = ReducedEdgesArray(j,1)
ReducedEdgesArray(i,2) = ReducedEdgesArray(j,2)
ReducedEdgesArray(i,3) = ReducedEdgesArray(j,3)
ReducedEdgesArray(j,1) = tmp_01
ReducedEdgesArray(j,2) = tmp_02
ReducedEdgesArray(j,3) = tmp_03
end if
end do
end do

!do i=1,k
!write(6,*) i, ReducedEdgesArray(i,1),ReducedEdgesArray(i,2),ReducedEdgesArray(i,3)
!end do

!----------------------------------------------------------------------------------------!
! (Step3) Kruskal's algorithm 

MST = 0 ! Initialization
NbEdges = 1
k = 1
MST(NbEdges,1) = ReducedEdgesArray(k,1) ! Min. Spanning Tree MST (First Edge)
MST(NbEdges,2) = ReducedEdgesArray(k,2) ! Min. Spanning Tree MST (First Edge)
MST(NbEdges,3) = ReducedEdgesArray(k,3) ! Min. Spanning Tree MST (First Edge)

LV = 0 ! Initialization
LV( int(MST(NbEdges,1)), int(MST(NbEdges,2)) ) = 1 ! Edges added in first passe
LV( int(MST(NbEdges,2)), int(MST(NbEdges,1)) ) = 1 ! Edges added in first passe
do i = 1, n
LV(i,i) = 1
end do

do while ( NbEdges < (n-1) ) 

k = k + 1

!write(6,*) 'k',k

!----- Check if Vertices are already connected 

FormACycle = .FALSE.
i1 = int( ReducedEdgesArray(k,1) )
j1 = int( ReducedEdgesArray(k,2) )
if( LV(i1,j1) == 1 ) FormACycle = .TRUE.

!----- Add an edge to the MST as long as 
!----- it does not form a cycle with edges added in previous passes

!write(6,*) FormACycle, ReducedEdgesArray(k,1), ReducedEdgesArray(k,2)

if( FormACycle .eqv. .FALSE. )then

  NbEdges = NbEdges + 1

  !----- Min. Spanning tree (Add an Edge)
  
  MST(NbEdges,1) = ReducedEdgesArray(k,1)
  MST(NbEdges,2) = ReducedEdgesArray(k,2)
  MST(NbEdges,3) = ReducedEdgesArray(k,3)
  
  !----- Edges added in previous passes
  
  LV( int(MST(NbEdges,1)), int(MST(NbEdges,2)) ) = 1 ! Edges added in previous passes
  LV( int(MST(NbEdges,2)), int(MST(NbEdges,1)) ) = 1 ! Edges added in previous passes

  FindNewEdge = .TRUE.
  do while ( FindNewEdge )
  FindNewEdge = .FALSE.
  do i = 1, n
	do j = 1, n
	if ( j .ne. i )then
	if ( LV(i,j) .eq. 1 )then
	  do j1 = 1, n
	  if( j1 .ne. j )then
	  if( LV(j,j1) .eq. 1 )then
	  if( LV(i,j1) .eq. 0 )then
	  LV(i,j1) = 1
	  FindNewEdge = .TRUE.
	  end if 
	  end if
	  end if
	  end do   
	end if
	end if
	end do
  end do
  end do

!do i = 1, n
!write(6,*) (LV(i,j),j=1,n)
!end do

end if

end do

write(6,*) "Minimum Spanning Tree (MST) Found: "

do i=1, n-1
write(6,*) i, MST(i,1), MST(i,2), MST(i,3)
end do

end program kruskal_example

Algorithm to find the Legendre polynomial expansion of a phase function (fortran 90)

					
program main

implicit none 

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k

real(pr) :: a

integer :: n

integer :: ierr, nb_lines

real(pr), dimension(:), allocatable :: angle,p11,p12,p22,p33,p43,p44

real(pr), dimension(:,:), allocatable :: pl

real(pr) :: pl_n2, pl_n1

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: cdr

real(pr) :: x

real(pr), dimension(:,:), allocatable :: MA, ML
real(pr), dimension(:), allocatable :: MB, MY, MX, p11_lp

cdr = pi / 180.0

!----------------------------------------------------------------------------------------!
! user input 

n = 120 ! Number of Legendre Polynomial Coefficients 

open(1,file='hexagomal_pristine_crystal_GOM.txt', &
    & form="formatted",iostat=ierr,status="old")
	do i = 1, 18
		read(1,*)
	end do
	if (ierr/=0) then
		write(*,*) "Impossible to open the file: modis_associated_files.txt"
		stop
	endif
	nb_lines = 0
	do while (ierr==0)
		read(1,*,iostat=ierr)
		if (ierr==0) nb_lines = nb_lines + 1
	enddo
close(1)

!write(6,*) 'nb_lines', nb_lines

allocate( angle( nb_lines ) )
allocate( p11( nb_lines ) )
allocate( p12( nb_lines ) )
allocate( p22( nb_lines ) )
allocate( p33( nb_lines ) )
allocate( p43( nb_lines ) )
allocate( p44( nb_lines ) )

open(1,file='hexagomal_pristine_crystal_GOM.txt')
	do i = 1, 18
		read(1,*)
	end do
	do i = 1 , nb_lines
		read(1,*) angle(i),p11(i),p12(i),p22(i),p33(i),p43(i),p44(i)
		!write(6,*) angle(i)
	enddo
close(1)

!----------------------------------------------------------------------------------------!
! Legendre Polynomials 

allocate( pl( n, nb_lines ) )

do i = 1, nb_lines
	x = dcos( angle(i) * cdr )
    pl(1,i) = 1.0
    pl(2,i) = x
    pl_n2 = pl(1,i)
    pl_n1 = pl(2,i)
    do j = 3 , n
        pl(j,i) = ( 1.0 / j ) * ( ( 2.0 * j - 1.0 ) * x * pl_n1 - ( j - 1.0 ) * pl_n2 )
        pl_n2 = pl_n1
        pl_n1 = pl(j,i)
    end do
end do

!----------------------------------------------------------------------------------------!
! A matrix

allocate( MA(n,n) )

do i = 1 , n
	do j = 1 , n
	a = 0.0
		do k = 1, nb_lines
			a = a + ( pl(i,k) * pl(j,k) ) / ( p11(k)**2 )
		enddo
	MA(i,j) = a
	enddo
enddo

!----------------------------------------------------------------------------------------!
! B matrix

allocate( MB(n) )

do i = 1 , n
	a = 0.0
	do k = 1, nb_lines
		a = a + pl(i,k) / p11(k)
	enddo
	MB(i) = a 
enddo

!----------------------------------------------------------------------------------------!
! Choleski Factorisation 

allocate( ML(n,n) )

ML(1,1) = sqrt( MA(1,1) )
do i = 2 , n
	ML(i,1) = MA(1,i) / ML(1,1)
enddo

do j = 2 , n
	do i = 1 , n
		if( i .lt. j ) then
			ML(i,j) = 0.0
		else
			if( i .eq. j) then
				a = 0.0
				do k = 1 , i -1
					a = a + ML(i,k)**2
				enddo
				ML(i,j) = sqrt( MA(i,i) - a )
			else
				a = 0.0
				do k = 1 , i - 1
					a = a + ML(j,k) * ML(i,k)
				enddo
				ML(i,j) = ( MA(j,i) - a ) / ML(j,j)
			endif
		endif
	enddo
enddo

!----------------------------------------------------------------------------------------!
! Solve linear system 1

allocate( MY(n) )

MY(1) = MB(1) / ML(1,1)
do i = 2 , n
	a = 0.0
	do k = 1 , i-1
		a = a + ML(i,k) * MY(k)
	enddo
	MY(i) = ( MB(i) - a ) / ML(i,i)
enddo

!----------------------------------------------------------------------------------------!
! Solve linear system 2

allocate( MX(n) )

MX(n) = MY(n) / ML(n,n)
do i = n-1, 1, -1
	a = 0.0
	do k = n, i+1, -1
		a = a + ML(k,i) * MX(k)
	enddo
	MX(i) = ( MY(i) - a ) / ML(i,i)
enddo

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

allocate( p11_lp(nb_lines) )

do i = 1, nb_lines
	a = 0
    do j = 1 , n
    	a = a + MX(j) * pl(j,i) 
    end do
	p11_lp(i) = a
end do

!----------------------------------------------------------------------------------------!
! Save results

open(1,file='phase_function_legendre_polynomial_expansion_01.txt')
do i = 1 , n
	MX(i) = MX(i) / ( 2.0d0 * i + 1.0d0 )
	write(1,*) MX(i) / MX(1)
	!write(6,*) MX(i) !/ MX(1)
enddo
close(1)

open(1,file='phase_function_legendre_polynomial_expansion_02.txt')
do i = 1 , nb_lines
	write(1,*) angle(i), p11_lp(i), p11(i)
enddo
close(1)

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

end program main

Plot CALIOP 01&05km CLay cloud phase vertical profile with python and matplotlib

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