Changeset 569 for trunk


Ignore:
Timestamp:
Mar 7, 2012, 4:56:38 PM (13 years ago)
Author:
acolaitis
Message:

PYTHON
######

=========
myplot.py
=========

  • Added missing value handling for API files (+1e36)
  • Added check_localtime function that re-order the time axis in localtime mode

=========
planetoplot.py
=========

  • Corrected a problem when trying to plot a contour or unidim plot in localtime along a cyclic time axis (i.e. when the day was changing).

    one day is now numbered from 0 to 24 and the rest is either negative (below 0) or greater than 24.

  • Localtimes are no longer integers, which was causing problems when using files with more than 1 localtime per hour.
  • Corrected a very wierd bug where the asked value for --lon and --lat was changing when plotting along different files. This was due to a spooky pointer problem where the value of a variable was changing when an other one was changed !!!
  • Corrected a bug with localtimes where there was a drift in initial localtimes when using different files and when not using the --operation cat option.
  • Added the possibility to ask for the variable 'UV' (or 'uv' or 'uvmet' if using API) which will automatically compute the norm of horizontal winds 'U' and 'V' (or 'u' and 'v' or 'Um' and 'Vm' if using API). This also work when using reference file. This will only compute the norm if no variable with this name already exists in the file.
  • Added the possibility to call overcontours in plot with a reference file (--fref option) using the usual -w ... syntax. The overcontour will also be displayed on the "operation" plot, without being processed.
  • Added the possibility of specifying manually labels for unidim curves (see below)

=========
myscript.py
=========

  • New option --labels "label 1","label 2",......,"label N" will label in this order each curve of a multi-file unidim plot, with "label N" corresponding to file N.

=========
pp.py
=========

  • Added calls to API for reference file when using --fref so that the reference file is also processed.
Location:
trunk/UTIL/PYTHON
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/myplot.py

    r568 r569  
    3535    ltst = ltst % 24
    3636    return ltst
     37
     38## Author: AC
     39def check_localtime(time):
     40    a=-1
     41    for i in range(len(time)-1):
     42       if time[i] > time[i+1]: a=i
     43    if a >= 0:
     44       print "Sorry, time axis is not regular."
     45       print "Contourf needs regular axis... recasting"
     46       for i in range(a+1):
     47          time[i]=time[i]-24.
     48    return time
    3749
    3850## Author: AS, AC, JL
    … …  
    7183       mask = np.ma.getmask(masked)
    7284       if (True in np.array(mask)):out=masked
    73        else:out=field
     85       else:
     86       # missing values from api are 1e36
     87          masked=np.ma.masked_where(field > 1e35,field)
     88          masked.set_fill_value([np.NaN])
     89          mask = np.ma.getmask(masked)
     90          if (True in np.array(mask)):out=masked
     91          else:out=field
    7492    return out
    7593

  • trunk/UTIL/PYTHON/myscript.py

    r562 r569  
    3838    parser.add_option('--xlabel',       action='store',dest='xlab',       type="string",  default=None, help='customize the x-axis label')
    3939    parser.add_option('--ylabel',       action='store',dest='ylab',       type="string",  default=None, help='customize the y-axis label')
     40    parser.add_option('--labels',       action='store',dest='labels',     type="string",  default=None, help='customize 1D curve labels. String separated by commas. [None]')
    4041
    4142    ### SPECIFIC FOR MAPPING [MAPMODE 1]

  • trunk/UTIL/PYTHON/planetoplot.py

    r562 r569  
    6262           seevar=False,\
    6363           xlab=None,\
    64            ylab=None):
     64           ylab=None,\
     65           lbls=None):
    6566
    6667    ####################################################################################################################
    … …  
    7374                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
    7475                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
    75                        getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds
     76                       getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds
    7677    from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
    7778    import matplotlib as mpl
    … …  
    139140      ### ... VAR
    140141      varname=var[vvv]
    141       if varname not in nc.variables:
     142      if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet']))):
    142143          if len(varinfile) == 1:   varname = varinfile[0]
    143144          else:                     varname = False
    … …  
    196197          ###### principle: calculate correct indices then repopulate slon and slat
    197198          if slon is not None or slat is not None:
     199              if slon is not None: sslon = np.zeros([1,2])
     200              if slat is not None: sslat = np.zeros([1,2])
     201
    198202              if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile:   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
    199203              else:                                                                           iwantawhereplot = None   #do not show anything, just select indices
    … …  
    208212                 indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot) 
    209213                 lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] )
    210                  #print vlon, lonp, vlat, latp
    211214              for iii in range(numlon):
    212215               for jjj in range(numlat):
    213                  if slon is not None: slon[iii][0] = indices[iii,0,1] ; slon[iii][1] = indices[iii,0,1]  #...this is idx
    214                  if slat is not None: slat[jjj][0] = indices[0,jjj,0] ; slat[jjj][1] = indices[0,jjj,0]  #...this is idy
     216                 if slon is not None: sslon[iii][0] = indices[iii,0,1] ; sslon[iii][1] = indices[iii,0,1]  #...this is idx
     217                 if slat is not None: sslat[jjj][0] = indices[0,jjj,0] ; sslat[jjj][1] = indices[0,jjj,0]  #...this is idy
    215218              lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
    216219          ######
    … …  
    235238              elif "Time" in nc.variables:  time = count + np.arange(0,len(nc.variables["Time"]),1)
    236239              else:                         time = count + np.arange(0,1,1)
    237               if nnn > 0:  count = time[-1] + 1  ## so that a cat is possible with simple subscripts
    238               else:        count = 0
     240              if ope in ["cat"]:
     241                 if nnn > 0:  count = time[-1] + 1  ## so that a cat is possible with simple subscripts
     242                 else:        count = 0
     243              else: count = 0
    239244          if axtime in ["lt"]:
    240               for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
     245              ftime = np.zeros(len(time))
     246              for i in range(len(time)):
     247                 ftime[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
     248              time=ftime
     249              time=check_localtime(time)
    241250              print "LOCAL TIMES.... ", time
     251
    242252          ###
    243253          if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
    … …  
    271281      else:
    272282      ##### GENERAL STUFF HERE
    273           all_var[k] = getfield(nc,varname)
    274      
     283          if ((varname in ['UV','uv','uvmet']) and (typefile in ['meso']) and (varname not in nc.variables)):
     284             if "U" in varinfile: zu=getfield(nc,'U')
     285             elif "Um" in varinfile: zu=getfield(nc,'Um')
     286             if "V" in varinfile: zv=getfield(nc,'V')
     287             elif "Vm" in varinfile: zv=getfield(nc,'Vm')
     288             znt,znz,zny,znx = np.array(zu).shape
     289             if "U" in varinfile:znx=znx-1
     290             zuint = np.zeros([znt,znz,zny,znx])
     291             zvint = np.zeros([znt,znz,zny,znx])
     292             if "U" in varinfile:
     293                for xx in np.arange(znx):
     294                    zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
     295                for yy in np.arange(zny):
     296                    zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
     297             else:
     298                    zuint=zu
     299                    zvint=zv
     300             all_var[k] = np.sqrt(zuint**2 + zvint**2)
     301          else:
     302             all_var[k] = getfield(nc,varname)
     303
    275304      print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
    276305      #### End of for namefile in namefiles
    … …  
    283312             if len(var) > 1: errormess("for this operation... please set only one var !")
    284313             if ope in ["-","+","-%"]:
    285                 if fileref is not None:   all_var[k] = getfield(Dataset(fileref),all_varname[k-1]) ; all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
    286                 else:                     errormess("fileref is missing!")
     314                if fileref is not None:   
     315
     316                   if ((all_varname[k-1] in ['UV','uv','uvmet']) and (typefile in ['meso']) and (all_varname[k-1] not in Dataset(fileref).variables)):
     317                      if "U" in varinfile: zu=getfield(Dataset(fileref),'U')
     318                      elif "Um" in varinfile: zu=getfield(Dataset(fileref),'Um')
     319                      if "V" in varinfile: zv=getfield(Dataset(fileref),'V')
     320                      elif "Vm" in varinfile: zv=getfield(Dataset(fileref),'Vm')
     321                      znt,znz,zny,znx = np.array(zu).shape
     322                      if "U" in varinfile:znx=znx-1
     323                      zuint = np.zeros([znt,znz,zny,znx])
     324                      zvint = np.zeros([znt,znz,zny,znx])
     325                      if "U" in varinfile:
     326                         for xx in np.arange(znx):
     327                             zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
     328                         for yy in np.arange(zny):
     329                             zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
     330                      else:
     331                         zuint=zu
     332                         zvint=zv
     333                      all_var[k] = np.sqrt(zuint**2 + zvint**2)
     334                   else:
     335                      all_var[k] = getfield(Dataset(fileref),all_varname[k-1])
     336         
     337                   all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1]
     338
     339                else: errormess("fileref is missing!")
    287340                if ope == "-":     all_var[k+1]= all_var[k-1] - all_var[k]
    288341                elif ope == "+":   all_var[k+1]= all_var[k-1] + all_var[k]
    … …  
    291344                    masked.set_fill_value([np.NaN])
    292345                    all_var[k+1]= 100.*(all_var[k-1] - masked)/masked
    293                 all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; numplot = numplot+2
     346                all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; all_var2[k+1] = all_var2[k] ; numplot = numplot+2
    294347             elif ope in ["cat"]:
    295348                tabtime = all_time[0];tab = all_var[0];k = 1
    … …  
    329382       ## get all indexes to be taken into account for this subplot and then reduce field
    330383       ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
    331        indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
    332        indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
     384       indexlon  = getsindex(sslon,(nplot-1)%nlon,lon)
     385       indexlat  = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
    333386       indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
    334387       if ope is not None:
    … …  
    456509
    457510                    if which == "unidim":
     511                        #lbls = ["<TH+RiSL+MY4> tau=0.5","<Convadj> tau=0.5","<TH+RiSL+MY4> tau=1","<Convadj> tau=1"]
    458512                        lbl = ""
    459513                        if indexlat is not None:  lbl = lbl + " ix" + str(indexlat[0])
    … …  
    462516                        if indextime is not None: lbl = lbl + " it" + str(indextime[0])
    463517                        if lbl == "": lbl = namefiles[index_f]
     518
     519                        if lbls is not None: lbl=lbls[index_f]
     520
    464521                        if mrate is not None: x = y  ## because swapaxes...
    465522                        #what_I_plot_frame = np.diff(what_I_plot_frame, n=1) ; x = x[1:]
    … …  
    473530                        if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot*1000+imov)+'.txt')
    474531
    475                     elif which == "regular": 
     532                    elif which == "regular":
    476533                        if hole:         what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
    477534                        else:            what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
    … …  
    570627       if mrate is not None: basename = "movie_" + basename
    571628       if typefile in ['meso']:
    572             if slon is not None: basename = basename + "_lon_" + str(int(round(lonp)))
    573             if slat is not None: basename = basename + "_lat_" + str(int(round(latp)))
     629            if sslon is not None: basename = basename + "_lon_" + str(int(round(lonp)))
     630            if sslat is not None: basename = basename + "_lat_" + str(int(round(latp)))
    574631            plottitle = basename
    575632            ### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)

  • trunk/UTIL/PYTHON/pp.py

    r562 r569  
    8686            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
    8787
    88         #########################################################
    89         ### Call Fortran routines for vertical interpolations ###     
    9088        #########################################################
    9189        if opt.itp is not None:
    … …  
    114112                if fff == 0: zetab = newname
    115113                else:        zetab = np.append(zetab,newname)
     114            if interpref:
     115                reffile = api_onelevel (  path_to_input   = '', \
     116                                               input_name      = opt.fref, \
     117                                               fields          = zefields, \
     118                                               interp_method   = opt.itp, \
     119                                               interp_level    = ze_interp_levels, \
     120                                               onelevel        = zelevel, \
     121                                               nocall          = opt.nocall )
    116122            zefiles = zetab #; print zefiles
    117123            zelevel = 0 ## so that zelevel could play again the role of nvert
    … …  
    158164                blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\
    159165                mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\
    160                 redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab)
     166                redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels))
    161167        print 'DONE: '+name
    162168        system("rm -f to_be_erased")
Note: See TracChangeset for help on using the changeset viewer.