atb_calipso

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from pyhdf.SD import SD, SDC
def extraction(fichier):
	datasd = SD(fichier, SDC.READ)
	lg = datasd.select('Longitude')[:]
	lt = datasd.select('Latitude')[:]
	ti = datasd.select('Profile_UTC_Time')[:,:]
	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('lidaralt.asc', 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('lidaralt.asc', 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

	atb = atb[idx,:] #erreur, si pas : le tableau devient (_,), si : le tableau (_,583), probleme absent de goccp
	n = len(lt6090)
	return n, lt6090, lg6090, ti6090,atb 
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)

Pour commencer : c'est CALIPSO avec un I, pas un Y. :-)

Le probleme etait un de ceux qu'on avait deja vus : les vecteurs lus dans un fichier HDF ont une shape de type (56205, 1), au lieu de (56205,) tout court. Il y a une dimension en trop.

Du coup, quand tu fais idx = ( lt < -60 ) & ( lt >= -90 ), le vecteur idx a la meme shape (56205,1). Ca ne pose pas de probleme pour indexer lon, lat et temps qui ont cette shape-la aussi, mais ca merde quand tu essayes d'indexer atb qui a des dimensions differentes.

Version qui marche : https://gist.github.com/vnoel/11401464

Le probleme n'arrive pas quand on lit les fichiers GOCCP qui sont en netCDF : quand tu lis un vecteur dans un fichier netCDF il n'y a pas de dimension en trop, le vecteur a les bonnes dimensions (56205,).