Skip to content

Instantly share code, notes, and snippets.

@stageipsl

stageipsl/hist1Dregion

Created May 21, 2014 12:35
Show Gist options
  • Save stageipsl/91e67c6ec7d431d51b36 to your computer and use it in GitHub Desktop.
Save stageipsl/91e67c6ec7d431d51b36 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
hist1Dregion
#!/usr/bin/env python
import numpy as np
import numpy.matlib
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'][:,:]
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):
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 profil(n,lg,lt,ti,atb,alt):
plt.semilogx(atb[n,:], alt)
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, sol_alt, SR,lt):
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
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_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=(9, 6), dpi=100)
plt.plot(lt, sol_alt,'w', alpha=0.4)
plt.pcolormesh(lt,alt,atb.T)
plt.ylim(-1, 20)
plt.xlim(lt[0], lt[-1])
plt.clim(0, 4)
plt.xlabel('Latitude')
plt.ylabel('Altitude')
plt.colorbar()
plt.show()
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) / 10
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(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 histo_1D_orbite(SR,alt,dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = int(nbr_fichier / 5 )
#u,v = 56295,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
#alt_array = np.zeros((u,v))
#alt_array[:][:] = alt
alt_array = np.zeros_like(SR)
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)
plt.hist(total, bins=np.arange(0,14,0.5))
plt.title('Histogramme d epaisseur de la zone de transition Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier))
plt.show()
def histo_2D(alt, SR):
alt_array = np.zeros_like(SR)
alt_array[:][:] = alt
plt.hist2d(SR , alt_array)
plt.colorbar()
plt.show()
def histo_1D_region(SR,alt,dossier):
nbr_fichier,fichier = lecture_dossier(dossier)
nbr_fichier = 5 #int(nbr_fichier / 5 )
u,v = 56295,40 # passe manuellement car certains SR ont des tailles inferieur a 56295
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()
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]
for region in latbands.keys():
#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 < latbands[region][1]) & (lt > latbands[region][0])
temp = alt_SR5[idx]-alt_max[idx,0]
if i == 0: dic_hist[region] = temp
else : dic_hist[region] = np.concatenate([dic_hist[region], temp])
for region in latbands.keys():
plt.hist(dic_hist[region], bins=np.arange(0,14,0.5))
plt.title('Histogramme d epaisseur de la zone de transition Du ' + fichier[0][17:27] +' au ' + fichier[nbr_fichier][17:27] + ' Nombre de fichiers :' + str(nbr_fichier))
path = region
path += 'test.png'
plt.savefig(path)
plt.close()
#plt.show()
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)
#profil_SR(n,lg,lt,ti,atb,alt,fichier)
#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)
#histo_2D(alt,SR)
#alt_min_SR5(alt,SR)
histo_1D_region(SR,alt,dossier)
#histo_1D(SR,alt)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment