source: trunk/UTIL/PYTHON/meso.py @ 388

#!/usr/bin/env python

### A. Spiga

###########################################################################################
###########################################################################################
### What is below relate to running the file as a command line executable (very convenient)
if __name__ == "__main__":
    import sys
    from optparse import OptionParser    ### to be replaced by argparse
    from api_wrapper import api_onelevel
    from netCDF4 import Dataset
    from myplot import getlschar, separatenames, readslices, adjust_length
    from os import system
    from planetoplot import planetoplot
    from myscript import getparseroptions
    import numpy as np

    #############################
    ### Get options and variables
    parser = OptionParser()
    getparseroptions(parser)
    (opt,args) = parser.parse_args()

    if opt.file is None: 
        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
        exit()
    if opt.var is None and opt.anomaly is True:
        print "Cannot ask to compute anomaly if no variable is set"
        exit()
    print "Options:", opt

    listvar = '' 
    if opt.var is None:
        zerange = [-999999]
    else:
        zelen = len(opt.var)
        zerange = range(zelen)
        #if zelen == 1: listvar = opt.var[0] + ','
        #else         : 
        for jjj in zerange: listvar += opt.var[jjj] + ','
        listvar = listvar[0:len(listvar)-1]
        vmintab = adjust_length (opt.vmin, zelen)  
        vmaxtab = adjust_length (opt.vmax, zelen)


    ################################
    ### General check
 
    if opt.fref is not None:
       if opt.operat is None:
          print "you must specify an operation when using a reference file"
          exit()
    if opt.operat in ["+","-"]:
       if opt.fref is None:
          print "you must specifiy a reference file when using inter-file operations"
          exit()

    interpref=False
    if opt.fref is not None:
       if opt.operat is not None:
          if opt.itp is not None:
             interpref=True

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


    print "file, length", opt.file, len(opt.file)

    zeslat  = readslices(opt.slat)
    zeslon  = readslices(opt.slon)
    zesvert = readslices(opt.svert)
    zestime = readslices(opt.stime)
    print "slat,zeslat", opt.slat, zeslat
    print "slon,zeslon", opt.slon, zeslon
    print "svert,zesvert", opt.svert, zesvert
    print "stime,zestime", opt.stime, zestime

    for i in range(len(opt.file)):

        zefiles = separatenames(opt.file[i])
        print "zefiles", zefiles

        #zefile = opt.file[i]
        #print zefile    
        zefile = zefiles[0]
              
        #zelevel = opt.lvl   
        stralt = None
        [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
    
        inputnvert = separatenames(opt.lvl)
        if np.array(inputnvert).size == 1:
            zelevel = float(inputnvert[0])
            ze_interp_levels = [-9999.]
        else:
            zelevel = -99.
            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
        print 'level: ', zelevel
        print 'interp_levels: ',ze_interp_levels

        #####################################################
        ### Call Fortran routines for vertical interpolations        
        if opt.itp is not None:
            if zelevel == 0. and opt.itp == 4:  zelevel = 0.010
            ### winds or no winds
            if opt.winds            :  zefields = 'uvmet'
            else                    :  zefields = ''
            ### var or no var
            #if opt.var is None      :  pass
            if zefields == ''       :  zefields = listvar 
            else                    :  zefields = zefields + "," + listvar 
            if opt.var2 is not None : zefields = zefields + "," + opt.var2  
            print zefields
            zefile = api_onelevel (  path_to_input   = '', \
                                     input_name      = zefile, \
                                     fields          = zefields, \
                                     interp_method   = opt.itp, \
                                     interp_level    = ze_interp_levels, \
                                     onelevel        = zelevel, \
                                     nocall          = opt.nocall )
            print zefile
            zelevel = 0 ## so that zelevel could play again the role of nvert

        if opt.var is None: zerange = [-999999]
        else:               zerange = range(zelen) 








        if interpref:
           interpolated_ref=""
           interpolated_ref=call_zrecast(interp_mode=opt.itp,\
                    input_name=[opt.fref],\
                    fields=zevars)

           reffile=interpolated_ref[0]
        else:
           reffile=opt.fref
# Divers ####################################################
    
        zexaxis=[opt.xmin,opt.xmax]
        zeyaxis=[opt.ymin,opt.ymax]

        for jjj in zerange:
            if jjj == -999999: 
                argvar  = None
                zevmin = None
                zevmax = None
            else:
                argvar = opt.var[jjj]
                if vmintab[jjj] != -999999:  zevmin = vmintab[jjj]
                else:                        zevmin = None
                if vmaxtab[jjj] != -999999:  zevmax = vmaxtab[jjj] 
                else:                        zevmax = None

            #############
            ### Main call
            name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=argvar,\
                numplot=opt.num,colorb=opt.clb,winds=opt.winds,\
                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
                itstep=opt.it,hole=opt.hole,save=opt.save,\
                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.frt,\
                mult=opt.mult,zetitle=opt.zetitle,\
                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
                outputname=opt.out,resolution=opt.res,\
                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
                blat=opt.blat,tsat=opt.tstat)
        print 'Done: '+name
        system("rm -f to_be_erased")
  
    #########################################################
    ### Generate a .sh file with the used command saved in it
    command = ""
    for arg in sys.argv: command = command + arg + ' '
    f = open(name+'.sh', 'w')
    f.write(command)