source: trunk/UTIL/PYTHON/pp.py @ 641

Last change on this file since 641 was 612, checked in by acolaitis, 13 years ago

PYTHON
======

New options:

i) --mark lon,lat will superimpose a cross at lon,lat (floats) on a longitude latitude plot (mapmode 1)
ii) --lstyle let you customize linestyles for 1D plots. Works in the same way as --labels. Exemple: --lstyle "--r","-r","--b","-b"
iii) One can now ask for the variable -v slopexy or SLOPEXY to plot (or overplot with --var2 slopexy) the amplitude of the slope in mesoscale simulations
iv) One can now ask for the variable -v deltat or DELTAT to plot the temperature difference between surface temperature and an arbitrary model level. This automatically calls API with "tk,tsurf". One should use this in conjunction with -i 4 -l xxx so that the difference is made between 0m (surface) and xxx m.

Minor corrections:

Verification plot showing longitude latitude of a slice in a HGT map now also requires "display" to be true.
Winds can now be overplotted in --operation - maps.

  • Property svn:executable set to *
File size: 9.4 KB
RevLine 
[349]1#!/usr/bin/env python
2
[392]3### A. Spiga + T. Navarro + A. Colaitis
[349]4
5###########################################################################################
6###########################################################################################
7### What is below relate to running the file as a command line executable (very convenient)
8if __name__ == "__main__":
9    import sys
10    from optparse import OptionParser    ### to be replaced by argparse
11    from api_wrapper import api_onelevel
[464]12    from gcm_transformations import call_zrecast
[349]13    from netCDF4 import Dataset
[392]14    from myplot import getlschar, separatenames, readslices, adjust_length, whatkindfile, errormess
[349]15    from os import system
16    from planetoplot import planetoplot
[377]17    from myscript import getparseroptions
[478]18    import glob
[349]19    import numpy as np
20
[478]21
[349]22    #############################
23    ### Get options and variables
[392]24    parser = OptionParser() ; getparseroptions(parser) ; (opt,args) = parser.parse_args()
25    if opt.file is None:                                errormess("I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h")
26    if opt.var is None and opt.anomaly is True:         errormess("Cannot ask to compute anomaly if no variable is set")
27    if opt.fref is not None and opt.operat is None:     errormess("you must specify an operation when using a reference file")
28    if opt.operat in ["+","-"] and opt.fref is None:    errormess("you must specifiy a reference file when using inter-file operations")
29    if opt.fref is not None and opt.operat is not None and opt.itp is not None: interpref=True
30    else:   interpref=False
[451]31    if opt.rate is not None:      opt.save = "avi"
32    elif opt.save == "avi":       opt.rate = 8   ## this is a default value for -S avi
33    if opt.save == "html":        opt.rate = -1  ## this is convenient because everything is done in planetoplot with mrate
[349]34
[392]35    #############################
36    ### Get infos about slices
37    zeslat  = readslices(opt.slat) ; zeslon  = readslices(opt.slon) ; zesvert = readslices(opt.svert) ; zestime = readslices(opt.stime)
38    reffile = opt.fref
39    zexaxis = [opt.xmin,opt.xmax] ; zeyaxis=[opt.ymin,opt.ymax]
[349]40
[392]41    #############################
[478]42    ### Catch multiple files
[479]43    if "*" in opt.file[0] or "?" in opt.file[0]: 
44        yeah = glob.glob(opt.file[0]) ; yeah.sort()
45        opt.file[0] = yeah[0] 
[489]46        for file in yeah[1:]: opt.file[0] = opt.file[0] + "," + file
[478]47
48    #############################
[392]49    ### 1. LOOP ON FILE LISTS TO BE PUT IN DIFFERENT FIGURES
50    for i in range(len(opt.file)):
[349]51
[392]52      zefiles = separatenames(opt.file[i])
[380]53
[426]54      typefile = whatkindfile(Dataset(zefiles[0])) ; stralt = None
[562]55      if typefile in ["meso"]:         
[392]56          [lschar,zehour,zehourin] = getlschar ( zefiles[0] )
[426]57          if opt.var is None:  opt.var = ["HGT"] ; opt.clb = "nobar"
[429]58      elif typefile in ["geo"]:
[426]59          [lschar,zehour,zehourin] = ["",0,0]
60          if opt.var is None:  opt.var = ["HGT_M"] ; opt.clb = "nobar"
[392]61      else:                                       
62          [lschar,zehour,zehourin] = ["",0,0]
[426]63          if opt.var is None:  opt.var = ["phisinit"] ; opt.clb = "nobar"
[380]64
[392]65      if opt.vmin is not None : zevmin  = opt.vmin[min(i,len(opt.vmin)-1)]
66      else:                     zevmin = None
67      if opt.vmax is not None : zevmax  = opt.vmax[min(i,len(opt.vmax)-1)]
68      else:                     zevmax = None
69      #print "vmin, zevmin", opt.vmin, zevmin ; print "vmax, zevmax", opt.vmax, zevmax
[380]70
[392]71      #############################
72      ### 2. LOOP ON VAR LISTS TO BE PUT IN DIFFERENT FIGURES
73      for j in range(len(opt.var)):
[380]74
[392]75        zevars = separatenames(opt.var[j])
[380]76
[377]77        inputnvert = separatenames(opt.lvl)
[351]78        if np.array(inputnvert).size == 1:
79            zelevel = float(inputnvert[0])
80            ze_interp_levels = [-9999.]
[610]81        elif np.array(inputnvert).size > 2:
[489]82            zelevel = -99.
[610]83            start = float(inputnvert[0])
84            stop = float(inputnvert[1])
85            if np.array(inputnvert).size == 2:  numsample = 20
86            else:                               numsample = float(inputnvert[2])
87            if stop > start:   
88               ## altitude coordinates
89               ze_interp_levels = np.linspace(start,stop,numsample)
90            else:
91               ## pressure coordinates
92               ze_interp_levels = np.logspace(np.log10(start),np.log10(stop),numsample)
[351]93
[392]94        #########################################################
[377]95        if opt.itp is not None:
[562]96         if opt.itp > 0:
[392]97          #####
98          ##### MESOSCALE : written by AS
99          #####
[562]100          if typefile in ["meso"]:
[377]101            if zelevel == 0. and opt.itp == 4:  zelevel = 0.010
[349]102            ### winds or no winds
103            if opt.winds            :  zefields = 'uvmet'
[612]104            elif opt.var[j] in ['deltat','DELTAT'] : zefields = 'tk,TSURF'
[349]105            else                    :  zefields = ''
106            ### var or no var
[392]107            if zefields == ''       :  zefields = opt.var[j] 
108            else                    :  zefields = zefields + "," + opt.var[j]
109            if opt.var2 is not None :  zefields = zefields + "," + opt.var2 
110            ### call fortran routines
111            for fff in range(len(zefiles)):
112                newname = api_onelevel (  path_to_input   = '', \
113                                               input_name      = zefiles[fff], \
114                                               fields          = zefields, \
115                                               interp_method   = opt.itp, \
116                                               interp_level    = ze_interp_levels, \
117                                               onelevel        = zelevel, \
118                                               nocall          = opt.nocall )
119                if fff == 0: zetab = newname
120                else:        zetab = np.append(zetab,newname)
[569]121            if interpref:
122                reffile = api_onelevel (  path_to_input   = '', \
123                                               input_name      = opt.fref, \
124                                               fields          = zefields, \
125                                               interp_method   = opt.itp, \
126                                               interp_level    = ze_interp_levels, \
127                                               onelevel        = zelevel, \
128                                               nocall          = opt.nocall )
[392]129            zefiles = zetab #; print zefiles
[349]130            zelevel = 0 ## so that zelevel could play again the role of nvert
[392]131          #####
132          ##### GCM : written by AC
133          #####
134          elif typefile == "gcm":
[475]135            inputvar = zevars
136            if opt.var2 is not None : inputvar = np.append(inputvar,opt.var2)
[392]137            interpolated_files=""
138            interpolated_files=call_zrecast(interp_mode=opt.itp,\
139                    input_name=zefiles,\
[475]140                    fields=inputvar,\
[489]141                    limites = ze_interp_levels,\
[464]142                    predefined=opt.intas)
[349]143
[392]144            zefiles=interpolated_files
145            if interpref:
146               interpolated_ref=""
147               interpolated_ref=call_zrecast(interp_mode=opt.itp,\
[380]148                    input_name=[opt.fref],\
[409]149                    fields=zevars,\
[464]150                    predefined=opt.intas)
[380]151
[392]152               reffile=interpolated_ref[0]
153          else:
[429]154            print "type not supported"
[392]155            exit()
[377]156
[392]157        #############
158        ### Main call
159        name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
160                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevars,\
[424]161                colorb=opt.clb,winds=opt.winds,\
[349]162                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
163                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
[424]164                hole=opt.hole,save=opt.save,\
165                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,\
[578]166                mult=opt.mult,zetitle=separatenames(opt.zetitle),\
[359]167                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
[377]168                outputname=opt.out,resolution=opt.res,\
[380]169                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
[385]170                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
[453]171                blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\
[483]172                mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\
[612]173                redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels),\
174                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark))
[392]175        print 'DONE: '+name
[349]176        system("rm -f to_be_erased")
177 
178    #########################################################
179    ### Generate a .sh file with the used command saved in it
[392]180    command = "" 
[349]181    for arg in sys.argv: command = command + arg + ' '
[424]182    #if typefile not in ["meso","mesoapi"]: name = 'pycommand'
[451]183    if opt.save == "gui":    name = 'pycommand'
184    elif opt.save == "avi":  system("mv -f movie*.avi "+name+".avi")
185    elif opt.save == "html": system("cat $PYTHONPATH/header.html > anim.html ; cat zepics >> anim.html ; cat $PYTHONPATH/body.html >> anim.html ; rm -rf zepics "+name+" ; mkdir "+name+" ; mv anim.html image*png "+name) 
[349]186    f = open(name+'.sh', 'w')
187    f.write(command)
[392]188
[468]189    #print "********** OPTIONS: ", opt
[392]190    print "********************************************************** END"
Note: See TracBrowser for help on using the repository browser.