Changeset 3833


Ignore:
Timestamp:
Jul 7, 2025, 1:39:13 PM (41 hours ago)
Author:
afalco
Message:

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

Location:
trunk/LMDZ.PLUTO/util/script_figures
Files:
8 added
1 deleted
21 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/util/script_figures/FV3_utils.py

    r3823 r3833  
    33pi=math.pi
    44import matplotlib
    5 matplotlib.use('TKAgg')
    6 
    7 name="../Xhistins2015_short"
     5matplotlib.use('Agg')
     6from input import *
     7from display_netcdf import *
    88
    99def fmtsci(x, pos):
    … …  
    12011201
    12021202def printvar(infile):
    1203     # Get all variable names for a netCDF file ---
     1203    '''Get all variable names for a netCDF file '''
    12041204    print("Variables:")
    12051205    variableNames = list(infile.variables.keys());
    … …  
    12081208
    12091209def getind(myloc,field):
    1210     # get a specific index in the lat, lon, time or pfull 1D field
     1210    '''get a specific index in the lat, lon, time or pfull 1D field'''
    12111211    res=[]
    12121212    for loc in np.atleast_1d(myloc):
    … …  
    12201220# getinds=np.vectorize(getind, excluded="field")
    12211221
    1222 def getvar(nc1,var,times=None,tim=None,longitudes=None,lon=None,t_mean=False,l_mean=False):
     1222def getvar(nc1,var,times=None,longitudes=None,latitudes=None,altitudes=None,t_mean=False,l_mean=False):
     1223    """Get variables from netcdf file and slice it up!
     1224
     1225    Args:
     1226        nc1 : netcdf identifier
     1227        var (str): variable name
     1228        times (list, optional): times to select. Defaults to None.
     1229        longitudes (list, optional): longitudes to select. Defaults to None.
     1230        latitudes (list, optional): latitudes to select. Defaults to None.
     1231        t_mean (bool, optional): Apply mean on time dimension. Defaults to False.
     1232        l_mean (bool, optional): Apply mean on longitude dimension. Defaults to False.
     1233
     1234    Returns:
     1235        array: variable data array
     1236    """
     1237   
    12231238    if var == "latitude" and var not in nc1.variables: var="lat"
    12241239    if var == "longitude" and var not in nc1.variables: var="lon"
    … …  
    12271242    if var == "phisinit" and var not in nc1.variables: var="phisfi"
    12281243
    1229     myvar = nc1.variables[var][:]
    1230     try:
    1231         t=getind(times,tim)
    1232         if len(t)==2:
    1233             myvar = myvar[t[0]:t[1]+1]
    1234         else:
    1235             myvar = myvar[t]
    1236     except:
    1237         # print("no time selected")
    1238         pass
    1239     try:
    1240         l=getind(longitudes,lon)
    1241         if len(l)==2:
    1242             myvar = myvar[...,l[0]:l[1]+1]
    1243         else:
    1244             myvar = myvar[...,l]
    1245     except:
    1246         # print("no longitude selected")
    1247         pass
    1248     if t_mean: myvar=np.mean(myvar,axis=0) # temporal mean
    1249     if l_mean: myvar=np.mean(myvar,axis=-1) # longitudinal mean
     1244    var_data = nc1.variables[var]
     1245    dims = var_data.dimensions
     1246    myvar = var_data[:]
     1247
     1248    if times or longitudes or latitudes:
     1249        slicer = [slice(None)] * len(dims)
     1250
     1251        time_it = (times,TIME_DIMS)
     1252        alt_it = (altitudes,ALT_DIMS)
     1253        lon_it = (longitudes,LON_DIMS)
     1254        lat_it = (latitudes,LAT_DIMS)
     1255
     1256        # iterate over dimensions to pick indices to select
     1257        for values, dim in (time_it, lon_it, lat_it, alt_it):
     1258            try:
     1259                dim_idx = find_dim_index(dims, dim)
     1260                tmp_var = find_coord_var(nc1, dim)
     1261                file_values = nc1.variables[tmp_var][:]
     1262                i = getind(values, file_values)
     1263                print(dim, i)
     1264                if len(i)==2: # range of indices
     1265                    slicer[dim_idx] = slice(min(i),max(i)+1)
     1266                else: # specific index
     1267                    slicer[dim_idx] = i
     1268            except:
     1269                continue
     1270       
     1271        # slice the variable according to the indices
     1272        myvar = myvar[tuple(slicer)]
     1273       
     1274    if l_mean: # longitudinal mean
     1275        myvar=np.mean(myvar,axis=find_dim_index(dims, LON_DIMS))
     1276    if t_mean: # temporal mean
     1277        myvar=np.mean(myvar,axis=find_dim_index(dims, TIME_DIMS))
    12501278    print(np.shape(myvar))
    12511279    try:

  • trunk/LMDZ.PLUTO/util/script_figures/ch4cloudsection.py

    r3823 r3833  
    2323nc1=Dataset(filename1)
    2424
    25 lat=nc1.variables["lat"][:]
    26 lon=nc1.variables["lon"][:]
    27 alt=nc1.variables["altitude"][:]
     25lat=getvar(nc1,"latitude")
     26lon=getvar(nc1,"longitude")
     27alt=getvar(nc1,"altitude")
    2828print(('alt=',alt))
    2929############################

  • trunk/LMDZ.PLUTO/util/script_figures/ch4vmrsection.py

    r3823 r3833  
    88import matplotlib.colors as mcolors
    99from FV3_utils import *
     10from input import *
    1011
    1112############################
    12 folder="../"
    13 filename1=folder+"diagfi_mean_A.nc"
    1413var="ch4_gas" #variable
    1514tint=[30,32] #Time must be as written in the input file
    1615xarea="-180,179"
    1716
    18 nc1=Dataset(filename1)
     17nc1=Dataset(name+"_A.nc")
    1918
    20 lat=nc1.variables["lat"][:]
    21 lon=nc1.variables["lon"][:]
    22 alt=nc1.variables["altitude"][:]
    23 tim=nc1.variables["time_counter"][:]
     19lat=getvar(nc1,"latitude")
     20lon=getvar(nc1,"longitude")
     21alt=getvar(nc1,"altitude")
     22tim=getvar(nc1,"Time")
    2423############################
    2524
    26 # def getvar(filename,var,tint,xarea):
    27 #     myvar = pp(file=filename,var=var,t=tint,x=xarea,compute="mean").getf()
    28 #     print((shape(myvar)))
    29     # return myvar
    30 
    31 myvar=getvar(nc1,var,tint,tim)
     25myvar=getvar(nc1,var,tint, l_mean=True, t_mean=True)
    3226myvar=myvar*28/16.*100.
    3327
    … …  
    4236xticks=[-90,-60,-30,0,30,60,90]
    4337#yticks=np.linspace(0,240,9)
    44 alt=alt/1000.
    4538
    4639mymin=0.1

  • trunk/LMDZ.PLUTO/util/script_figures/map_alb_ini.py

    r3823 r3833  
    11#! /usr/bin/env python
    2 from ppclass import pp
    32from    netCDF4               import    Dataset
    43from    numpy                 import    *
    … …  
    109import matplotlib.colors as mcolors
    1110from mpl_toolkits.basemap import Basemap
     11from FV3_utils import *
    1212
    1313############################
    14 filename="diagfi2015.nc"
     14filename=name+".nc"
    1515nb_dataset=1  #change vect as well
    16 tint=["0,0.001"] #Time must be as written in the input file
     16tint=[0,0.001] #Time must be as written in the input file
    1717#zint=["0"] #alt in km
    1818#xarea="-180,179"
    1919#yarea="-90,90"
    2020#step=5  #step between each bin
    21 var="albedo" #variable
     21var="ALB" #variable
    2222var2="phisinit"
    2323############################
    2424
    2525nc=Dataset(filename)
    26 lat=nc.variables["lat"][:]
    27 lon=nc.variables["lon"][:]
     26lat=getvar(nc,"latitude")
     27lon=getvar(nc,"longitude")
     28tim=getvar(nc,"Time")
    2829
    29 myvar = pp(file=filename,var=var,t=tint,compute="mean").getf()  # get data to be changed according to selected variable
    30 myvar2 = pp(file=filename,var=var2,t=tint,compute="mean").getf()  # get data to be changed according to selected variable
     30myvar = getvar(nc,var,tint,t_mean=True)
     31myvar2 = getvar(nc,var2)
    3132print(('shape var1 =',shape(myvar)))
    3233print(('shape lat = ',shape(lat)))
    … …  
    9091           myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)]
    9192
    92 lon=lon+180.
     93if lon[0]<0:
     94    lon=lon+180.
    9395
    9496

  • trunk/LMDZ.PLUTO/util/script_figures/mapch4cloud.py

    r3823 r3833  
    2626nc2=Dataset(filename2)
    2727
    28 lat=nc1.variables["lat"][:]
    29 lon=nc1.variables["lon"][:]
     28lat=getvar(nc1,"latitude")
     29lon=getvar(nc1,"longitude")
    3030
    3131# altitdue file 2

  • trunk/LMDZ.PLUTO/util/script_figures/mapmeanwinds.py

    r3823 r3833  
    88from matplotlib import ticker
    99import matplotlib.colors as colors
    10 import datetime
    1110from mpl_toolkits.basemap import Basemap, shiftgrid
    1211from FV3_utils import *
     12from input import *
    1313
    14 ############################
    1514fa='sans-serif'
    1615hfont = {'fontname':'Arial'}
    … …  
    2524lvls=np.linspace(37,57,21)
    2625
     26altitude=1 # in km
     27
    2728### Data
    28 name2=name+'.nc'
    29 # name=name+'_A.nc'
    30 name=name+'.nc'
    31 print("Plot "+name)
    32 nc1=Dataset(name)
    33 nc2=Dataset(name2)
     29filename1=name+'_A.nc'
     30# name=name+'.nc'
     31filename2=name+'.nc'
     32print("Plot "+filename1)
     33nc1=Dataset(filename1)
     34nc2=Dataset(filename2)
    3435ts=nc2.variables["temperature"][:,0,:,:]
    3536# ts=nc2.variables["tsurf"][:,:,:]
    3637u=nc1.variables["u"][:,:,:,:]
    3738v=nc1.variables["v"][:,:,:,:]
    38 lat=nc1.variables["lat"][:]
    39 alt=nc1.variables["altitude"][:]
    40 lon=nc1.variables["lon"][:]
     39lat=getvar(nc1,"latitude")
     40alt=getvar(nc1,"altitude")
     41lon=getvar(nc1,"longitude")
    4142ps=nc2.variables["ps"][:,:]
    4243# ps=nc2.variables["phisinit"][:,:]/0.6169/1000.  # altitude km
    4344
    44 numalt=getind(1,alt)
    45 # numalt=getind(1000,alt)
     45numalt=getind(altitude,alt)
    4646print('numalt =',numalt,'altitude=',alt[numalt])
    4747u=u[:,numalt,:,:]
    … …  
    5252ps=np.mean(ps,axis=0)
    5353
    54 ts=switchlon(ts)
    55 u=switchlon(u)
    56 v=switchlon(v)
    57 # ps=switchlon(ps)
    58 # topo=switchlon(topo)
    59 # lon=lon+180.
     54if lon[0]<0:
     55    ts=switchlon(ts)
     56    u=switchlon(u)
     57    v=switchlon(v)
     58    ps=switchlon(ps)
     59    # topo=switchlon(topo)
     60    lon=lon+180.
    6061
    6162### Figure
    … …  
    9697mpl.show()
    9798
    98 #######################
    99 

  • trunk/LMDZ.PLUTO/util/script_figures/maptemp.py

    r3823 r3833  
    11#! /usr/bin/env python
     2import os
    23from    netCDF4               import    Dataset
    34from    numpy                 import    *
    … …  
    1112from mpl_toolkits.basemap import Basemap, shiftgrid
    1213from matplotlib.cm import get_cmap
    13 from FV3_utils import * # import name
     14from FV3_utils import *
     15from input import * # import name
     16print("Running "+os.path.basename(__file__))
    1417
    1518############################
    … …  
    2629
    2730### Data
    28 # name='../diagfi2015.nc' # read from FV3_util
    29 print(name)
    3031try:
     32    print("Plotting "+name+"_A.nc")
    3133    nc1=Dataset(name+"_A.nc")
    3234except:
     35    print("Plotting "+name+".nc")
    3336    nc1=Dataset(name+".nc")
    34 alt=nc1.variables["altitude"][:]
    35 lat=nc1.variables["lat"][:]
    36 lon=nc1.variables["lon"][:]
     37alt=getvar(nc1,"altitude")
     38lat=getvar(nc1,"latitude")
     39lon=getvar(nc1,"longitude")
    3740# temp=switchlon(temp)
    38 lon=lon+180.
    3941
    40 def plot_alt(altitude = 1000):
     42if lon[0]<0:
     43    lon=lon+180.
     44
     45def plot_alt(altitude = 1):
    4146    temp=nc1.variables["temperature"][:,:,:,:]
    4247    numalt=getind(altitude,alt)
    … …  
    6368
    6469    mpl.grid()
    65     mpl.title(f"Temperatures @ z={altitude/1000}km",fontsize=font)
     70    mpl.title(f"Temperatures @ z={altitude}km",fontsize=font)
    6671    mpl.ylabel(r'Latitude',labelpad=10,fontsize=font)
    6772    mpl.xlabel('Longitude',labelpad=10, fontsize=font)
    … …  
    7277    mpl.yticks(yticks,fontsize=font)
    7378    mpl.xticks(xticks,fontsize=font)
    74     mpl.savefig(f"maptemp{altitude}",bbox_inches='tight',dpi=70)
     79    output=f"maptemp{altitude}"
     80    mpl.savefig(output,bbox_inches='tight',dpi=70)
     81    print(f"Saved {output}")
    7582    #mpl.show()
    7683
    7784
    78 plot_alt(1000)
    79 plot_alt(5000)
    80 plot_alt(20000)
    81 plot_alt(50000)
    82 plot_alt(100000)
     85plot_alt(1)
     86plot_alt(5)
     87plot_alt(20)
     88plot_alt(50)
     89plot_alt(100)

  • trunk/LMDZ.PLUTO/util/script_figures/meanmeridwind.py

    r3823 r3833  
    99import colorsys
    1010from FV3_utils import *
     11from input import *
    1112
    1213############################
    … …  
    6364rvb1=diverge_map(h,l)
    6465
    65 myvar=getvar(nc1,var,tint,tim,xarea,lon,t_mean=True,l_mean=True)
     66myvar=getvar(nc1,var,tint,xarea,lon,t_mean=True,l_mean=True)
    6667
    6768mpl.figure(figsize=(20, 10))
    … …  
    7677xticks=[-90,-60,-30,0,30,60,90]
    7778#yticks=np.linspace(0,240,9)
    78 alt=alt/1000.
    7979
    8080CF=mpl.contourf(lat,alt,myvar,lev,cmap=pal,extend='both')
    … …  
    8989mpl.clabel(CS, inline=1, fontsize=20, fmt='%1.1f',inline_spacing=1)
    9090
    91 #mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
     91mpl.title('Meridional wind (m/s)',fontsize=font)
    9292mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
    9393mpl.xlabel('Latitude (deg)',labelpad=10, fontsize=font)

  • trunk/LMDZ.PLUTO/util/script_figures/meanzonalwind.py

    r3823 r3833  
    99import colorsys
    1010from FV3_utils import *
     11from input import *
    1112
    1213############################
    … …  
    6465rvb1=diverge_map(h,l)
    6566
    66 myvar=getvar(nc1,var,tint,tim,xarea,lon,t_mean=True,l_mean=True)
     67myvar=getvar(nc1,var,tint,xarea,lon,t_mean=True,l_mean=True)
    6768
    6869mpl.figure(figsize=(20, 10))
    … …  
    7778xticks=[-90,-60,-30,0,30,60,90]
    7879#yticks=np.linspace(0,240,9)
    79 alt=alt/1000.
    8080
    8181CF=mpl.contourf(lat,alt,myvar,lev,cmap=pal,extend='both')
    … …  
    9090mpl.clabel(CS, inline=1, fontsize=20, fmt='%1.0f',inline_spacing=1)
    9191
    92 #mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
     92mpl.title('Zonal wind (m/s)',fontsize=font)
    9393mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
    9494mpl.xlabel('Latitude (deg)',labelpad=10, fontsize=font)

  • trunk/LMDZ.PLUTO/util/script_figures/meanzonalwind_zoom.py

    r3823 r3833  
    1818nc1=Dataset(filename1)
    1919
    20 lat=nc1.variables["lat"][:]
    21 lon=nc1.variables["lon"][:]
    22 alt=nc1.variables["altitude"][:]
     20lat=getvar(nc1,"latitude")
     21lon=getvar(nc1,"longitude")
     22alt=getvar(nc1,"altitude")
    2323############################
    2424

  • trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/FV3_utils.py

    r3832 r3833  
    33pi=math.pi
    44import matplotlib
    5 matplotlib.use('TKAgg')
    6 
    7 name="../Xhistins2015_short"
     5matplotlib.use('Agg')
    86
    97def fmtsci(x, pos):

  • trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/vmr_ch4.py

    r3832 r3833  
    1010import datetime
    1111from mpl_toolkits.basemap import Basemap, shiftgrid
    12 from FV3_utils import *
     12from FV3_utils import * # import name
     13from input import * # import 'name'
    1314
    1415############################
    15 # basename="../../Xhistins2072"
    16 filename1="../"+name+"_A.nc"
     16filename1="../"+name+"_A.nc" # define name in input.py
    1717filename2="../"+name+".nc"
    18 filename3="../../phisinit.nc"
    19 var="tsurf" #variable
    20 phisinit="phisinit" #variable
     18var="vmr_ch4"
     19var2="phisinit" #variable
    2120vari="u"
    2221varj="v"
    23 tint=[30,32] #Time must be as written in the input file
     22font=26
     23tint=[30,31] #Time must be as written in the input file
    2424#tintstr=["03:00","09:00","15:00","21:00"] #Time must be as written in the input file
    2525
    26 font=26
     26# altitude in km
     27altitude=1
    2728
    2829nc1=Dataset(filename1)
    2930nc2=Dataset(filename2)
    30 nc3=Dataset(filename3)
    3131
    3232lat=getvar(nc1,"latitude")
    … …  
    3535tim=getvar(nc1,"Time")
    3636############################
     37
     38def swinglon(myvar):
     39   # changer les longitudes pour mettre TR au centre
     40   vec=shape(myvar)
     41   myvarbis=np.zeros(vec,dtype='f')
     42# i lat : pas de changement
     43# j lon :
     44   for i in range(vec[0]):
     45       for j in range(vec[1]):
     46           if j < int(vec[1]/2.) :
     47             myvarbis[i,j]=myvar[i,j+int(vec[1]/2.)]
     48             # myvar2bis[i,j]=myvar2[i,j+int(vec[1]/2)]
     49           else:
     50             myvarbis[i,j]=myvar[i,j-int(vec[1]/2.)]
     51             # myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)]
     52   return myvarbis
    3753
    3854def getwinds(lon,lat,vecx,vecy):
    … …  
    4460          color='black'    # arrow color
    4561          pivot='mid'      # arrow around middle of box. Alternative : tip
    46           scale=100     # scale arrow : plus grand = fleche plus petite
     62          scale=200     # scale arrow : plus grand = fleche plus petite
    4763          width=0.002      # width arrow
    48           linewidths=0.5   # epaisseur contour arrow
     64          linewidths=0.1   # epaisseur contour arrow
    4965          edgecolors='k'   # couleur contour arrow
    5066
    … …  
    86102def getfigvar(i):
    87103    pal=get_cmap(name="jet")
    88     lev=np.linspace(37,51,15)
    89     # newlon=lon+180
    90     newlon=lon
     104    lev=np.linspace(0.2,0.8,50)
     105    newlon=lon+180
    91106    CF=mpl.contourf(newlon, lat, myvarbis,lev,cmap=pal,extend='both')
    92107    yticks=[-90,-60,-30,0,30,60,90]
    93108    xticks=[0,60,120,180,240,300,360]
    94     cbar=mpl.colorbar(CF,shrink=1, format="%1.0f")
    95     cbar.ax.set_title("Tsurf [K]",y=1.04,fontsize=font)
     109    cbar=mpl.colorbar(CF,shrink=1, format="%1.2f")
     110    cbar.ax.set_title("VMR CH4 [%]",y=1.04,fontsize=font)
    96111
    97112    for t in cbar.ax.get_yticklabels():
    98113        t.set_fontsize(font)
    99114
    100     c=mpl.contour(newlon, lat, phisinit2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
     115    c=mpl.contour(newlon, lat, myvar2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
    101116    mpl.clabel(c, fmt='%2.1f',inline=1, colors='k', fontsize=23,inline_spacing=1)
    102117    #mpl.title('Local Time at Sputnik Planum='+str(i*3)+'H00',fontsize=font)
    … …  
    105120    mpl.xticks(xticks,fontsize=font)
    106121    mpl.yticks(yticks,fontsize=font)
    107     getwinds(newlon,lat,u,v)
     122    #getwinds(newlon,lat,u,v)
    108123
    109 # def getnumalt(choicealt,alt):
    110 #     numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
    111 #     numalt=numalt[0][0]
    112 #     return numalt
     124def getnumalt(choicealt,alt):
     125    numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
     126    numalt=numalt[0][0]
     127    return numalt
    113128
    114129#######################
    115 # numalt=getnumalt(30,alt)
    116 numalt=np.searchsorted(alt[...],30)
     130numalt=getnumalt(altitude,alt)
    117131print(('numalt =',numalt,'altitude=',alt[numalt]))
    118132uini=getvar(nc1,vari,tint,tim)[:,numalt]
    119133vini=getvar(nc1,varj,tint,tim)[:,numalt]
    120134myvar=getvar(nc2,var,tint,tim)
    121 phisinit2=getvar(nc3,phisinit) # phisinitmyvar2=phisinit2/0.6169/1000.  # altitude km
     135myvar2=getvar(nc2,var2) # phisinit
     136myvar2=myvar2/0.6169/1000.  # altitude km
    122137nbfig=uini.shape[0]
    123138print(("nbfig=",nbfig))
    124 phisinit2bis=np.copy(phisinit2)
     139myvar2bis=swinglon(myvar2)
    125140
    126141for i in range(nbfig):
    … …  
    129144   v2=vini[i,:,:]
    130145   myv=myvar[i,:,:]
    131    u=np.copy(u2)
    132    v=np.copy(v2)
    133    myvarbis=np.copy(myv)
    134    print(i)
     146   u=swinglon(u2)
     147   v=swinglon(v2)
     148   myvarbis=swinglon(myv)
     149   print(i,"/",nbfig)
    135150   getfigvar(i)
    136151   mpl.savefig('mapwinds'+str('{0:03}'.format(i))+'.eps',dpi=200)
    … …  
    146161#mpl.subplots_adjust(hspace = .1)
    147162
     163
    148164#mpl.show()
    149165

  • trunk/LMDZ.PLUTO/util/script_figures/movie_winds/tsurf_winds.py

    r3823 r3833  
    1616filename1="../"+name+"_A.nc"
    1717filename2="../"+name+".nc"
    18 filename3="../../phisinit.nc"
     18filename3="../"+name+".nc"
     19# filename3="../../phisinit.nc"
    1920var="tsurf" #variable
    2021phisinit="phisinit" #variable
    … …  
    8788    pal=get_cmap(name="jet")
    8889    lev=np.linspace(37,51,15)
    89     # newlon=lon+180
    90     newlon=lon
     90    newlon=lon+180
    9191    CF=mpl.contourf(newlon, lat, myvarbis,lev,cmap=pal,extend='both')
    9292    yticks=[-90,-60,-30,0,30,60,90]
    … …  
    9898        t.set_fontsize(font)
    9999
    100     c=mpl.contour(newlon, lat, phisinit2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
     100    c=mpl.contour(newlon, lat, myvar2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
    101101    mpl.clabel(c, fmt='%2.1f',inline=1, colors='k', fontsize=23,inline_spacing=1)
    102102    #mpl.title('Local Time at Sputnik Planum='+str(i*3)+'H00',fontsize=font)
    … …  
    107107    getwinds(newlon,lat,u,v)
    108108
    109 # def getnumalt(choicealt,alt):
    110 #     numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
    111 #     numalt=numalt[0][0]
    112 #     return numalt
     109def getnumalt(choicealt,alt):
     110    numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
     111    numalt=numalt[0][0]
     112    return numalt
    113113
    114114#######################
    115 # numalt=getnumalt(30,alt)
    116 numalt=np.searchsorted(alt[...],30)
     115numalt=getnumalt(30,alt)
    117116print(('numalt =',numalt,'altitude=',alt[numalt]))
    118117uini=getvar(nc1,vari,tint,tim)[:,numalt]
    119118vini=getvar(nc1,varj,tint,tim)[:,numalt]
    120119myvar=getvar(nc2,var,tint,tim)
    121 phisinit2=getvar(nc3,phisinit) # phisinitmyvar2=phisinit2/0.6169/1000.  # altitude km
     120myvar2=getvar(nc3,phisinit) # phisinitmyvar2=myvar2/0.6169/1000.  # altitude km
    122121nbfig=uini.shape[0]
    123122print(("nbfig=",nbfig))
    124 phisinit2bis=np.copy(phisinit2)
     123myvar2bis=switchlon(myvar2)
    125124
    126125for i in range(nbfig):
    … …  
    129128   v2=vini[i,:,:]
    130129   myv=myvar[i,:,:]
    131    u=np.copy(u2)
    132    v=np.copy(v2)
    133    myvarbis=np.copy(myv)
    134    print(i)
     130   u=switchlon(u2)
     131   v=switchlon(v2)
     132   myvarbis=switchlon(myv)
     133   print(i,"/",nbfig)
    135134   getfigvar(i)
    136135   mpl.savefig('mapwinds'+str('{0:03}'.format(i))+'.eps',dpi=200)

  • trunk/LMDZ.PLUTO/util/script_figures/proftempNH.py

    r3823 r3833  
    11#! /usr/bin/env python
    2 from ppclass import pp
    32from    netCDF4               import    Dataset
    43from    numpy                 import    *
    … …  
    1110import datetime
    1211from mpl_toolkits.basemap import Basemap, shiftgrid
     12from FV3_utils import *
    1313
    1414############################
    15 zefile="diagfi2015_A.nc"
    16 d1="restart_ref/"
    17 
    18 f1=d1+zefile
    19 
    20 
    2115var="temperature" #variable
    2216
    23 nc1=Dataset(f1)
    24 lat=nc1.variables["lat"][:]
    25 lon=nc1.variables["lon"][:]
    26 alt=nc1.variables["altitude"][:]
    27 tim=nc1.variables["time_counter"][:]
     17nc1=Dataset(name+"_A.nc")
     18lat=getvar(nc1,"latitude")
     19lon=getvar(nc1,"longitude")
     20alt=getvar(nc1,"altitude")
     21tim=getvar(nc1,"Time")
    2822
    2923print(('Time = ',tim))
    … …  
    4337    print(('Point =',p,' Time=',tim[indt1]))
    4438    return indt1
    45  
    46 def getvar(filename,var):
    47     myvar = pp(file=filename,var=var,compute="nothing").getf()
    48     return myvar
    4939
    5040def getindex(lat,lon,mylat,mylon):
    … …  
    5747    indt=findindextime(tini,lt0,ltchoice,p)
    5848    indlat,indlon=getindex(lat,lon,p[1],p[0])
    59     myvar=getvar(f1,var)[indt,:,indlat,indlon]
     49    myvar=getvar(nc1,var)[indt,:,indlat,indlon]
    6050    return myvar
    6151############################
    … …  
    9282font=23
    9383
    94 alt=alt/1000.
    95 
    9684mpl.plot(myvar1,alt,'r',label='Entry')
    9785mpl.plot(myvar2,alt,'g',label='Entry-like on edge')

  • trunk/LMDZ.PLUTO/util/script_figures/proftempNH_SP.py

    r3823 r3833  
    1717
    1818nc1=Dataset(filename1)
    19 lat=nc1.variables["lat"][:]
    20 lon=nc1.variables["lon"][:]
    21 alt=nc1.variables["altitude"][:]
    22 tim=nc1.variables["time_counter"][:]
     19lat=getvar(nc1,"latitude")
     20lon=getvar(nc1,"longitude")
     21alt=getvar(nc1,"altitude")
     22tim=getvar(nc1,"Time")
    2323
    2424print(('Time = ',tim))

  • trunk/LMDZ.PLUTO/util/script_figures/proftempNH_obs.py

    r3823 r3833  
    2222
    2323nc1=Dataset(f1)
    24 lat=nc1.variables["lat"][:]
    25 lon=nc1.variables["lon"][:]
    26 alt=nc1.variables["altitude"][:]
    27 tim=nc1.variables["time_counter"][:]
     24lat=getvar(nc1,"latitude")
     25lon=getvar(nc1,"longitude")
     26alt=getvar(nc1,"altitude")
     27tim=getvar(nc1,"Time")
    2828
    2929print(('Time = ',tim))

  • trunk/LMDZ.PLUTO/util/script_figures/temp_section.py

    r3823 r3833  
    1212from mpl_toolkits.basemap import Basemap, shiftgrid
    1313from FV3_utils import *
    14 matplotlib.use('TKAgg')
     14from input import *
    1515
    1616############################
    … …  
    3333mpl.figure(figsize=(20, 10))
    3434
    35 myvar=getvar(nc1,var,tint,tim,l_mean=True,t_mean=True)
     35myvar=getvar(nc1,var,tint,l_mean=True,t_mean=True)
    3636font=26
    3737

  • trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom.py

    r3823 r3833  
    11#! /usr/bin/env python
    2 from ppclass import pp
    32from    netCDF4               import    Dataset
    43from    numpy                 import    *
    … …  
    1110import datetime
    1211from mpl_toolkits.basemap import Basemap, shiftgrid
     12from FV3_utils import *
    1313
    1414############################
    15 filename1="diagfi2015_S.nc"
     15filename1=name+"_S.nc"
    1616var="temperature" #variable
    17 xarea="-169,-165"
    18 yarea="-19,-15"
     17longitude=[-169,-165]
     18# longitude=-165
     19latitude=[-19,-15]
     20# latitude=-19
    1921
    2022# local time de la longitude consideree a t=0
    2123loct=12
    22 sol0=30
     24sol0=12
    2325t0=1./24*loct
    24 t1=t0+1
    25 tint=[str(sol0+t0)+','+str(sol0+t1)] #Time must be as written in the input file
     26t1=t0+12
     27tint=[sol0+t0,sol0+t1] #Time must be as written in the input file
    2628print(tint)
    2729nc1=Dataset(filename1)
    2830
    29 lat=nc1.variables["lat"][:]
    30 lon=nc1.variables["lon"][:]
    31 alt=nc1.variables["altitude"][:]
     31lat=getvar(nc1,"latitude")
     32lon=getvar(nc1,"longitude")
     33alt=getvar(nc1,"altitude")
     34tim=getvar(nc1,"Time",times=tint)
    3235############################
    33 
    34 def getvar(filename,var,tint,xarea,yarea):
    35     myvar = pp(file=filename,var=var,t=tint,x=xarea,y=yarea,compute="nothing").getf()
    36     print(('shape myvar = ',shape(myvar)))
    37     return myvar
    38 
    3936
    4037mpl.figure(figsize=(18, 10))
    4138
    42 
    43 myvar=getvar(filename1,var,tint,xarea,yarea)[:,:,0,0]
     39myvar=getvar(nc1,var,times=tint,longitudes=longitude,latitudes=latitude)[:,:,0,0]
    4440font=23
    45 tim=np.linspace(0,24,9)
     41# tim=np.linspace(0,24,9)
    4642print(("tim=",tim))
    4743print(('on prend les premiers indice, shape (tmps, alt, var) =',shape(tim), shape(alt), shape(myvar)))
    … …  
    4945pal=get_cmap(name="Spectral_r")
    5046lev=np.linspace(40,50,10)
    51 xticks=[0,2,4,6,8,10,12,14,16,18,20,22,24]
     47# xticks=[0,2,4,6,8,10,12,14,16,18,20,22,24]
    5248print(('hello:',np.linspace(0,24,13)))
    5349#yticks=np.linspace(0,240,9)
    54 alt=alt/1000.
    5550
    5651
    … …  
    7267mpl.xlabel('Local Time (h)',labelpad=10,fontsize=font)
    7368mpl.ylabel('Altitude (km)',labelpad=10, fontsize=font)
    74 mpl.xticks(xticks,fontsize=font)
    75 #mpl.xticks(fontsize=font)
     69# mpl.xticks(xticks,fontsize=font)
     70mpl.xticks(fontsize=font)
    7671#mpl.yticks(yticks,fontsize=font)
    7772mpl.yticks(fontsize=font)

  • trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom_1.py

    r3823 r3833  
    2323
    2424nc1=Dataset(f1)
    25 lat=nc1.variables["lat"][:]
    26 lon=nc1.variables["lon"][:]
    27 alt=nc1.variables["altitude"][:]
    28 tim=nc1.variables["time_counter"][:]
     25lat=getvar(nc1,"latitude")
     26lon=getvar(nc1,"longitude")
     27alt=getvar(nc1,"altitude")
     28tim=getvar(nc1,"Time")
    2929
    3030def findindextime(tini,lt0,lt1,p):

  • trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom_2.py

    r3823 r3833  
    2323
    2424nc1=Dataset(f1)
    25 lat=nc1.variables["lat"][:]
    26 lon=nc1.variables["lon"][:]
    27 alt=nc1.variables["altitude"][:]
    28 tim=nc1.variables["time_counter"][:]
     25lat=getvar(nc1,"latitude")
     26lon=getvar(nc1,"longitude")
     27alt=getvar(nc1,"altitude")
     28tim=getvar(nc1,"Time")
    2929
    3030def findindextime(tini,lt0,lt1,p):

  • trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom_3.py

    r3823 r3833  
    2323
    2424nc1=Dataset(f1)
    25 lat=nc1.variables["lat"][:]
    26 lon=nc1.variables["lon"][:]
    27 alt=nc1.variables["altitude"][:]
    28 tim=nc1.variables["time_counter"][:]
     25lat=getvar(nc1,"latitude")
     26lon=getvar(nc1,"longitude")
     27alt=getvar(nc1,"altitude")
     28tim=getvar(nc1,"Time")
    2929
    3030def findindextime(tini,lt0,lt1,p):

  • trunk/LMDZ.PLUTO/util/script_figures/temp_time_zoom_ref.py

    r3823 r3833  
    2323
    2424nc1=Dataset(f1)
    25 lat=nc1.variables["lat"][:]
    26 lon=nc1.variables["lon"][:]
    27 alt=nc1.variables["altitude"][:]
    28 tim=nc1.variables["time_counter"][:]
     25lat=getvar(nc1,"latitude")
     26lon=getvar(nc1,"longitude")
     27alt=getvar(nc1,"altitude")
     28tim=getvar(nc1,"Time")
    2929
    3030def findindextime(tini,lt0,lt1,p):

  • trunk/LMDZ.PLUTO/util/script_figures/tempglobmean.py

    r3823 r3833  
    1111from mpl_toolkits.basemap import Basemap, shiftgrid
    1212from FV3_utils import *
     13from input import *
    1314
    1415############################
    … …  
    5253mpl.figure(figsize=(8, 10))
    5354
    54 myvar=getvar(nc1,var,tint,tim,t_mean=True)
    55 #myvar2=getvar(nc1,var,tint,tim)
     55myvar=getvar(nc1,var,tint,t_mean=True)
     56#myvar2=getvar(nc1,var,tint)
    5657aire = getvar(nc2,"aire")
    5758totarea=zetotarea(aire)
    … …  
    7071#mpl.plot(meantemp2,alt2,'b--')
    7172
    72 #mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
     73mpl.title('Global mean temperature',fontsize=font)
    7374mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
    7475mpl.xlabel('Temperature (K)',labelpad=10, fontsize=font)
    … …  
    7879mpl.yticks(fontsize=font)
    7980mpl.grid()
    80 mpl.legend(["Ref","Alt"],prop={'size':27},loc='upper left')
     81# mpl.legend(["Ref","Alt"],prop={'size':27},loc='upper left')
    8182
    8283left  = 0.2  # the left side of the subplots of the figure
Note: See TracChangeset for help on using the changeset viewer.