source: trunk/LMDZ.PLUTO/util/script_figures/waves/temp_anomaly.py @ 3831

#! /usr/bin/env python
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
from FV3_utils import *

############################
step_begin = 2000
step_end = 2015

var="temperature" #variable
# tint=["30,37"] #Time must be as written in the input file
tint=None #Time must be as written in the input file
xarea="0,1"
yarea="0,1"

prefix="Xhistins"
suffix="_A.nc"
filename=f"../../{prefix}{step_begin}{suffix}"
nc1=Dataset(filename)

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

myvar=getvar(nc1,var,tint,xarea,yarea)
myvar=myvar[:,:,0,0]
print((shape(myvar)))

# nb of day
nbday=3

# nb of time:
nbtime=size(myvar[:,0])
print(nbtime)

# nb time step /day
nbstep=50

# nb alt
nbalt=size(alt)
print(nbalt)

#meanvar=np.zeros((nbday,nbalt),dtype='f')
meanvar=np.zeros((nbtime-2*nbstep,nbalt),dtype='f')
anovar=np.zeros((nbtime-2*nbstep,nbalt),dtype='f')

# pour chaque jour : calcul moyenne diurne
for i in range(nbtime-2*nbstep):
    #i=i+nbstep/2
    meanvar[i,:]=np.mean(myvar[nbstep//2+i:nbstep+nbstep//2+i,:],axis=0)
#for i in range(nbday):
#    meanvar[i,:]=np.mean(myvar[0:8,:],axis=0)

# pour chaque time : calcul anomaly le dernier time est pour autre jour
for i in range(nbtime-2*nbstep):
    #index=int(i/nbstep)
    anovar[i,:]=myvar[nbstep+i,:]-meanvar[i,:]

print((meanvar[:,:]))
print((myvar[0,:]))
# print((myvar[8,:]))


mpl.figure(figsize=(20, 10))

font=26

#pal=get_cmap(name="RdYlBu_r")
# pal=get_cmap(name="Spectral_r")
lev=np.linspace(-0.1,0.1,10)

#xticks=[-90,-60,-30,0,30,60,90]
#yticks=np.linspace(0,240,9)
tim=np.arange(nbtime-2*nbstep)/floor(nbstep)
#tim=np.arange(nbtime)/nbstep

print((shape(tim), shape(alt),shape(anovar)))

CF=mpl.contourf(tim,alt,np.transpose(anovar),lev,cmap="coolwarm",extend='both')
cbar=mpl.colorbar(CF,shrink=1, format="%1.2f")
cbar.ax.set_title("[K]",y=1.04,fontsize=font)
for t in cbar.ax.get_yticklabels():
      t.set_fontsize(font)

#vect=lev
#CS=mpl.contour(lat,alt,myvar,vect,colors='k',linewidths=0.5)
#### inline=1 : values over the line
#mpl.clabel(CS, inline=1, fontsize=20, fmt='%1.0f',inline_spacing=1)

#mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
mpl.xlabel('Time (Pluto days)',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,150000])
# pylab.ylim([0,np.max(alt)])

mpl.savefig('tempanom.eps',dpi=200)
mpl.savefig('tempanom.png',dpi=200)
mpl.show()