Changeset 382 for trunk/UTIL/PYTHON


Ignore:
Timestamp:
Nov 14, 2011, 4:14:44 PM (13 years ago)
Author:
acolaitis
Message:

PYTHON.

Added new option: --yint, to switch the way vertical coordinates are treated when --vert z1,z2 is specified.

Without specifying --yint, a simple average of the field is done along the vertical coordinate:
mean(vert[z1,z2])

When specifying --yint, a trapezoidal integration of the field is performed along the vertical coordinate.
It is worth mentionning that this operation greatly depends on the "quality" of the z-axis data. By default, z-axis values
stored in the netcdf file for GCM is an approximation. We recommend using zrecast on the file to get trusty values of altitude.
This can be simply done by specifying -i 4 in the plot command (for gcm).

Location:
trunk/UTIL/PYTHON
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/gcm.py

    r381 r382  
    175175                outputname=opt.out,resolution=opt.res,\
    176176                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
    177                 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy)
     177                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.yint)
    178178        print 'Done: '+name
    179179        system("rm -f to_be_erased")

  • trunk/UTIL/PYTHON/meso.py

    r380 r382  
    172172                outputname=opt.out,resolution=opt.res,\
    173173                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
    174                 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy)
     174                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.yint)
    175175        print 'Done: '+name
    176176        system("rm -f to_be_erased")

  • trunk/UTIL/PYTHON/myplot.py

    r376 r382  
    5656    return field
    5757
    58 ## Author: AS + TN
    59 def reducefield (input,d4=None,d3=None,d2=None,d1=None):
     58## Author: AS + TN + AC
     59def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None):
    6060    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
    6161    ### it would be actually better to name d4 d3 d2 d1 as t z y x
    … …  
    9797        elif max(d2) >= shape[2]: error = True
    9898        elif max(d1) >= shape[3]: error = True
    99         elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: 
    100             output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
     99        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
     100             output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
    101101        elif d4 is not None and d3 is not None and d2 is not None:
    102             output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0)
     102            output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0)
    103103        elif d4 is not None and d3 is not None and d1 is not None:
    104             output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[:,d1],axis=1)
     104            output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[:,d1],axis=1)
    105105        elif d4 is not None and d2 is not None and d1 is not None:
    106106            output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
    107107        elif d3 is not None and d2 is not None and d1 is not None:
    108             output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
    109         elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0)
     108            output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
     109        elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3])
    110110        elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
    111111        elif d4 is not None and d1 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2)
    112         elif d3 is not None and d2 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1)
    113         elif d3 is not None and d1 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0)
     112        elif d3 is not None and d2 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1)
     113        elif d3 is not None and d1 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,:,d1],axis=0)
    114114        elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
    115115        elif d1 is not None:                     output = mean(input[:,:,:,d1],axis=3)
    116116        elif d2 is not None:                     output = mean(input[:,:,d2,:],axis=2)
    117         elif d3 is not None:                     output = mean(input[:,d3,:,:],axis=1)
     117        elif d3 is not None:                     output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt[d3])
    118118        elif d4 is not None:                     output = mean(input[d4,:,:,:],axis=0)
    119119    dimension = np.array(output).ndim
    … …  
    121121    print 'dim,shape: ',dimension,shape
    122122    return output, error
     123
     124## Author: AC
     125
     126def reduce_zaxis (input,ax=None,yint=False,vert=None):
     127    from mymath import max,mean
     128    from scipy import integrate
     129    if yint:
     130       output = integrate.trapz(input,x=vert,axis=ax)
     131    else:
     132       output = mean(input,axis=ax)
     133    return output
    123134
    124135## Author: AS + TN
    … …  
    578589             "ALBBARE":      "%.2f",\
    579590             "TAU":          "%.1f",\
     591             "CO2":          "%.2f",\
    580592             #### T.N.
    581593             "TEMP":         "%.0f",\
    … …  
    617629             "ALBBARE":      "spectral",\
    618630             "TAU":          "YlOrBr_r",\
     631             "CO2":          "YlOrBr_r",\
    619632             #### T.N.
    620633             "MTOT":         "Jet",\

  • trunk/UTIL/PYTHON/myscript.py

    r380 r382  
    4343    parser.add_option('--lon',          action='append',dest='slon',   type="string",  default=None, help='slices along lon. 2 comma-separated values: averaging')
    4444    parser.add_option('--vert',         action='append',dest='svert',  type="string",  default=None, help='slices along vert. 2 comma-separated values: averaging')
     45    parser.add_option('--yint',         action='store_true',dest='yint',               default=False,help='change the function of --vert z1,z2 from MEAN to INTEGRATE along z')
    4546    parser.add_option('--time',         action='append',dest='stime',  type="string",  default=None, help='slices along time. 2 comma-separated values: averaging')
    4647    parser.add_option('--xmax',         action='store',dest='xmax',    type="float",   default=None, help='max value for x-axis in contour-plots [max(xaxis)]')

  • trunk/UTIL/PYTHON/planetoplot.py

    r380 r382  
    4848           xaxis=[None,None],\
    4949           yaxis=[None,None],\
    50            ylog=False):
     50           ylog=False,\
     51           yintegral=False):
    5152
    5253
    … …  
    265266       #### Contour plot
    266267       if var2:
    267            what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
     268           what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, yint=yintegral, alt=vert)
    268269           #what_I_plot = what_I_plot*mult
    269270           if not error:
    … …  
    288289       #### Shaded plot
    289290       if var:
    290            what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
     291           what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert , yint=yintegral, alt=vert)
    291292           what_I_plot = what_I_plot*mult
    292293           if not error:
    … …  
    315316                     if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
    316317                     #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
    317                      zelevels = np.linspace(zevmin,zevmax)#,num=ndiv+1)
     318                     zelevels = np.linspace(zevmin,zevmax,num=2.*ndiv+1)
    318319                     #contourf(what_I_plot, zelevels, cmap = palette )
    319320
    … …  
    365366       ### Vector plot --- a adapter
    366367       if winds:
    367            vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert )
    368            vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert )
     368           vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
     369           vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
    369370           #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
    370371           if not error:
Note: See TracChangeset for help on using the changeset viewer.