source: trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom_ref.py

#! /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_ref/"

f1=d1+fi

var="temperature" #variable

p1=[-148,16]

nc1=Dataset(f1)
lat=getvar(nc1,"latitude")
lon=getvar(nc1,"longitude")
alt=getvar(nc1,"altitude")
tim=getvar(nc1,"Time")

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_ref_'+str(p1[0])+'_'+str(p1[1])+'.eps',dpi=200)
mpl.savefig('temploctime_ref_'+str(p1[0])+'_'+str(p1[1])+'.png',dpi=200)
#mpl.show()