source: trunk/LMDZ.PLUTO/util/script_figures/proftempNH_obs.py @ 3823

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

############################
zefile="diagfi2015_A.nc"
d1="restart_ref/"

f1=d1+zefile


var="temperature" #variable

nc1=Dataset(f1)
lat=nc1.variables["lat"][:]
lon=nc1.variables["lon"][:]
alt=nc1.variables["altitude"][:]
tim=nc1.variables["time_counter"][:]

print(('Time = ',tim))
print(('Long = ',lon))
print(('Lat = ',lat))

############################

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]))
    return indt1
 
def getvar(filename,var):
    myvar = pp(file=filename,var=var,compute="nothing").getf()
    return myvar

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 main(tini,lt0,ltchoice,p):
    indt=findindextime(tini,lt0,ltchoice,p)
    indlat,indlon=getindex(lat,lon,p[1],p[0])
    myvar=getvar(f1,var)[indt,:,indlat,indlon]
    return myvar
############################

#points: entre -180 et 180 
p1=[-168.75, -18.75] # Entry
p2=[-168.75, -30]    # sud est SP
p2=[-164., -18.75]    # Entry un peu plus sur bord du bassin
p3=[135, 45]    # Burney crater
p4=[180, 45]    # SP
p5=[-124, 2]    # BTD au fond
p6=[-148,16]    # TR
p7=[-5.625,86.25]    # Lowell
p8=[-5.625,-86.25]    # south pole
p7=[129.3,45]    # pente Burney
p8=[166,-18]    # low SP west 
p9=[15,0]    # Exit 

tini=32  # choix du jour dans le diagfi
lt0=0   # local time a t=0 et longitude=0
ltchoice=16.5 # local time choisi
ltchoice2=4 # local time choisi

myvar1=main(tini,lt0,ltchoice,p1)
#myvar2=main(tini,lt0,ltchoice,p2)
#myvar3=main(tini,lt0,ltchoice,p3)
myvar4=main(tini,lt0,ltchoice,p4)
#myvar5=main(tini,lt0,ltchoice,p5)
#myvar6=main(tini,lt0,ltchoice,p6)
#myvar7=main(tini,lt0,ltchoice,p7)
#myvar8=main(tini,lt0,ltchoice,p8)
#myvar9=main(tini,lt0,ltchoice2,p9)


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

font=23

alt=alt/1000.

mpl.plot(myvar1,alt,'r',label='Entry')
#mpl.plot(myvar2,alt,'g',label='Entry-like on edge')
#mpl.plot(myvar3,alt,'k',label='Burney crater')
mpl.plot(myvar4,alt,'k',label='Sputnik Planitia')
#mpl.plot(myvar5,alt,'c',label='BTD')
#mpl.plot(myvar6,alt,'b',label='TR')
#mpl.plot(myvar7,alt,'y',label='Lowell')
#mpl.plot(myvar8,alt,'orange',label='southpole')
#mpl.plot(myvar7,alt,'y',label='Slope Burney')
#mpl.plot(myvar8,alt,'orange',label='South-West SP')
#mpl.plot(myvar9,alt,'b',label='Exit')

mpl.legend(prop={'size':20},loc='upper left')
#mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
mpl.xlabel('Temperature (K)',labelpad=10, fontsize=font)
#mpl.xticks(xticks,fontsize=font)
mpl.xticks(fontsize=font)
#mpl.yticks(yticks,fontsize=font)
mpl.yticks(fontsize=font)
mpl.grid()
#mpl.legend(["Ref","Alt"])

pylab.ylim([-3,2])
pylab.xlim([35,50])

#mpl.savefig('proftempNH_'+str(ltchoice)+'h.eps',dpi=200)
#mpl.savefig('proftempNH_'+str(ltchoice)+'h.png',dpi=200)
mpl.show()