#! /usr/bin/env python from ppclass import pp from netCDF4 import Dataset from numpy import * import numpy as np import matplotlib.pyplot as mpl from matplotlib.cm import get_cmap import pylab from matplotlib import ticker import matplotlib.colors as colors import datetime from mpl_toolkits.basemap import Basemap, shiftgrid ############################ fi="diagfi2015_S.nc" d1="../restart_2/" f1=d1+fi var="temperature" #variable p1=[-148,16] nc1=Dataset(f1) lat=nc1.variables["lat"][:] lon=nc1.variables["lon"][:] alt=nc1.variables["altitude"][:] tim=nc1.variables["time_counter"][:] def findindextime(tini,lt0,lt1,p): lt180=(lt0+12)%24 # a t=0 et longitude=180 lt0p1=lt180*((p[0]+360)%360)/180 # local time a t=0 a la longitude p1 diff=lt1-lt0p1 # diff du local time a p1 a t=0 avec celui recherche indp1=tini+diff/24. # on adapte lindice pour tomber sur le bon moment # on cherche dans Time l'indice le plus proche : indt1=np.where(abs(tim[:]-indp1)==min(abs(tim[:]-indp1)))[0][0] print(('Point =',p,' Time=',tim[indt1],' diff = ',diff)) return indt1 def getindex(lat,lon,mylat,mylon): indlat=np.where(abs(lat[:]-mylat)==min(abs(lat[:]-mylat)))[0][0] indlon=np.where(abs(lon[:]-mylon)==min(abs(lon[:]-mylon)))[0][0] print((lon[indlon],lat[indlat])) return indlat,indlon def getvar(filename,var): myvar = pp(file=filename,var=var,compute="nothing").getf() print(('shape myvar = ',shape(myvar))) return myvar tini=32 # choix du jour dans le diagfi lt0=0 # local time a t=0 et longitude=0 ltchoice0=0 # local time initial ltchoice1=25 # local time final indt0=findindextime(tini,lt0,ltchoice0,p1) indt1=findindextime(tini,lt0,ltchoice1,p1) print((indt0,indt1)) indlat,indlon=getindex(lat,lon,p1[1],p1[0]) mpl.figure(figsize=(18, 10)) myvar=getvar(f1,var)[indt0:indt1+1,:,indlat,indlon] tim=tim[indt0:indt1+1] print(("time is : ",tim)) font=26 axtim=(tim-tim[0])*24 print(("axtim=",axtim)) print(('on prend les premiers indice, shape (tmps, alt, var) =',shape(axtim), shape(alt), shape(myvar))) pal=get_cmap(name="Spectral_r") lev=np.linspace(34,54,19) xticks=[0,2,4,6,8,10,12,14,16,18,20,22,24] alt=alt/1000. CF=mpl.contourf(axtim,alt,np.transpose(myvar),lev,cmap=pal,extend='both') cbar=mpl.colorbar(CF,shrink=1, format="%1.0f") cbar.ax.set_title("Temp [K]",y=1.04,fontsize=font) for t in cbar.ax.get_yticklabels(): t.set_fontsize(font) vect=lev CS=mpl.contour(axtim,alt,np.transpose(myvar),vect,colors='k',linewidths=0.5) #inline=1 : values over the line #mpl.clabel(CS, inline=2, fontsize=10, fmt='%1.1e') lab=mpl.clabel(CS, inline=1, fontsize=15, fmt='%1.0f',inline_spacing=1) for l in lab: l.set_rotation(0) #mpl.title('Latitude ='+str(tintstr[i]),fontsize=font) mpl.xlabel('Local Time (h)',labelpad=10,fontsize=font) mpl.ylabel('Altitude (km)',labelpad=10, fontsize=font) mpl.xticks(xticks,fontsize=font) #mpl.xticks(fontsize=font) #mpl.yticks(yticks,fontsize=font) mpl.yticks(fontsize=font) pylab.ylim([0,4]) pylab.xlim([0,24]) mpl.savefig('temploctime_2_'+str(p1[0])+'_'+str(p1[1])+'.eps',dpi=200) mpl.savefig('temploctime_2_'+str(p1[0])+'_'+str(p1[1])+'.png',dpi=200) #mpl.show()