Changeset 388


Ignore:
Timestamp:
Nov 15, 2011, 7:27:55 PM (14 years ago)
Author:
acolaitis
Message:

PYTHON.
#######
~
Added new option and routines that go along.
~
#######

M 387 meso.py
M 387 gcm.py
M 387 myscript.py
-------------------- Added the --tsat option. Allows the computation and plotting of Tsat-T instead of T when the z-axis is a pressure coordinate and the variable to plot is temperature. This is usefull to study saturation. Be carefull not to use the --column option, which uses a z-axis in meters, and would output wrong results with an axis in pressure. This could however be worked out in a future commit.

M 387 mymath.py
-------------------- Added several new functions, that may not be at their best place in mymath...:

  • get_tsat(pressure,temp=None,zlon=None,zlat=None,zalt=None,ztime=None) -->when given pressure only, the routine outputs a vertical profile of condensation temperature. -->when given all the inputs, the routine returns Tcond-T by taking care of all dimension problems, automatically finding the vertical axis in the temperature field, and taking care of eventual missing values by masking the one present in the input by the value 1e20 in the output.
  • get_dim(zlon,zlat,zalt,ztime,zvar): --> returns a dictionnary giving the position of each of the dimensions in zvar in the following form:

{'latitude': 2, 'altitude': 1, 'longitude': 3, 'Time': 0}

--> this works for any shape for zvar between 1 and 4 dimensions, and takes care of tricky cases when variable's dimensions are not in the right order, and when two dimensions have the same number of elements.

  • get_key(self,value) --> this function returns the key(s) associated to a value in a dictionnary. This is used to be able to access dictionnaries in the opposite way than the natural one.

M 387 planetoplot.py
-------------------- Added changes relative to the --tsat option

