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

Last change on this file since 673 was 647, checked in by tnavarro, 13 years ago

Corrected bug in reducefield for masked arrays with grid area. Possibility to see stream function for a lat/vert slice with --stream option. Output of both data and axis with save txt option in 1D. Small bug corrected: title is the file name for multiple plots with multiple plots. Cleanup in myplot.py

  • 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),\
[647]174                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),streamflag=opt.stream)
[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.