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

Last change on this file was 3833, checked in by afalco, 36 hours ago

Pluto: updated plots scripts.
Fixed some issues with reading XIOS, etc.
Included display_netcdf.py tool from Mars PCM.
AF

  • Property svn:executable set to *
File size: 3.4 KB
Line 
1#! /usr/bin/env python
2from    netCDF4               import    Dataset
3from    numpy                 import    *
4import  numpy                 as        np
5import  matplotlib.pyplot     as        mpl
6from matplotlib.cm import get_cmap
7import pylab
8from matplotlib import ticker
9import matplotlib.colors as colors
10import datetime
11from mpl_toolkits.basemap import Basemap, shiftgrid
12from FV3_utils import *
13
14############################
15var="temperature" #variable
16
17nc1=Dataset(name+"_A.nc")
18lat=getvar(nc1,"latitude")
19lon=getvar(nc1,"longitude")
20alt=getvar(nc1,"altitude")
21tim=getvar(nc1,"Time")
22
23print(('Time = ',tim))
24print(('Long = ',lon))
25print(('Lat = ',lat))
26
27############################
28
29def findindextime(tini,lt0,lt1,p):
30    lt180=(lt0+12)%24   # a t=0 et longitude=180
31    lt0p1=lt180*((p[0]+360)%360)/180   # local time a t=0 a la longitude p1
32    diff=lt1-lt0p1                 # diff du local time a p1 a t=0 avec celui recherche
33    indp1=tini+diff/24.                 # on adapte lindice pour tomber sur le bon moment
34
35    # on cherche dans Time l'indice le plus proche :
36    indt1=np.where(abs(tim[:]-indp1)==min(abs(tim[:]-indp1)))[0][0]
37    print(('Point =',p,' Time=',tim[indt1]))
38    return indt1
39
40def getindex(lat,lon,mylat,mylon):
41    indlat=np.where(abs(lat[:]-mylat)==min(abs(lat[:]-mylat)))[0][0]
42    indlon=np.where(abs(lon[:]-mylon)==min(abs(lon[:]-mylon)))[0][0]
43    print((lon[indlon],lat[indlat]))
44    return indlat,indlon
45
46def main(tini,lt0,ltchoice,p):
47    indt=findindextime(tini,lt0,ltchoice,p)
48    indlat,indlon=getindex(lat,lon,p[1],p[0])
49    myvar=getvar(nc1,var)[indt,:,indlat,indlon]
50    return myvar
51############################
52
53#points: entre -180 et 180
54p1=[-168.75, -18.75] # Entry
55p2=[-168.75, -30]    # sud est SP
56p2=[-164., -18.75]    # Entry un peu plus sur bord du bassin
57p3=[135, 45]    # Burney crater
58p4=[180, 45]    # SP
59p5=[-124, 2]    # BTD au fond
60p6=[-148,16]    # TR
61p7=[-5.625,86.25]    # Lowell
62p8=[-5.625,-86.25]    # south pole
63p7=[129.3,45]    # pente Burney
64p8=[166,-18]    # low SP west
65
66tini=32  # choix du jour dans le diagfi
67lt0=0   # local time a t=0 et longitude=0
68ltchoice=16.5 # local time choisi
69
70myvar1=main(tini,lt0,ltchoice,p1)
71myvar2=main(tini,lt0,ltchoice,p2)
72myvar3=main(tini,lt0,ltchoice,p3)
73myvar4=main(tini,lt0,ltchoice,p4)
74myvar5=main(tini,lt0,ltchoice,p5)
75myvar6=main(tini,lt0,ltchoice,p6)
76myvar7=main(tini,lt0,ltchoice,p7)
77myvar8=main(tini,lt0,ltchoice,p8)
78
79
80mpl.figure(figsize=(10, 10))
81
82font=23
83
84mpl.plot(myvar1,alt,'r',label='Entry')
85mpl.plot(myvar2,alt,'g',label='Entry-like on edge')
86mpl.plot(myvar3,alt,'k',label='Burney crater')
87mpl.plot(myvar4,alt,'m',label='SP')
88mpl.plot(myvar5,alt,'c',label='BTD')
89mpl.plot(myvar6,alt,'b',label='TR')
90#mpl.plot(myvar7,alt,'y',label='Lowell')
91#mpl.plot(myvar8,alt,'orange',label='southpole')
92mpl.plot(myvar7,alt,'y',label='Slope Burney')
93mpl.plot(myvar8,alt,'orange',label='South-West SP')
94
95mpl.legend(prop={'size':20},loc='upper left')
96#mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
97mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
98mpl.xlabel('Temperature (K)',labelpad=10, fontsize=font)
99#mpl.xticks(xticks,fontsize=font)
100mpl.xticks(fontsize=font)
101#mpl.yticks(yticks,fontsize=font)
102mpl.yticks(fontsize=font)
103mpl.grid()
104#mpl.legend(["Ref","Alt"])
105
106pylab.ylim([-4,6])
107pylab.xlim([30,60])
108
109mpl.savefig('proftempNH_'+str(ltchoice)+'h.eps',dpi=200)
110mpl.savefig('proftempNH_'+str(ltchoice)+'h.png',dpi=200)
111mpl.show()
112
113
Note: See TracBrowser for help on using the repository browser.