Location:
trunk/UTIL/PYTHON
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/gcm.py

    r385 r388  
    176176                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
    177177                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
    178                 blat=opt.blat)
     178                blat=opt.blat,tsat=opt.tsat)
    179179        print 'Done: '+name
    180180        system("rm -f to_be_erased")

  • trunk/UTIL/PYTHON/meso.py

    r385 r388  
    173173                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
    174174                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
    175                 blat=opt.blat)
     175                blat=opt.blat,tsat=opt.tstat)
    176176        print 'Done: '+name
    177177        system("rm -f to_be_erased")

  • trunk/UTIL/PYTHON/mymath.py

    r349 r388  
    2525    return
    2626
     27
     28# A.C. routine to compute saturation temperature
     29def get_tsat(pressure,temp=None,zlon=None,zlat=None,zalt=None,ztime=None):
     30    import math as mt
     31    import numpy as np
     32    acond=3.2403751E-04
     33    bcond=7.3383721E-03
     34    # if temp is not in input, the routine simply outputs the vertical profile
     35    # of Tsat
     36    if temp is None:
     37      # Identify dimensions in temperature field
     38      output=np.zeros(np.array(pressure).shape)
     39      if len(np.array(pressure).shape) is 1:
     40         #pressure field is a 1d column, (i.e. the altitude coordinate)
     41         #temperature has to have a z-axis
     42         i=0
     43         for pp in pressure:
     44            output[i]=1./(bcond-acond*mt.log(.0095*pp))         
     45            i=i+1
     46      else:
     47         #pressure field is a field present in the file. Unhandled
     48         #by this routine for now, which only loads unique variables.
     49         print "3D pressure field not handled for now, exiting in tsat"
     50         print "Use a vertical pressure coordinate if you want to compute Tsat"
     51         exit()
     52    # if temp is in input, the routine computes Tsat-T by detecting where the
     53    # vertical axis is in temp
     54    else:
     55      output=np.zeros(np.array(temp).shape)
     56      vardim=get_dim(zlon,zlat,zalt,ztime,temp)
     57#      m=np.ma.masked_where(temp[temp == -1],temp,copy=False)
     58#      print temp
     59#      mindex=np.ma.nonzero(m.mask)
     60#      print mindex[0]
     61      if 'altitude' not in vardim.keys():
     62         print 'no altitude coordinate in temperature field for Tsat computation'
     63         exit()
     64      zdim=vardim['altitude']
     65      ndim=len(np.array(temp).shape)
     66      print '--- in tsat(). vardim,zdim,ndim: ',vardim,zdim,ndim
     67      i=0
     68      for pp in pressure:
     69        if ndim is 1:
     70           output[i]=1./(bcond-acond*mt.log(.0095*pp))-temp[i]
     71        elif ndim is 2:
     72           if zdim is 0:
     73              output[i,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[i,:]
     74           elif zdim is 1:
     75              output[:,i]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,i]
     76           else:
     77              print "stop in get_tsat: zdim: ",zdim
     78              exit()
     79        elif ndim is 3:
     80           if zdim is 0:
     81              output[i,:,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[i,:,:]
     82           elif zdim is 1:
     83              output[:,i,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,i,:]
     84           elif zdim is 2:
     85              output[:,:,i]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,:,i]
     86           else:
     87              print "stop in get_tsat: zdim: ",zdim
     88              exit()
     89        elif ndim is 4:
     90           if zdim is 0:
     91              output[i,:,:,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[i,:,:,:]
     92           elif zdim is 1:
     93              output[:,i,:,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,i,:,:]
     94           elif zdim is 2:
     95              output[:,:,i,:]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,:,i,:]
     96           elif zdim is 3:
     97              output[:,:,:,i]=1./(bcond-acond*mt.log(.0095*pp))-temp[:,:,:,i]
     98           else:
     99              print "stop in get_tsat: zdim: ", zdim
     100              exit()
     101        else:
     102           print "stop in get_tsat: ndim: ",ndim
     103           exit()
     104        i=i+1
     105    m=np.ma.masked_where(temp == 0,temp,copy=False)
     106    zoutput=np.ma.array(output,mask=m.mask)#,fill_value=-999999)
     107    zout=zoutput.filled()
     108    return zout
     109
     110# A.C. Dirty routine to determine where are the axis of a variable
     111def get_dim(zlon,zlat,zalt,ztime,zvar):
     112   import numpy as np
     113   nx,ny,nz,nt=0,0,0,0
     114   if zlon is not None:
     115      nx=len(zlon)
     116   if zlat is not None:
     117      ny=len(zlat)
     118   if zalt is not None:
     119      nz=len(zalt)
     120   if ztime is not None:
     121      nt=len(ztime)
     122   zdims={}
     123   zdims['longitude']=nx
     124   zdims['latitude']=ny
     125   zdims['altitude']=nz
     126   zdims['Time']=nt
     127   zvardim=np.array(zvar).shape
     128   ndim=len(zvardim)
     129   zzvardim=[[]]*ndim
     130   j=0
     131   output={}
     132   for dim in zvardim:
     133       if dim not in zdims.values():
     134          print "WARNING -----------------------------"
     135          print "Dimensions given to subroutine do not match variables dimensions :"
     136          exit()
     137       else:
     138          a=get_key(zdims,dim)
     139          if len(a) is not 1:
     140             if j is 0:                ##this should solve most conflicts with Time
     141                zzvardim[j]=a[1]
     142             else:
     143                zzvardim[j]=a[0]
     144          else:
     145              zzvardim[j]=a[0]
     146          output[zzvardim[j]]=j
     147          j=j+1
     148   return output
     149
     150# A.C. routine that gets keys from a dictionnary value
     151def get_key(self, value):
     152    """find the key(s) as a list given a value"""
     153    return [item[0] for item in self.items() if item[1] == value]
     154

  • trunk/UTIL/PYTHON/myscript.py

    r385 r388  
    6060    parser.add_option('--titleref',     action='store',dest='titref',  type="string",  default="fill",  help='title for the reference file. [title of fig (1)]')
    6161
     62    ### SPECIAL
     63    parser.add_option('--tsat',         action='store_true',dest='tsat',               default=False,help='convert temperature field T in Tsat-T using pressure')
     64
    6265    return parser

  • trunk/UTIL/PYTHON/planetoplot.py

    r387 r388  
    5050           ylog=False,\
    5151           yintegral=False,\
    52            blat=None):
     52           blat=None,\
     53           tsat=False):
    5354
    5455
    … …  
    6364                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
    6465                       getname,localtime,polarinterv,getsindex,define_axis,determineplot
    65     from mymath import deg,max,min,mean
     66    from mymath import deg,max,min,mean,get_tsat
    6667    import matplotlib as mpl
    6768    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
    … …  
    104105         errormess("Not the same kind of files !", [typefile0, typefile])
    105106      ### ... VAR
     107      varname=var
    106108      if var not in nc.variables: var = False
    107109      ### ... WINDS
    … …  
    165167
    166168      print "var, var2: ", var, var2
    167       if var: all_var[k] = getfield(nc,var)
     169      if varname in ["temp","t","T_nadir_nit","T_nadir_day"] and var and tsat:
     170          print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
     171          tt=getfield(nc,var)
     172          tt.set_fill_value([0])
     173          all_var[k]=get_tsat(vert,tt,zlon=lon,zlat=lat,zalt=vert,ztime=time)
     174      else:
     175         if var: all_var[k] = getfield(nc,var)
    168176      if var2: all_var2[k] = getfield(nc,var2)
    169177     
    … …  
    320328                     #contourf(what_I_plot, zelevels, cmap = palette )
    321329
    322                      if mapmode == 0:
    323                        contourf( x , y, what_I_plot, zelevels, cmap = palette)
    324                      elif mapmode == 1:
     330                     if mapmode == 1:
    325331                       m.contourf( x, y, what_I_plot, zelevels, cmap = palette)
     332                     elif mapmode == 0:
     333                       contourf( x, y, what_I_plot, zelevels, cmap = palette)
    326334                     zxmin, zxmax = xaxis
    327335                     zymin, zymax = yaxis
    … …  
    351359                     else:
    352360                        colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
    353                                            ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
     361                                           ticks=np.linspace(zevmin,zevmax,min([ndiv+1,10])),\
    354362                                           extend='neither',spacing='proportional')
    355363
Note: See TracChangeset for help on using the changeset viewer.