Skip to content

Instantly share code, notes, and snippets.

@stageipsl
Last active August 29, 2015 14:01
Show Gist options
  • Save stageipsl/4edfee9fee8a31c4b621 to your computer and use it in GitHub Desktop.
Save stageipsl/4edfee9fee8a31c4b621 to your computer and use it in GitHub Desktop.
#!/usr/bin/env python
import numpy as np
import numpy.matlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from netCDF4 import Dataset
import os as os
import resource
from mpl_toolkits.basemap import Basemap
from scipy.stats import norm
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'][:,:]
datago.close()
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):
lt = lt[n1:n2]
sol_alt = sol_alt[n1:n2]
u,v = np.shape(atb[n1:n2,:])
trait = lt[int (((n2-n1)/2))]
plt.plot ([trait,trait] , [alt[0],alt[-1]],'w--',linewidth=2)
plt.plot(lt, sol_alt,'w',linewidth=2)
plt.pcolormesh(lt,alt,atb[n1:n2,:].T)
plt.ylim(-1, 16)
plt.xlim(lt[0],lt[-1])
#plt.gca().invert_xaxis()
plt.clim(0, 0.005)
plt.xlabel('Latitude')
plt.ylabel('Altitude [km]')
plt.colorbar()
plt.show()
def afficheSR(lg,lt,ti,atb,alt,sol_alt):
n1 = 0
n2,x = np.shape(atb)
lt = lt[n1:n2]
sol_alt = sol_alt[n1:n2]
u,v = np.shape(atb[n1:n2,:])
sol_alt = sol_alt[n1:n2]
fig = plt.figure(figsize=(9, 6), dpi=100)
plt.plot(lt, sol_alt,'w',linewidth=2)
plt.pcolormesh(lt,alt,atb[n1:n2,:].T)
plt.gca().invert_xaxis()
plt.ylim(-1, 16)
plt.clim(1e-5, 1e-1)
plt.xlabel('Latitude')
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 draw():
m = Basemap()
lon1 ,lon2,lat1,lat2 = -75,-90,0,15 # shallowcu
#lon1 ,lon2,lat1,lat2 = 180,70,-15,15 # warmpool
fig = plt.figure(figsize=(14, 7), dpi=100)
m.plot([lon1,lon1], [lat1,lat2], lw = 2, color = 'Gold')
m.plot([lon2,lon2], [lat1,lat2], lw = 2, color = 'Gold')
m.plot([lon1,lon2], [lat1,lat1], lw = 2, color = 'Gold')
m.plot([lon1,lon2], [lat2,lat2], lw = 2, color = 'Gold')
m.bluemarble()
parallels = np.arange(-90.,90,20.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10,color = 'Chocolate')
meridians = np.arange(-180.,180.,20.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10, color = 'Chocolate')
#m.drawcoastlines()
#m.fillcontinents()
#titre = 'Warmpool'
titre = 'Shallowcu'
path = titre + '_carte'
plt.title(titre)
plt.savefig(path)
#plt.show()
plt.close()
def profil(n,lg,lt,ti,atb,alt,fichier):
datago = Dataset(fichier, 'r')
atmol = datago.variables['ATB_mol'][:]
p1, = plt.semilogx(atb[n,:], alt)
p2, = plt.semilogx(atmol[n,:], alt)
plt.legend([p1,p2], ["ATB mesure","ATB moleculaire"])
#plt.xlim(1e-7, 1e-1)
plt.ylim(0,20)
plt.xlabel('ATB [km-1 sr-1]')
plt.ylabel('Altitude [km]')
plt.title('Profil ' + str(lt[n]))
plt.show()
def profil_SR(n,lg,lt,ti,atb,alt,fichier):
plt.plot(atb[n,:], alt)
plt.xlim(-0.1, 50)
plt.ylim(0,20)
plt.xlabel('SR')
plt.ylabel('Altitude [km]')
plt.title('Profil ' + str(lt[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, n):
u, v = np.shape(SR)
j = - 1
if (SR[n,0] > 0.01 or SR[n,0] < -10) or (SR[n,1] > 0.01 or SR[n,1] < -10):
return j
i = 2
while i < v:
if (SR[n,i] > -10) and (SR[n,i] < 0.01):
j = i
elif (SR[n,i] > 0.01) and (SR[n,i-1] > 0.01):
return j
i = i + 1
return j
def alt_max_orbite(alt, SR):
p,v = np.shape(SR)
tab1 = np.zeros((p))
tab2 = np.zeros((p))
for o in range (0,p):
j = alt_max(SR,alt,o)
if j != -1 :
tab1[o] = alt[j]
#if tab1[o] > 5.: print o
else : tab1[o] = -10
tab2[o] = o
return tab1
'''
tab1 = np.ma.masked_where(tab1 <= -1, tab1)
fig = plt.figure(figsize=(15, 10), dpi=300)
plt.plot(lt,tab1)
plt.plot(lt, sol_alt,'g',linewidth=2)
plt.xlim(np.min(lt), np.max(lt))
plt.ylabel('Altitude')
plt.xlabel('indice')
plt.show()
'''
def alt_max_orbite_flag(alt, SR):
u,v= np.shape(SR)
tab1 = np.zeros((u,v))
for o in range (0,u):
j = alt_max(SR,alt,o)
if j != -1 :tab1[:][o] = alt[j]
else : tab1[:][o] = -10
return tab1
def alt_min_SR5(alt_array,SR,altmax):
u,v = np.shape(SR)
alt_SR = np.zeros((u,v))
idx = (SR > 5)
alt_copy = alt_array.copy()
alt_copy[~idx] = 9999.
alt_copy[altmax < 0] = 9999.
alt_SR = np.min(alt_copy, axis=1)
alt_SR[alt_SR == 9999] = -1.
return alt_SR
def alt_min_SR5_orbite(alt_array,SR,altmax):
u,v = np.shape(SR)
alt_SR = np.zeros((u,v))
alt_tab = np.zeros((u,v))
alt_total = np.zeros((u,v))
idx = (SR > 5)
alt_copy = alt_array.copy()
alt_copy[~idx] = 9999.
alt_copy[altmax < 0] = 9999.
alt_SR = np.min(alt_copy, axis=1)
alt_SR[alt_SR == 9999] = -1.
for i in range(0,v-1):
if i == 0: alt_tab = alt_SR
alt_tab = np.vstack((alt_tab,alt_SR))
alt_tab = alt_tab.T
print alt_tab
print np.shape(alt_tab)
return alt_tab
def alt_max_profil(alt,SR,n,fichier):
j = alt_max(SR,alt,n)
s = np.shape(alt)
if j != -1 :alt_ext= alt[j]
else : alt_ext= -10
print alt_ext
profil_SR(n,lg,lt,ti,SR,alt,fichier)
def affiche_flag(lg,lt,atb,alt,sol_alt):
n1 = 0
n2, x = np.shape(atb)
lt = lt[n1:n2]
sol_alt = sol_alt[n1:n2]
u,v = np.shape(atb[n1:n2,:])
fig = plt.figure(figsize=(14, 7), dpi=100)
plt.plot(lt, sol_alt,'Maroon', alpha=1)
cmap = cm.colors.ListedColormap(['White', 'Navy', 'Lime', 'Chocolate','Gold'])
plot = plt.pcolormesh(lt,alt,atb.T,cmap=cmap)
plt.ylim(0, 18)
plt.xlim(-80, 80)
plt.clim(0, 4)
plt.xlabel('Latitude')
plt.ylabel('Altitude en km')
cb =plt.colorbar(ticks=[0, 1.3, 2, 2.75,3.5])
cb.ax.set_yticklabels(['Clear Sky','Nuages','Faux eteint','Zone de transition','Vrai eteint'])
path = 'SR_orbite_vrai_faux'
plt.savefig(path)
plt.close()
def affiche_flag_grl(lg,lt,atb,alt,sol_alt):
n1 = 0
n2, x = np.shape(atb)
lt = lt[n1:n2]
sol_alt = sol_alt[n1:n2]
u,v = np.shape(atb[n1:n2,:])
fig = plt.figure(figsize=(14, 7), dpi=100)
plt.plot(lt, sol_alt,'Maroon', alpha=1)
cmap = cm.colors.ListedColormap(['White', 'Navy', 'Chocolate'])
plot = plt.pcolormesh(lt,alt,atb.T,cmap=cmap)
plt.ylim(0, 18)
plt.xlim(-80, 80)
plt.clim(0, 2)
plt.xlabel('Latitude')
plt.ylabel('Altitude en km')
cb =plt.colorbar(ticks=[0, 1.3, 2])
cb.ax.set_yticklabels(['Clear Sky','Nuages','Eteint'])
path = 'SR_orbite_GRL'
plt.savefig(path)
plt.close()
def flag_orbite_grl(alt,SR):
u,v = np.shape(SR)
flag = np.zeros((u,v))
alt_array = np.zeros((u,v))
alt_SR = np.zeros((u,v))
alt_array[:][:] = alt
alt_max = alt_max_orbite_flag(alt,SR)
alt_SR = np.matlib.repmat(alt_min_SR5(alt_array,SR,alt_max),v,1).T
print np.shape(alt_SR)
#Nuage
idx = (SR > 5)
flag[idx] = 1
#Vrai eteint
idx4 = (SR < 0.01) & (SR> -10)
flag[idx4] = 2
#Zone de transition
#idx3 = (alt_array > alt_max) & (alt_array < alt_SR)
#flag[idx3] = 3
return flag
def flag_orbite(alt,SR):
u,v = np.shape(SR)
flag = np.zeros((u,v))
alt_array = np.zeros((u,v))
alt_SR = np.zeros((u,v))
alt_array[:][:] = alt
alt_max = alt_max_orbite_flag(alt,SR)
alt_SR = np.matlib.repmat(alt_min_SR5(alt_array,SR,alt_max),v,1).T
print np.shape(alt_SR)
#Nuage
idx = (SR > 5)
flag[idx] = 1
#Vrai eteint
idx4 = (alt_array < alt_max)&(SR < 0.01) & (SR> -10)
flag[idx4] = 4
#Faux eteint
idx2 = ((alt_array > alt_max) & (alt_max!= -10)&(SR < 0.01) & (SR> -10)) | ((alt_max == -10) & (SR < 0.01) & (SR> -10))
flag[idx2] = 2
#Zone de transition
idx3 = (alt_array > alt_max) & (alt_array < alt_SR)
flag[idx3] = 3
return flag
def lecture_dossier(repertoire):
u = np.max(np.shape(os.listdir(repertoire)))
return u, os.listdir(repertoire)
def occur(flag,alt):
clouds = (flag == 1).sum(axis=0)
faux_eteint = (flag == 2).sum(axis=0)
cloud_opa = (flag == 3).sum(axis=0)
vrai_eteint = (flag == 4).sum(axis=0)
return clouds, faux_eteint, cloud_opa, vrai_eteint
def freq_occur(flag,alt):
u,v = np.shape(flag)
clouds, faux_eteint, cloud_opa, vrai_eteint = occur(flag,alt)
clouds, faux_eteint, cloud_opa, vrai_eteint = (1.)*clouds / u, (1.)*faux_eteint / u, (1.)*cloud_opa / u, (1.)*vrai_eteint / u
fig = plt.figure(figsize=(9, 15), dpi=100)
p1, = plt.plot(vrai_eteint,alt)
sum1 = vrai_eteint + cloud_opa
p2, = plt.plot(sum1,alt)
sum2 = vrai_eteint + cloud_opa + clouds
p3, = plt.plot(sum2,alt)
sum3 = vrai_eteint + cloud_opa + clouds + faux_eteint
p4, = plt.plot(sum3,alt)
#plt.xlim(0.0001, 0.1)
plt.legend([p1,p2,p3,p4], ["Nuages","Faux eteint","Zone de transition","Vrai eteint"], loc=1)
plt.xlabel("Frequence d'occurence")
plt.ylabel("Altitude en km")
plt.show()
def occur_in_region(flag,alt,lt,latrange):
idx = (lt < latrange[1]) & (lt > latrange[0])
flag = flag[idx]
clouds = (flag == 1).sum(axis=0)
faux_eteint = (flag == 2).sum(axis=0)
cloud_opa = (flag == 3).sum(axis=0)
vrai_eteint = (flag == 4).sum(axis=0)
return clouds, faux_eteint, cloud_opa, vrai_eteint
def freq_occur_region(dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
dic_opa_cloud = dict()
dic_cloud = dict()
dic_vrai = dict()
dic_faux = dict()
#nbr_fichier=1
nbr_fichier=np.max(nbr_fichier)
latbands = {'tropique' :[-20,20],'pole_nord' :[60,90],'pole_sud' : [-90,-60],'mid_sud' :[-50,-30],'mid_nord' :[30, 50 ]}
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
flag = flag_orbite(alt,SR)
u,v = np.shape(flag)
print fichier[i]
for region, tab_reg in latbands.items():
#flag = flag[tab_reg,:]
clouds, faux_eteint, cloud_opa, vrai_eteint = occur_in_region(flag, alt, lt, latbands[region])
if region in dic_opa_cloud:
dic_opa_cloud[region] += 1.* cloud_opa/ (u * nbr_fichier)
dic_cloud[region] += 1.* clouds/ (u * nbr_fichier)
dic_vrai[region] += 1.* vrai_eteint/ (u * nbr_fichier)
dic_faux[region] += 1.* faux_eteint/ (u * nbr_fichier)
else:
dic_opa_cloud[region] = 1.* cloud_opa/ (u * nbr_fichier)
dic_cloud[region] = 1.* clouds/ (u * nbr_fichier)
dic_vrai[region] = 1.* vrai_eteint/ (u * nbr_fichier)
dic_faux[region] = 1.* faux_eteint/ (u * nbr_fichier)
for region in latbands.keys():
affiche_freq_separee(alt,dic_opa_cloud[region],dic_cloud[region],dic_vrai[region],dic_faux[region],region,fichier,nbr_fichier)
def affiche_freq(alt,op,cl,vr,fa,region,fichier,nbr_fichier):
fig = plt.figure(figsize=(9, 15), dpi=100)
p1, = plt.plot(op,alt)
sum1 = vr + op
p2, = plt.plot(sum1,alt)
sum2 = sum1 + fa
p3, = plt.plot(sum2,alt)
sum3 = sum2 + cl
p4, = plt.plot(sum3,alt)
plt.legend([p1], ["Nuages opaques"], loc=1)
plt.legend([p1,p2,p3,p4], ["Zone de transition","Vrais eteints","Faux eteints","Nuages"], loc=1)
plt.title(region + ' Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier))
plt.xlabel("Frequence d'occurence")
plt.ylabel("Altitude en km")
path = region
path += '.png'
plt.savefig(path)
plt.close()
def affiche_freq_separee(alt,op,cl,vr,fa,region,fichier,nbr_fichier):
fig = plt.figure(figsize=(9, 15), dpi=100)
p1, = plt.plot(op,alt)
p2, = plt.plot(cl,alt)
p3, = plt.plot(vr,alt)
plt.legend([p1,p2,p3], ["Zone de transition","Nuages","Vrai eteint"], loc=1)
plt.title(region + ' Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier))
plt.xlabel("Frequence d'occurence")
plt.ylabel("Altitude en km")
path = region
path += 'occurence_zone_de_transition_nuage_separ.png'
plt.savefig(path)
plt.close()
p3, = plt.plot(vr,alt)
p4, = plt.plot(fa,alt)
plt.legend([p3,p4], ["Vrais eteints","Faux eteints"], loc=1)
plt.title(region + ' Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier))
plt.xlabel("Frequence d'occurence")
plt.ylabel("Altitude en km")
path = region
path += 'occurence_vrai_faux_separ.png'
plt.savefig(path)
plt.close()
def histo_1D_orbite(dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier )
lg, lt, ti, atb,alt, sol_alt, SR = extraction(dossier + fichier[0])
u,v = 60000,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
print np.shape(alt)
alt_array = np.zeros((u,v))
alt_array[:][:] = alt
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
print fichier3
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
u,v = np.shape(SR)
alt_array_copy = alt_array[0:u,0:v]
alt_max = alt_max_orbite_flag(alt,SR)
alt_SR5 = alt_min_SR5(alt_array_copy,SR,alt_max)
idx = ( alt_SR5 > 0 )
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: total = temp
else : total = np.concatenate([total, temp])
print np.shape(total)
fig = plt.figure(figsize=(14, 7), dpi=100)
plt.hist(total, bins=np.arange(0,14,0.5))
(mu, sigma) = norm.fit(total)
plt.title('Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier)+ ' \n Moyenne : ' + str(mu) + ' Ecart type :' + str(sigma))
plt.xlabel('Taille de la zone de transition en km')
plt.ylabel("Nombre d'occurences")
plt.show()
path = 'globale_distri_taille_transition.png'
plt.savefig(path)
plt.close()
def histo_2D(alt, SR,lt,nbr_profil,profil_debut):
n1 = profil_debut
n2 = profil_debut + nbr_profil -1
u,v = np.shape(SR)
alt_array = np.zeros_like(SR)
alt_array[:][:] = alt
fig = plt.figure(figsize=(14, 7))
SRbin=np.arange(-1,40)
altbin = np.arange(-0.240,18.71+0.240,0.480)
x = np.concatenate([SR[n1,:], SR[n2,:]], axis=0)
y = np.concatenate([alt_array[n1,:], alt_array[n2,:]], axis=0)
#x = SR.ravel()
#y = alt_array.ravel()
print 'hehe',np.shape(x),n1,n2, x
plt.hist2d(x, y,bins=(SRbin,altbin))
if nbr_profil == 2 :
plt.plot(SR[n1,:], alt_array[0,:],'-*k',linewidth=2)
plt.plot(SR[n2,:], alt_array[0,:],'-+y',linewidth=2)
plt.xlim(-5,40)
plt.xlabel('SR')
#plt.clim(0,300)
plt.ylabel('Altitude en km')
if nbr_profil ==1 : plt.title('Latitude :'+ str( lt[n1]))
else : plt.title('Latitude '+ str(lt[n1])+' a '+ str(lt[n2]))
plt.colorbar()
plt.show()
def histo_2D_orbite(alt, SR,o,p, test):
alt_array = np.zeros_like(SR)
alt_array[:][:] = alt
u, v = np.shape(SR)
if test == 0:
n1 = 0
n2 =u * (v)
elif test == 1:
n1 = o * (v)
n2 = p * (v)
alt_max = alt_max_orbite_flag(alt, SR)
alt_minSR =alt_min_SR5_orbite(alt_array,SR,alt_max)
idx = (alt_max != -10)#(alt_array <= alt_max) & (alt_array >= alt_minSR) & (alt_minSR != -1)
SR = SR[idx]
alt_array = alt_array[idx]
print 'SR',np.shape(SR)
fig = plt.figure(figsize=(14, 7))
#SRb=np.arange(-1,40,1)
SRb=np.concatenate([np.arange(-1,3,0.2),np.arange(3,40,1)])
binwidth = np.max(np.shape(SRb))
altbin = np.arange(-0.240,18.71+0.240,0.480)
#altbin = np.linspace(-0.240,19,binwidth)
#x = np.concatenate([SR[0,:], SR[u-1,:]], axis=0)
#y = np.concatenate([alt_array[0,:], alt_array[u-1,:]], axis=0)
x = SR.ravel()
y = alt_array.ravel()
plt.hist2d(x[n1:n2], y[n1:n2],bins=(SRb,altbin))
plt.xlabel('SR')
plt.clim(0,800)
plt.ylabel('Altitude en km')
titre = ('Histogramme 2D de la latitude '+str(lt[o])+' a '+str(lt[p]))
plt.title(titre)
plt.colorbar()
plt.show()
def histo_1D_region2(dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier)
lg, lt, ti, atb,alt, sol_alt, SR = extraction(dossier + fichier[0])
u,v = 60000,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
print np.shape(alt)
alt_array = np.zeros((u,v))
alt_array[:][:] = alt
latbands = {'tropique' :[-20,20],'pole_nord' :[60,90],'pole_sud' : [-90,-60],'mid_sud' :[-50,-30],'mid_nord' :[30, 50 ]}
dic_hist = dict()
jour = np.zeros((nbr_fichier),dtype = int)
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
jour[i] = int(fichier[i][25:27])
print jour[i], fichier[i]
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
u,v = np.shape(SR)
alt_array_copy = alt_array[0:u,0:v]
print np.shape(alt_array_copy),np.shape(SR)
print i, '/', nbr_fichier
#idx1 = (lt < latbands[region][1]) & (lt > latbands[region][0])
alt_max = alt_max_orbite_flag(alt,SR)
alt_SR5 = alt_min_SR5(alt_array_copy,SR,alt_max)
idx = ( alt_SR5 > 0 ) & (lt < 20) & (lt > -20)
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist['tropique'] = temp
else : dic_hist['tropique'] = np.concatenate([dic_hist['tropique'], temp])
idx = ( alt_SR5 > 0 ) & (lt < 90) & (lt > 60)
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist['pole_nord'] = temp
else : dic_hist['pole_nord'] = np.concatenate([dic_hist['pole_nord'], temp])
idx = ( alt_SR5 > 0 ) & (lt < -60) & (lt >-90)
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist['pole_sud'] = temp
else : dic_hist['pole_sud'] = np.concatenate([dic_hist['pole_sud'], temp])
idx = ( alt_SR5 > 0 ) & (lt < -30) & (lt > -50)
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist['mid_sud'] = temp
else : dic_hist['mid_sud'] = np.concatenate([dic_hist['mid_sud'] , temp])
idx = ( alt_SR5 > 0 ) & (lt < 50) & (lt > 30)
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist['mid_nord'] = temp
else : dic_hist['mid_nord'] = np.concatenate([dic_hist['mid_nord'] , temp])
for region, tab in dic_hist.items():
print tab
(mu, sigma) = norm.fit(tab)
plt.hist(tab, bins=np.arange(0,14,0.5))
plt.rcParams["axes.titlesize"] = 8
plt.title(region+' Du ' + fichier[0][17:25]+str(np.min(jour)) +' au ' + fichier[nbr_fichier][17:25] +str(np.max(jour))+ ' Nombre de fichiers :' + str(nbr_fichier)+ '\n Moyenne : '+ str(mu) +' Ecart type : '+ str(sigma))
plt.xlabel('Taille de la zone de transition en km')
plt.ylabel("Nombre d'occurences")
path = region
path += '_histo_1D_moins.png'
plt.savefig(path)
plt.close()
#plt.show()
def histo_2D_region(SR,alt,dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier / 1.)
latbands = {'tropique' :[-20,20],'pole_nord' :[60,90],'pole_sud' : [-90,-60],'mid_sud' :[-50,-30],'mid_nord' :[30, 50 ]}
dic_hist = dict()
jour = np.zeros((nbr_fichier),dtype = int)
print nbr_fichier
u,v = 100000*(nbr_fichier)/5,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
print np.shape(alt)
alt_array = np.zeros((u,v))
alt_array[:][:] = alt
alt_array = alt_array.ravel()
SRb=[0,0.01,1.2,3,5,7,10,15,20,25,30,40,50,60,80]
altbin = np.arange(-0.240,18.71+0.240,0.480)
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
jour[i] = int(fichier[i][25:27])
print jour[i], fichier[i]
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
print i, '/', nbr_fichier
#idx1 = (lt < latbands[region][1]) & (lt > latbands[region][0])
idx = (lt < 20) & (lt > -20)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['tropique'] = temp
else : dic_hist['tropique'] = dic_hist['tropique'] + temp
idx = (lt < 90) & (lt > 60)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['pole_nord'] = temp
else : dic_hist['pole_nord'] = dic_hist['pole_nord'] + temp
idx = (lt < -60) & (lt > -90)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['pole_sud'] = temp
else : dic_hist['pole_sud'] = dic_hist['pole_sud'] + temp
idx = (lt < -30) & (lt > -50)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['mid_sud'] = temp
else : dic_hist['mid_sud'] = dic_hist['mid_sud'] + temp
idx = (lt < 50) & (lt > 30)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['mid_nord'] = temp
else : dic_hist['mid_nord'] = dic_hist['mid_nord'] + temp
SR_copy = SR.ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['total'] = temp
else : dic_hist['total'] = dic_hist['total'] + temp
for region, tab in dic_hist.items():
x = np.arange(0,np.shape(SRb)[0])
plt.pcolormesh(x, altbin, dic_hist[region].T / np.sum(dic_hist[region]) )
plt.gca().set_xticks(x)
plt.gca().set_xticklabels(SRb)
plt.rcParams["axes.titlesize"] = 10
plt.title(region+' Du ' + fichier[0][17:25]+str(np.min(jour)) +' au ' + fichier[0][17:25] +str(np.max(jour))+ ' Nombre de fichiers :' + str(nbr_fichier))
plt.clim(0,0.001)
#plt.xlim(0,80)
plt.ylim(0,20)
plt.xlabel('SR')
plt.ylabel("Altitude en km")
plt.colorbar()
path = region
path += '_histo_2D.png'
plt.savefig(path)
plt.clf()
def hist_cloud_ext(fichier,date):
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier)
SRb=[0,0.01,1.2,3,5,7,10,15,20,25,30,40,50,60,80]
altbin = np.arange(-0.240,18.71+0.240,0.480)
u,v = np.shape(SR)
alt_array = np.zeros((u,v))
alt_array [:][:] = alt
alt_max = alt_max_orbite_flag(alt, SR)
#idx = ( alt_max != -10)
alt_max = alt_max_orbite_flag(alt, SR)
alt_minSR =alt_min_SR5_orbite(alt_array,SR,alt_max)
idx = (alt_max != -10) & ((alt_array <= alt_max) | (alt_array >= alt_minSR)) & (alt_minSR != -1)
SR = SR[idx].ravel()
alt_array= alt_array[idx].ravel()
fig = plt.figure(figsize=(14, 7), dpi=100)
#plt.hist2d(SR,alt_array,bins= (SRb,altbin))
histSR, x,y = np.histogram2d(SR ,alt_array,bins = (SRb, altbin))
x = np.arange(0,np.shape(SRb)[0])
plt.pcolormesh(x, altbin, histSR.T / np.sum(histSR))
plt.gca().set_xticks(x)
plt.gca().set_xticklabels(SRb)
plt.xlabel('SR')
plt.ylabel('Altitude en km')
plt.ylim(0,19)
plt.title('Histogramme 2D sur une orbite des nuages responsables des extinctions')
plt.clim(0,0.001)
plt.colorbar()
path = './lot/'+ date + '.png'
plt.show()
plt.savefig(path)
plt.clf
plt.close()
def hist_could_ext_dossier(dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier )
for i in range(0,nbr_fichier):
print '== fichier', fichier[i]
print i ,'/', nbr_fichier
calfile = dossier + fichier[i]
date = fichier[i][17:27]
hist_cloud_ext(calfile,date)
def histo_2D_region_eteint(SR,alt,dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier )
latbands = {'tropique' :[-20,20],'pole_nord' :[60,90],'pole_sud' : [-90,-60],'mid_sud' :[-50,-30],'mid_nord' :[30, 50 ]}
dic_hist = dict()
jour = np.zeros((nbr_fichier),dtype = int)
print nbr_fichier
u,v = 100000*(nbr_fichier)/5,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
print np.shape(alt)
alt_array = np.zeros((u,v))
alt_array[:][:] = alt
alt_array = alt_array.ravel()
SRb=[0,0.01,1.2,3,5,7,10,15,20,25,30,40,50,60,80]
altbin = np.arange(-0.240,18.71+0.240,0.480)
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
jour[i] = int(fichier[i][25:27])
print jour[i], fichier[i]
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
alt_max = alt_max_orbite(alt, SR)
print i, '/', nbr_fichier
#idx1 = (lt < latbands[region][1]) & (lt > latbands[region][0])
idx = (lt < 20) & (lt > -20) & (alt_max != -10 )
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['tropique'] = temp
else : dic_hist['tropique'] = dic_hist['tropique'] + temp
idx = (lt < 90) & (lt > 60) & (alt_max != -10 )
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['pole_nord'] = temp
else : dic_hist['pole_nord'] = dic_hist['pole_nord'] + temp
idx = (lt < -60) & (lt > -90) & (alt_max != -10 )
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['pole_sud'] = temp
else : dic_hist['pole_sud'] = dic_hist['pole_sud'] + temp
idx = (lt < -30) & (lt > -50) & (alt_max != -10 )
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['mid_sud'] = temp
else : dic_hist['mid_sud'] = dic_hist['mid_sud'] + temp
idx = (lt < 50) & (lt > 30) & (alt_max != -10 )
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: dic_hist['mid_nord'] = temp
else : dic_hist['mid_nord'] = dic_hist['mid_nord'] + temp
for region, tab in dic_hist.items():
x = np.arange(0,np.shape(SRb)[0])
plt.pcolormesh(x, altbin, dic_hist[region].T / np.sum(dic_hist[region]) )
plt.gca().set_xticks(x)
plt.gca().set_xticklabels(SRb)
plt.rcParams["axes.titlesize"] = 10
plt.title(region + ' Eteints' +' Du ' + fichier[0][17:25]+str(np.min(jour)) +' au ' + fichier[0][17:25] +str(np.max(jour))+ ' Nombre de fichiers :' + str(nbr_fichier))
plt.clim(0,0.001)
#plt.xlim(0,80)
plt.ylim(0,20)
plt.xlabel('SR')
plt.ylabel("Altitude en km")
plt.colorbar()
path = region
path += '_histo_2D_eteint.png'
plt.savefig(path)
plt.clf()
def histo_2D_warmpool_eteint(SR,alt,dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier)
jour = np.zeros((nbr_fichier),dtype = int)
print nbr_fichier
u,v = 100000*(nbr_fichier)/5,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
print np.shape(alt)
alt_array = np.zeros((u,v))
alt_array[:][:] = alt
alt_array = alt_array.ravel()
SRb=[0,0.01,1.2,3,5,7,10,15,20,25,30,40,50,60,80]
altbin = np.arange(-0.240,18.71+0.240,0.480)
for i in range(0,nbr_fichier):
fichier3 = dossier + fichier[i]
jour[i] = int(fichier[i][25:27])
print jour[i], fichier[i]
lg, lt, ti, atb,alt, sol_alt, SR = extraction(fichier3)
alt_max = alt_max_orbite(alt, SR)
print i, '/', nbr_fichier
# idx = (lt < 15) & (lt > -15) & (alt_max != -10 ) & ( lg<180 ) & (lg >70) #warmpool
idx = (lt < 0) & (lt > -15) & (alt_max != -10 ) & ( lg<-75 ) & (lg >-90)
SR_copy = SR[idx,:].ravel()
o = np.max(np.shape(SR_copy))
alt_array_copy = alt_array[0:o]
temp, xx , yy = np.histogram2d(SR_copy, alt_array_copy, bins = (SRb, altbin))
if i == 0: total = temp
else : total = total + temp
x = np.arange(0,np.shape(SRb)[0])
plt.pcolormesh(x, altbin, total.T / np.sum(total) )
plt.gca().set_xticks(x)
plt.gca().set_xticklabels(SRb)
plt.rcParams["axes.titlesize"] = 10
plt.title('Shallowcu ' +' Du ' + fichier[0][17:25]+str(np.min(jour)) +' au ' + fichier[0][17:25] +str(np.max(jour))+ ' Nombre de fichiers :' + str(nbr_fichier))
plt.clim(0,0.001)
#plt.xlim(0,80)
plt.ylim(0,20)
plt.xlabel('SR')
plt.ylabel("Altitude en km")
plt.colorbar()
#path = '_warmpool.png'
path = '_shallow2.png'
plt.savefig(path)
plt.clf()
def distri_alt_max(dossier) :
dic_distri = dict()
latbands = {'tropique' :[-20,20],'pole_nord' :[60,90],'pole_sud' : [-90,-60],'mid_sud' :[-50,-30],'mid_nord' :[30, 50 ]}
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier )
for i in range(0,nbr_fichier):
print i, '/',nbr_fichier
lg, lt, ti, atb,alt, sol_alt, SR = extraction(dossier + fichier[i])
alt_max = alt_max_orbite(alt, SR)
for region in latbands.keys():
idx = (lt < latbands[region][1]) & (lt > latbands[region][0])
temp = alt_max[idx]
if i == 0 : dic_distri[region] = temp
else : dic_distri[region] = np.concatenate([dic_distri[region],temp])
for region, tab in dic_distri.items():
tab = np.ma.masked_where(tab == -10, tab)
bin = np.arange(0,14,0.5)
#plt.hist(hist[0],hist[1],orientation = 'horizontal' )
idx = tab >=0
me = tab[idx]
plt.hist(tab, bins=bin,orientation = 'horizontal')
plt.rcParams["axes.titlesize"] = 8
(mu, sigma) = norm.fit(me)
plt.title(region+' Nombre de fichiers : '+str(nbr_fichier) + '\n Moyenne :' + str(mu) +' km Ecart type : ' + str(sigma) +' km')
plt.xlabel('Nombre d occurences')
plt.ylabel("Epaisseur de la zone de transition (km)")
plt.ylim(0,19)
path = region
path += 'distri_alt_max.png'
plt.savefig(path)
plt.close()
def hist_lonlat(dossier):
fig = plt.figure(figsize=(14, 7), dpi=100)
m = Basemap()
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier / 2)
for i in range(0,nbr_fichier):
lg, lt, ti, atb,alt, sol_alt, SR = extraction(dossier + fichier[i])
print np.shape(lt),np.shape(lg)
alt_max = alt_max_orbite(alt,SR)
idx = (alt_max != -10 )
alt_max = alt_max[idx]
lg = lg[idx]
lt = lt[idx]
lonb = np.arange(-90,90,2)
latb = np.arange(-180,180,2)
if i == 0:
latotal = lt
lgtotal = lg
else:
latotal = np.concatenate([latotal,lt])
lgtotal = np.concatenate([lgtotal,lg])
binlt = np.arange(-90,90,4)
binlg = np.arange(-180,180,4)
print np.shape(bin), np.shape(lt), np.shape(lg)
hist,xx,yy = np.histogram2d( lgtotal, latotal, bins = (binlg,binlt) )
#x = np.arange(0,np.shape(SRb)[0])
m.fillcontinents()
plt.pcolormesh(binlg, binlt, hist.T / np.sum(hist),alpha = 1 )
plt.rcParams["axes.titlesize"] = 8
plt.title(' Distribution de alt max - Nombre de fichiers : '+str(nbr_fichier))
plt.colorbar()
#m.bluemarble()
path = 'lon_lat2.png'
plt.xlim(-180,180)
plt.ylim(-90,90)
plt.savefig(path)
#plt.show()
plt.close()
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=0
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, sol_alt, SR,lt)
# Tableau de flag
#affiche_flag(lg, lt, flag,alt,sol_alt)
#freq_occur(flag,alt)
#n2,x = np.shape(SR)
dossier = '/bdd/CFMIP/CFMIP_OBS_LOCAL/GOCCP/instant_SR_CR_DR/grid_L40/2013/201308/night/'
#freq_occur_region(dossier)
#lonlat(n1, n2,lg, lt, ti, atb)
#afficheSR(lg, lt, ti, SR,alt,sol_alt)
n = 38418
#profil_SR(n,lg,lt,ti,atb,alt,fichier3)
#n = 41750
#profil_SR(n,lg,lt,ti,SR,alt,fichier3)
#freq_occur_region(dossier)
#alti = alt_min_SR5(alt,SR)
#flag = flag_orbite(alt,SR)
#affiche_flag(lg, lt, flag,alt,sol_alt)
#n1,n2 = 38352 , 42466
#histo_2D_orbite(alt,SR,n1,n2,0)#test : 1 = on definit les bornes 0 = sur l'orbite
#alt_min_SR5(alt,SR)
#histo_1D_region2(dossier)
#histo_1D(SR,alt)
#nbr = 2
#debut = 38418
#histo_2D(alt,SR,lt,nbr,debut)
#histo_2D_region(SR,alt,dossier)
#hist_cloud_ext(fichier3,'test_normalisation')
#hist_could_ext_dossier(dossier)
#distri_alt_max(dossier)
#histo_2D_warmpool_eteint(SR,alt,dossier)
#hist_lonlat(dossier)
#draw()
#listevar(fichier3)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment