Changeset 405


Ignore:
Timestamp:
Nov 21, 2011, 6:27:42 PM (13 years ago)
Author:
aslmd
Message:

GRAPHICS: 1. handling of anomaly in profiles as well. 2. fixed xmin ymin etc... that were not working with 1D plot. 3. fixed localtime for mesoscale files, still a bit of work though.

Location:
trunk/UTIL/PYTHON
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/myplot.py

    r400 r405  
    6767
    6868## Author: AS + TN + AC
    69 def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None):
     69def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False):
    7070    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
    7171    ### it would be actually better to name d4 d3 d2 d1 as t z y x
     72    ### ... note, anomaly is only computed over d1 and d2 for the moment
    7273    import numpy as np
    7374    from mymath import max,mean
     75    csmooth = 12 ## a fair amount of grid points
    7476    dimension = np.array(input).ndim
    7577    shape = np.array(input).shape
    7678    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
    7779    print 'IN  REDUCEFIELD dim,shape: ',dimension,shape
     80    if anomaly: print 'ANOMALY ANOMALY'
    7881    output = input
    7982    error = False
    … …  
    110113             output = mean(input[d4,:,:,:],axis=0)
    111114             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     115             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
    112116             output = mean(output[d2,:],axis=0)
    113117             output = mean(output[d1],axis=0)
    … …  
    115119             output = mean(input[d4,:,:,:],axis=0)
    116120             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     121             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
    117122             output = mean(output[d2,:],axis=0)
    118123        elif d4 is not None and d3 is not None and d1 is not None:
    119124             output = mean(input[d4,:,:,:],axis=0)
    120125             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     126             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
    121127             output = mean(output[:,d1],axis=1)
    122128        elif d4 is not None and d2 is not None and d1 is not None:
    123129             output = mean(input[d4,:,:,:],axis=0)
     130             if anomaly:
     131                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
    124132             output = mean(output[:,d2,:],axis=1)
    125133             output = mean(output[:,d1],axis=1)
     134             #noperturb = smooth1d(output,window_len=7)
     135             #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
     136             #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
    126137        elif d3 is not None and d2 is not None and d1 is not None:
    127138             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
     139             if anomaly:
     140                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
    128141             output = mean(output[:,d2,:],axis=1)
    129142             output = mean(output[:,d1],axis=1)
    … …  
    131144             output = mean(input[d4,:,:,:],axis=0)
    132145             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
     146             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
    133147        elif d4 is not None and d2 is not None: 
    134148             output = mean(input[d4,:,:,:],axis=0)
    … …  
    248262    from numpy import array
    249263    from string import rstrip
    250     namefile = rstrip( rstrip( namefile, chars="_z"), chars="_zabg")
     264    namefile = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
    251265    #### we assume that wrfout is next to wrfout_z and wrfout_zabg
    252266    nc  = Dataset(namefile)
    … …  
    394408        [lon2d,lat2d] = [lon,lat]
    395409    return lon2d,lat2d
     410
     411## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
     412def smooth1d(x,window_len=11,window='hanning'):
     413    import numpy
     414    """smooth the data using a window with requested size.
     415    This method is based on the convolution of a scaled window with the signal.
     416    The signal is prepared by introducing reflected copies of the signal
     417    (with the window size) in both ends so that transient parts are minimized
     418    in the begining and end part of the output signal.
     419    input:
     420        x: the input signal
     421        window_len: the dimension of the smoothing window; should be an odd integer
     422        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
     423            flat window will produce a moving average smoothing.
     424    output:
     425        the smoothed signal
     426    example:
     427    t=linspace(-2,2,0.1)
     428    x=sin(t)+randn(len(t))*0.1
     429    y=smooth(x)
     430    see also:
     431    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
     432    scipy.signal.lfilter
     433    TODO: the window parameter could be the window itself if an array instead of a string   
     434    """
     435    x = numpy.array(x)
     436    if x.ndim != 1:
     437        raise ValueError, "smooth only accepts 1 dimension arrays."
     438    if x.size < window_len:
     439        raise ValueError, "Input vector needs to be bigger than window size."
     440    if window_len<3:
     441        return x
     442    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
     443        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
     444    s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
     445    #print(len(s))
     446    if window == 'flat': #moving average
     447        w=numpy.ones(window_len,'d')
     448    else:
     449        w=eval('numpy.'+window+'(window_len)')
     450    y=numpy.convolve(w/w.sum(),s,mode='valid')
     451    return y
    396452
    397453## Author: AS
    … …  
    607663def fmtvar (whichvar="def"):
    608664    fmtvar    =     { \
    609              "tk":           "%.0f",\
     665             "TK":           "%.0f",\
    610666             "T_NADIR_DAY":  "%.0f",\
    611667             "T_NADIR_NIT":  "%.0f",\
    612668             "TEMP_DAY":     "%.0f",\
    613669             "TEMP_NIGHT":   "%.0f",\
    614              "tpot":         "%.0f",\
     670             "TPOT":         "%.0f",\
    615671             "TSURF":        "%.0f",\
    616672             "def":          "%.1e",\
    … …  
    623679             "VMR_ICE":      "%.1e",\
    624680             "MTOT":         "%.1f",\
    625              "anomaly":      "%.1f",\
     681             "ANOMALY":      "%.1f",\
    626682             "W":            "%.1f",\
    627683             "WMAX_TH":      "%.1f",\
    628684             "QSURFICE":     "%.0f",\
    629              "Um":           "%.0f",\
     685             "UM":           "%.0f",\
    630686             "ALBBARE":      "%.2f",\
    631687             "TAU":          "%.1f",\
    … …  
    650706             "def":          "spectral",\
    651707             "HGT":          "spectral",\
    652              "tk":           "gist_heat",\
     708             "TK":           "gist_heat",\
    653709             "TSURF":        "RdBu_r",\
    654710             "QH2O":         "PuBu",\
    … …  
    666722             "W":            "jet",\
    667723             "WMAX_TH":      "spectral",\
    668              "anomaly":      "RdBu_r",\
     724             "ANOMALY":      "RdBu_r",\
    669725             "QSURFICE":     "hot_r",\
    670726             "ALBBARE":      "spectral",\

  • trunk/UTIL/PYTHON/planetoplot.py

    r402 r405  
    258258           if len(var) != 1 or len(namefiles) != 1: indextime = first
    259259       else:
    260            indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) 
     260           indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
     261           if indextime is not None: ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) ) 
    261262
    262263       if fileref is not None:      index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2)  ## OK only 1 var, see test in the beginning
    … …  
    292293       varname = all_varname[index_f]
    293294       if varname:
    294            what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert , yint=yintegral, alt=vert)
     295           what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert , yint=yintegral, alt=vert, anomaly=anomaly)
    295296           what_I_plot = what_I_plot*mult
    296297           if not error:
    297                fvar = varname
    298                ###
    299                if anomaly:
    300                    what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,12)) - 1.)
    301                    fvar = 'anomaly'
    302                #if mult != 1:     
    303                #    fvar = str(mult) + "*" + var
    304                ###
     298               fvar = varname
     299               if anomaly: fvar = 'anomaly'
     300               ##### MAPMODE-SPECIFIC SETTINGS #####
    305301               if mapmode == 1:
    306302                   if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    … …  
    308304                   what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    309305                         indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
    310                if (fileref is not None) and (index_f is numplot-1):
    311                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
    312                else:
    313                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
     306                   zxmin, zxmax = xaxis ; zymin, zymax = yaxis
     307                   if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
     308                   if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
     309                   if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
     310                   if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
     311                   if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
     312                   if ylog:         mpl.pyplot.semilogy()
     313               if (fileref is not None) and (index_f is numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
     314               else:                                                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    314315               if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar.upper()))
    315316               elif (fileref is not None) and (index_f is numplot-1): palette = get_cmap(name="RdBu_r")
    … …  
    321322                     #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
    322323                     zelevels = np.linspace(zevmin,zevmax,num=ticks)
    323                      #contourf(what_I_plot, zelevels, cmap = palette )
    324 
    325324                     if mapmode == 1:       m.contourf( x, y, what_I_plot, zelevels, cmap = palette)
    326325                     elif mapmode == 0:     contourf( x, y, what_I_plot, zelevels, cmap = palette)
    327 
    328                      zxmin, zxmax = xaxis ; zymin, zymax = yaxis
    329                      if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
    330                      if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
    331                      if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
    332                      if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
    333                      if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
    334                      if ylog:         mpl.pyplot.semilogy()
    335 
    336326                 else:
    337327                     if hole:  what_I_plot = nolow(what_I_plot)
    338                      if mapmode == 1:
    339                         m.pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
    340                      elif mapmode == 0:
    341                         pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
    342 
     328                     if mapmode == 1:       m.pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
     329                     elif mapmode == 0:     pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
    343330                 if colorb != 'nobar' and varname != 'HGT' :       
    344331                     if (fileref is not None) and (index_f is numplot-1):
    … …  
    350337                                           ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),\
    351338                                           extend='neither',spacing='proportional')
    352 
    353339                                           # both min max neither
    354340               ##### 1D field
Note: See TracChangeset for help on using the changeset viewer.