atb_calipso.py

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from pyhdf.SD import SD, SDC

altfile = '/users/noel/pycode/staticdata/lidaralt.asc'

def extraction(fichier):
	datasd = SD(fichier, SDC.READ)
	lg = datasd.select('Longitude')[:,0]
	lt = datasd.select('Latitude')[:,0]
	ti = datasd.select('Profile_UTC_Time')[:,0]
	atb = datasd.select('Total_Attenuated_Backscatter_532')[:,:]
	return lg, lt, ti, atb

def affiche(n1,n2,lg,lt,ti,atb):
	#lg = lg[n1:n2]
	#lt = lt[n1:n2]
	print ' hehe' ,ti.shape
	ti = ti[n1:n2]
	alt = np.loadtxt(altfile, delimiter =',')
	print 'tailles des tableaux',ti.shape, alt.shape, atb.shape
	plt.pcolormesh(ti,alt,atb[n1:n2,:].T)
	plt.ylim(0, 20)
	plt.xlim(ti[0], ti[-1])
	plt.clim(0, 5e-3)
	plt.xlabel('Time')
	plt.ylabel('Altitude [km]')
	plt.colorbar()	  
	plt.show() 


def lonlat(n1, n2,lg,lt,ti,atb):
	from mpl_toolkits.basemap import Basemap
	m = Basemap()
	m.plot(lg, lt)
	lg = lg[n1:n2]
	lt = lt[n1:n2]
	m.plot(lg, lt, 'r')
	m.drawcoastlines()
	m.fillcontinents()
	plt.show()


def profil(n1,lg,lt,ti,atb):
	alt = np.loadtxt(altfile, delimiter =',')
	print 'tailles des tableaux',ti.shape, alt.shape, atb.shape
	plt.semilogx(atb[n1,:], alt)
	plt.xlim(1e-7, 1e-2)
	plt.ylim(0,20)
	plt.xlabel('ATB [km-1 sr-1]')
	plt.ylabel('Altitude [km]')
	plt.title('Profil %d' % (n1))
	plt.show()

def coord_pole(lg,lt,ti,atb):
	idx = ( lt < -60 ) & ( lt >= -90 )
	lt6090 = lt[idx]
	lg6090 = lg[idx]
	ti6090 = ti[idx]
	print 'atb en entree', atb.shape

	print atb.shape, idx.shape

	atb6090 = atb[idx,:] #erreur, si pas : le tableau devient (_,), si : le tableau (_,583), probleme absent de goccp
	n = len(lt6090)
	return n, lt6090, lg6090, ti6090,atb6090

if __name__=='__main__':

	fichier3="/bdd/CALIPSO/Lidar_L1/CAL_LID_L1.v3.30/2013/2013_08_24/CAL_LID_L1-ValStage1-V3-30.2013-08-24T23-19-30ZN.hdf"
	n1, n2 = 0, 100
	lg, lt, ti, atb = extraction(fichier3)
	#profil(n1,lg,lt,ti,atb)
	#affiche(n1,n2,lg, lt, ti, atb)
	#lonlat(n1,n2,lg, lt, ti, atb)
	n2, lt, lg, ti, atb = coord_pole(lg, lt, ti, atb)
	profil(n1,lg,lt,ti,atb)
	affiche(n1,n2,lg, lt, ti, atb)