autom

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import os as os

def listevar(fichier):
	datago = Dataset(fichier, 'r')
	x = datago.variables.keys()
	print x

def extraction(fichier):
	datago = Dataset(fichier, 'r')
	lg = datago.variables['longitude'][:]
	lt = datago.variables['latitude'][:]
	ti = datago.variables['time'][:]
	atb = datago.variables['ATB'][:,:]
	SR = datago.variables['instant_SR'][:,:]
	sol_alt = datago.variables['SE'][:]
	alt = datago.variables['alt_bound'][:,:]
	print np.shape(alt)
	u, v = np.shape(alt)
	alt_bas = alt[0,:]
	alt_haut = alt[1,:]
	alt = np.zeros((v,1))
	for i in range(0,v):
		alt[i] = alt_bas[i]
	#alt[v] = alt_haut[v-1]
	alt = alt[:,0]
	return lg, lt, ti, atb,alt, sol_alt, SR


def affiche(n1,n2,lg,lt,ti,atb,alt,sol_alt):
	#lg = lg[n1:n2]
	#lt = lt[n1:n2]
	ti = ti[n1:n2]
	sol_alt = sol_alt[n1:n2]
	u,v = np.shape(atb[n1:n2,:])
	trait = ti[int (((n2-n1)/2))] 
	plt.plot ([trait,trait] , [alt[0],alt[-1]],'w--',linewidth=2)
	plt.plot(ti, sol_alt,'w',linewidth=2)


	#print 'tailles des tableaux',ti.shape, alt.shape, atb.shape
	plt.pcolormesh(ti,alt,atb[n1:n2,:].T)
	plt.ylim(-1, 16)
	plt.xlim(ti[0], ti[-1])
	plt.clim(0, 0.03)
	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,n2,lg,lt,ti,atb,alt):
	#print 'tailles des tableaux',ti.shape, alt.shape, atb.shape
	n = int (((n2-n1)/2))
	plt.semilogx(atb[n,:], alt)
	#plt.xlim(1e-7, 1e-1)
	plt.ylim(-20,20)
	plt.xlabel('ATB [km-1 sr-1]')
	plt.ylabel('Altitude [km]')
	plt.title('Profil %d' % (n))
	plt.show()

def profil_SR(n1,n2,lg,lt,ti,atb,alt,fichier):
	#print 'tailles des tableaux',ti.shape, alt.shape, atb.shape
	n = int (((n2-n1)/2))
	plt.plot(atb[n,:], alt)
	plt.xlim(-0.1, 30)
	plt.ylim(-20,20)
	plt.xlabel('ATB [km-1 sr-1]')
	plt.ylabel('Altitude [km]')
	plt.title('Profil %d' % (n))
	#plt.show()
	path = fichier
	path += '.png'
	plt.savefig(path)
	plt.close()
	

def coord_pole(lg,lt,ti,atb, sol_alt):
	idx = ( lt < -60 ) & ( lt >= -90 )
	lt6090 = lt[idx]
	lg6090 = lg[idx]
	ti6090 = ti[idx]
	sol_alt = sol_alt[idx]
	atb = atb [idx,:]
	n = len(lt6090)
	return n, lt6090, lg6090, ti6090,atb, sol_alt

def alt_max(SR, alt, n1, n2):
	u, v = np.shape(SR)
	profil = int (((n2-n1)/2))
	#j = -1
	#for i in range(0,v):
	#	if SR[profil,i] < 0.1 :	
	#		if i > 0:
	#			if ( i - j == 1) :
	#				h_max = alt[i]
	#		j = i
	j = - 1
	i = 0
	while SR[profil,i] < 0.01 and i < v and SR[profil,i] > -10:
		j = i
		i = i + 1
		#if SR[profil,i] < 0 : print profil, i,SR[profil,i], alt[i]
	return j

def alt_max_orbite():
	lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
	p,v = np.shape(atb)
	tab1 = np.zeros((p))
	tab2 = np.zeros((p))
	for o in range (0,p):
		j = alt_max(SR,alt,0,2*o)
		if j != -1 :
			tab1[o] = alt[j]
			#if tab1[o] > 5.: print o
		else : tab1[o] = -10
		tab2[o] = o

	tab1 = np.ma.masked_where(tab1 < 1, tab1)
	plt.plot(tab2,tab1)
	plt.show()

def alt_max_profil(alt,SR,n1,n2,fichier):
	j = alt_max(SR,alt,n1,n2)
	s = np.shape(alt)
	if j != -1 :alt_ext= alt[j]
	else : alt_ext= -10
	#for i in range(0,40):
	#	alt[i] = alt[i] + alt_ext
	print alt_ext
	profil_SR(n1,n2,lg,lt,ti,SR,alt,fichier)
	
def tr_fl(alt,SR,n1,n2):
	j = alt_max(SR,alt,n1,n2)
	profil = (n2-n1)/2
	u = max(np.shape(alt))
	if j != -1 :
		j = j + 1
		b = 0
		while j < u :
			if SR[profil,j] < 0.01 and  SR[profil,j] > -10: 
				b = b +1
			j = j + 1
		if b != 0 : print profil , 'faux eteint'
		else : print profil,'vrai eteint'
	else : print profil,' ????? '

def tr_fl_orbite(alt,SR):	
	u,v = np.shape(SR)
	n1 = 0
	for i in range (0,max(u)):
		tr_fl(alt,SR,n1,2*i)

def lecture_dossier(repertoire):
	u = np.shape(os.listdir(repertoire))
	
	return u, os.listdir(repertoire)

if __name__=='__main__':

	#fichier3="/bdd/CFMIP/CFMIP_OBS_LOCAL/GOCCP/instant_SR_CR_DR/grid_L40/2013/201308/night/instant_SR_CR_DR_2013-08-24T23-19-30ZN_night_CFMIP2_2.70.nc"#38418 mid
	fichier3="/bdd/CFMIP/CFMIP_OBS_LOCAL/GOCCP/instant_SR_CR_DR/grid_L40/2013/201308/night/instant_SR_CR_DR_2013-08-23T19-18-36ZN_night_CFMIP2_2.70.nc" #20632 trop
	listevar(fichier3)
	n1 , n2 = 0  , 20632
	lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
 	#alt_max(SR, alt, n1,2*n2)
	#profil_SR(n1,2*n2,lg,lt,ti,SR,alt)
	#affiche(n1,n2,lg, lt, ti, SR,alt,sol_alt)
	'''
	n1 = 0
	n2, lt, lg, ti, atb, sol_alt = coord_pole(lg, lt, ti, atb,sol_alt)
	alt_max(SR, alt, n1,n2)
	profil(n1,n2,lg,lt,ti,atb,alt)	
	affiche(n1,n2,lg, lt, ti, atb,alt,sol_alt)
	lonlat(n1, n2,lg, lt, ti, atb)
	'''

	#lonlat(n1, n2,lg, lt, ti, atb)
	#affiche(n1,n2,lg, lt, ti, SR,alt,sol_alt)
	#alt_max_orbite()
	#alt_max_profil(alt,SR,n1,2*n2)
	dossier = '/bdd/CFMIP/CFMIP_OBS_LOCAL/GOCCP/instant_SR_CR_DR/grid_L40/2013/201308/night/'
	nbr_fichier,fichier = lecture_dossier(dossier)
	for i in range(0, 12):
		n1 , n2 = 0  , 20632
		fichier3 = dossier + fichier[i]
		if i > 0 :
			u,v = np.shape((SR))
			SR = np.zeros((u,v))
		lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)		
		alt_max_profil(alt,SR,n1,2*n2,fichier[i])