Ignore:
Timestamp:
Nov 4, 2011, 2:45:06 PM (13 years ago)
Author:
aslmd
Message:

PYTHON: updates from TN. planetoplot.py is a version of winds.py which works for GCM files and adds section, 1Dplot capabilities. new functions in myplot. checked compatibility with winds.py, OK.

Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py

    r310 r345  
    55def max (field,axis=None):
    66        import numpy as np
    7         return np.array(field).max(axis=axis)
     7        if field is None: return None
     8        else: return np.array(field).max(axis=axis)
    89
    910def mean (field,axis=None):

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py

    r317 r345  
     1## Author: AS
    12def errormess(text,printvar=None):
    23    print text
    … …  
    56    return
    67
     8## Author: AS
    79def getname(var=False,winds=False,anomaly=False):
    810    if var and winds:     basename = var + '_UV'
    … …  
    1315    return basename
    1416
     17## Author: AS
    1518def localtime(utc,lon):
    1619    ltst = utc + lon / 15.
    … …  
    1922    return ltst
    2023
     24## Author: AS
    2125def whatkindfile (nc):
    2226    if 'controle' in nc.variables:   typefile = 'gcm'
    … …  
    2529    elif 'U' in nc.variables:        typefile = 'meso'
    2630    elif 'HGT_M' in nc.variables:    typefile = 'geo'
    27     else:                            errormess("whatkindfile: typefile not supported.")
     31    #else:                            errormess("whatkindfile: typefile not supported.")
     32    else: typefile = 'gcm' # for lslin-ed files from gcm
    2833    return typefile
    2934
     35## Author: AS
    3036def getfield (nc,var):
    3137    ## this allows to get much faster (than simply referring to nc.variables[var])
    3238    dimension = len(nc.variables[var].dimensions)
     39    print "   Opening variable with", dimension, "dimensions ..."
    3340    if dimension == 2:    field = nc.variables[var][:,:]
    3441    elif dimension == 3:  field = nc.variables[var][:,:,:]
    … …  
    3643    return field
    3744
     45## Author: AS + TN
    3846def reducefield (input,d4=None,d3=None,d2=None,d1=None):
    3947    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
    4048    ### it would be actually better to name d4 d3 d2 d1 as t z y x
    4149    import numpy as np
     50    from mymath import max
    4251    dimension = np.array(input).ndim
    4352    shape = np.array(input).shape
    … …  
    4554    output = input
    4655    error = False
     56    ### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
     57    if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
     58    if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
     59    if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
     60    if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
     61    ### now the main part
    4762    if dimension == 2:
    4863        if   d2 >= shape[0]: error = True
    … …  
    5267        elif d2 is not None:         output = input[d2,:]
    5368    elif dimension == 3:
    54         if   d4 >= shape[0]: error = True
    55         elif d2 >= shape[1]: error = True
    56         elif d1 >= shape[2]: error = True
     69        if   max(d4) >= shape[0]: error = True
     70        elif max(d2) >= shape[1]: error = True
     71        elif max(d1) >= shape[2]: error = True
    5772        elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
    58         elif d4 is not None and d2 is not None:         output = input[d4,d2,:]
    59         elif d4 is not None and d1 is not None:         output = input[d4,:,d1]
    60         elif d2 is not None and d1 is not None:         output = input[:,d2,d1]
    61         elif d1 is not None:                output = input[:,:,d1]
    62         elif d2 is not None:                output = input[:,d2,:]
    63         elif d4 is not None:                output = input[d4,:,:]
     73        elif d4 is not None and d2 is not None:    output=np.mean(input[d4,:,:],axis=0); output=np.mean(output[d2,:],axis=0)
     74        elif d4 is not None and d1 is not None:    output=np.mean(input[d4,:,:],0); output=np.mean(output[:,d1],1)
     75        elif d2 is not None and d1 is not None:    output=np.mean(input[:,d2,:],axis=1); output=np.mean(output[:,d1],axis=1)
     76        elif d1 is not None:                output = np.mean(input[:,:,d1],axis=2)
     77        elif d2 is not None:                output = np.mean(input[:,d2,:],axis=1)
     78        elif d4 is not None:                output = np.mean(input[d4,:,:],axis=0)
    6479    elif dimension == 4:
    65         if   d4 >= shape[0]: error = True
    66         elif d3 >= shape[1]: error = True
    67         elif d2 >= shape[2]: error = True
    68         elif d1 >= shape[3]: error = True
     80        if   max(d4) >= shape[0]: error = True
     81        elif max(d3) >= shape[1]: error = True
     82        elif max(d2) >= shape[2]: error = True
     83        elif max(d1) >= shape[3]: error = True
    6984        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:  output = input[d4,d3,d2,d1]
    7085        elif d4 is not None and d3 is not None and d2 is not None:         output = input[d4,d3,d2,:]
    … …  
    7287        elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
    7388        elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
    74         elif d4 is not None and d3 is not None:                output = input[d4,d3,:,:]
    75         elif d4 is not None and d2 is not None:                output = input[d4,:,d2,:]
    76         elif d4 is not None and d1 is not None:                output = input[d4,:,:,d1]
    77         elif d3 is not None and d2 is not None:                output = input[:,d3,d2,:]
    78         elif d3 is not None and d1 is not None:                output = input[:,d3,:,d1]
    79         elif d2 is not None and d1 is not None:                output = input[:,:,d2,d1]
     89        elif d4 is not None and d3 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[d3,:,:],0)
     90        elif d4 is not None and d2 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,d2,:],1)
     91        elif d4 is not None and d1 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,:,d1],2)
     92        elif d3 is not None and d2 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,d2,:],1)
     93        elif d3 is not None and d1 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,:,d1],0)
     94        elif d2 is not None and d1 is not None:  output = np.mean(input[:,:,d2,:],2); output = np.mean(output[:,:,d1],2)
    8095        elif d1 is not None:                       output = input[:,:,:,d1]
    8196        elif d2 is not None:                       output = input[:,:,d2,:]
    … …  
    87102    return output, error
    88103
     104## Author: AS + TN
    89105def definesubplot ( numplot, fig ):
    90106    from matplotlib.pyplot import rcParams
    91107    rcParams['font.size'] = 12. ## default (important for multiple calls)
    92     if   numplot == 4:
    93         sub = 221
    94         fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    95         rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     108    if numplot <= 0:
     109        subv = 99999
     110        subh = 99999
     111    elif numplot == 1:
     112        subv = 99999
     113        subh = 99999
    96114    elif numplot == 2:
    97         sub = 121
     115        subv = 1
     116        subh = 2
    98117        fig.subplots_adjust(wspace = 0.35)
    99118        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
    100119    elif numplot == 3:
    101         sub = 131
     120        subv = 3
     121        subh = 1
    102122        fig.subplots_adjust(wspace = 0.5)
    103123        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    104     elif numplot == 6:
    105         sub = 231
     124    elif   numplot == 4:
     125        subv = 2
     126        subh = 2
     127        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
     128        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     129    elif numplot <= 6:
     130        subv = 2
     131        subh = 3
    106132        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
    107133        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    108     elif numplot == 8:
    109         sub = 331 #241
     134    elif numplot <= 8:
     135        subv = 2
     136        subh = 4
    110137        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    111138        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    112     elif numplot == 9:
    113         sub = 331
     139    elif numplot <= 9:
     140        subv = 3
     141        subh = 3
    114142        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    115143        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    116     elif numplot == 1:
    117         sub = 99999
    118     elif numplot <= 0:
    119         sub = 99999
     144    elif numplot <= 12:
     145        subv = 3
     146        subh = 4
     147        fig.subplots_adjust(wspace = 0.1, hspace = 0.1)
     148        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     149    elif numplot <= 16:
     150        subv = 4
     151        subh = 4
     152        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
     153        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    120154    else:
    121         print "supported: 1,2,3,4,6,8,9"
     155        print "number of plot supported: 1 to 16"
    122156        exit()
    123     return sub
    124 
     157    return subv,subh
     158
     159## Author: AS
    125160def getstralt(nc,nvert):
    126161    typefile = whatkindfile(nc)
    … …  
    140175    return stralt
    141176
     177## Author: AS
    142178def getlschar ( namefile ):
    143179    from netCDF4 import Dataset
    … …  
    169205    return lschar, zehour, zehourin
    170206
     207## Author: AS
    171208def getprefix (nc):
    172209    prefix = 'LMD_MMM_'
    … …  
    175212    return prefix
    176213
     214## Author: AS
    177215def getproj (nc):
    178216    typefile = whatkindfile(nc)
    … …  
    200238    return proj   
    201239
     240## Author: AS
    202241def ptitle (name):
    203242    from matplotlib.pyplot import title
    … …  
    205244    print name
    206245
     246## Author: AS
    207247def polarinterv (lon2d,lat2d):
    208248    import numpy as np
    … …  
    212252    return [wlon,wlat]
    213253
     254## Author: AS
    214255def simplinterv (lon2d,lat2d):
    215256    import numpy as np
    216257    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
    217258
     259## Author: AS
    218260def wrfinterv (lon2d,lat2d):
    219261    nx = len(lon2d[0,:])-1
    … …  
    231273    return [wlon,wlat]
    232274
     275## Author: AS
    233276def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
    234277    import  matplotlib.pyplot as plt
    … …  
    246289    return
    247290
     291## Author: AS
    248292def dumpbdy (field,n,stag=None):
    249293    nx = len(field[0,:])-1
    … …  
    253297    return field[n:ny-n,n:nx-n]
    254298
     299## Author: AS
    255300def getcoorddef ( nc ):   
    256301    import numpy as np
    … …  
    267312    return lon2d,lat2d   
    268313
     314## Author: AS
    269315def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
    270316    import numpy as np
    … …  
    279325    return lon2d,lat2d
    280326
     327## Author: AS
    281328def smooth (field, coeff):
    282329        ## actually blur_image could work with different coeff on x and y
    … …  
    285332        return result
    286333
     334## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
    287335def gauss_kern(size, sizey=None):
    288336        import numpy as np
    289         ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth     
    290337        # Returns a normalized 2D gauss kernel array for convolutions
    291338        size = int(size)
    … …  
    298345        return g / g.sum()
    299346
     347## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
    300348def blur_image(im, n, ny=None) :
    301349        from scipy.signal import convolve
    302         ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
    303350        # blurs the image by convolving with a gaussian kernel of typical size n.
    304351        # The optional keyword argument ny allows for a different size in the y direction.
    … …  
    307354        return improc
    308355
     356## Author: AS
    309357def getwinddef (nc):   
    310358    ## getwinds for predefined types
    … …  
    322370    return uchar,vchar,metwind
    323371
     372## Author: AS
    324373def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
    325374    ## scale regle la reference du vecteur
    … …  
    344393    return
    345394
     395## Author: AS
    346396def display (name):
    347397    from os import system
    … …  
    349399    return name
    350400
     401## Author: AS
    351402def findstep (wlon):
    352403    steplon = int((wlon[1]-wlon[0])/4.)  #3
    … …  
    361412    return step
    362413
    363 def define_proj (char,wlon,wlat,back=None):
     414## Author: AS
     415def define_proj (char,wlon,wlat,back=None,blat=False):
    364416    from    mpl_toolkits.basemap            import Basemap
    365417    import  numpy                           as np
    … …  
    368420    meanlon = 0.5*(wlon[0]+wlon[1])
    369421    meanlat = 0.5*(wlat[0]+wlat[1])
    370     if   wlat[0] >= 80.:   blat =  40.
    371     elif wlat[1] <= -80.:  blat = -40.
    372     elif wlat[1] >= 0.:    blat = wlat[0]
    373     elif wlat[0] <= 0.:    blat = wlat[1]
     422    if not blat:
     423        if   wlat[0] >= 80.:   blat =  40.
     424        elif wlat[1] <= -80.:  blat = -40.
     425        elif wlat[1] >= 0.:    blat = wlat[0]
     426        elif wlat[0] <= 0.:    blat = wlat[1]
    374427    print "blat ", blat
    375428    h = 50.  ## en km
    … …  
    407460    return m
    408461
     462## Author: AS
    409463#### test temporaire
    410464def putpoints (map,plot):
    … …  
    425479    return
    426480
     481## Author: AS
    427482def calculate_bounds(field,vmin=None,vmax=None):
    428483    import numpy as np
    … …  
    448503    return zevmin, zevmax
    449504
     505## Author: AS
    450506def bounds(what_I_plot,zevmin,zevmax):
    451507    from mymath import max,min,mean
    … …  
    461517    return what_I_plot
    462518
     519## Author: AS
    463520def nolow(what_I_plot):
    464521    from mymath import max,min
    … …  
    468525    return what_I_plot
    469526
     527## Author: AS
    470528def zoomset (wlon,wlat,zoom):
    471529    dlon = abs(wlon[1]-wlon[0])/2.
    … …  
    476534    return wlon,wlat
    477535
     536## Author: AS
    478537def fmtvar (whichvar="def"):
    479538    fmtvar    =     { \
    … …  
    489548             "TAU_ICE":      "%.2f",\
    490549             "VMR_ICE":      "%.1e",\
    491              "MTOT":         "%.0f",\
     550             "MTOT":         "%.1f",\
    492551             "anomaly":      "%.1f",\
    493552             "W":            "%.1f",\
    … …  
    497556             "ALBBARE":      "%.2f",\
    498557             "TAU":          "%.1f",\
     558             #### T.N.
     559             "TEMP":         "%.0f",\
     560             "VMR_H2OICE":   "%.0f",\
     561             "VMR_H2OVAP":   "%.0f",\
     562             "TAUTES":       "%.2f",\
     563             "TAUTESAP":     "%.2f",\
     564
    499565                    }
    500566    if whichvar not in fmtvar:
    … …  
    502568    return fmtvar[whichvar]
    503569
     570## Author: AS
    504571####################################################################################################################
    505572### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
    … …  
    514581             "HFX":          "RdYlBu",\
    515582             "ICETOT":       "YlGnBu_r",\
    516              "MTOT":         "PuBu",\
     583             #"MTOT":         "PuBu",\
     584             "CCNQ":         "YlOrBr",\
     585             "CCNN":         "YlOrBr",\
     586             "TEMP":         "Jet",\
    517587             "TAU_ICE":      "Blues",\
    518588             "VMR_ICE":      "Blues",\
    … …  
    523593             "ALBBARE":      "spectral",\
    524594             "TAU":          "YlOrBr_r",\
     595             #### T.N.
     596             "MTOT":         "Jet",\
     597             "H2O_ICE_S":    "RdBu",\
     598             "VMR_H2OICE":   "PuBu",\
     599             "VMR_H2OVAP":   "PuBu",\
    525600                     }
    526601#W --> spectral ou jet
    … …  
    531606    return whichcolorb[whichone]
    532607
     608## Author: AS
    533609def definecolorvec (whichone="def"):
    534610        whichcolor =    { \
    … …  
    549625        return whichcolor[whichone]
    550626
     627## Author: AS
    551628def marsmap (whichone="vishires"):
    552629        from os import uname
    … …  
    588665#       return whichlink
    589666
     667## Author: AS
    590668def latinterv (area="Whole"):
    591669    list =    { \
    … …  
    609687    return olon,olat
    610688
     689## Author: TN
     690def separatenames (name):
     691  from numpy import concatenate
     692  # look for comas in the input name to separate different names (files, variables,etc ..)
     693  if name is None:
     694     names = None
     695  else:
     696    names = []
     697    stop = 0
     698    currentname = name
     699    while stop == 0:
     700      indexvir = currentname.find(',')
     701      if indexvir == -1:
     702        stop = 1
     703        name1 = currentname
     704      else:
     705        name1 = currentname[0:indexvir]
     706      names = concatenate((names,[name1]))
     707      currentname = currentname[indexvir+1:len(currentname)]
     708  return names
     709
     710## Author: TN [old]
     711def getopmatrix (kind,n):
     712  import numpy as np
     713  # compute matrix of operations between files
     714  # the matrix is 'number of files'-square
     715  # 1: difference (row minus column), 2: add
     716  # not 0 in diag : just plot
     717  if n == 1:
     718    opm = np.eye(1)
     719  elif kind == 'basic':
     720    opm = np.eye(n)
     721  elif kind == 'difference1': # show differences with 1st file
     722    opm = np.zeros((n,n))
     723    opm[0,:] = 1
     724    opm[0,0] = 0
     725  elif kind == 'difference2': # show differences with 1st file AND show 1st file
     726    opm = np.zeros((n,n))
     727    opm[0,:] = 1
     728  else:
     729    opm = np.eye(n)
     730  return opm
     731 
     732## Author: TN [old]
     733def checkcoherence (nfiles,nlat,nlon,ntime):
     734  if (nfiles > 1):
     735     if (nlat > 1):
     736        errormess("what you asked is not possible !")
     737  return 1
     738
     739## Author: TN
     740def readslices(saxis):
     741  from numpy import empty
     742  if saxis == None:
     743     zesaxis = None
     744  else:
     745     zesaxis = empty((len(saxis),2))
     746     for i in range(len(saxis)):
     747        a = separatenames(saxis[i])
     748        if len(a) == 1:
     749           zesaxis[i,:] = float(a[0])
     750        else:
     751           zesaxis[i,0] = float(a[0])
     752           zesaxis[i,1] = float(a[1])
     753           
     754  return zesaxis
     755
     756## Author: TN
     757def  getsindex(saxis,index,axis):
     758# input  : all the desired slices and the good index
     759# output : all indexes to be taken into account for reducing field
     760  import numpy as np
     761  if saxis is None:
     762      zeindex = None
     763  else:
     764      aaa = int(np.argmin(abs(saxis[index,0] - axis)))
     765      bbb = int(np.argmin(abs(saxis[index,1] - axis)))
     766      [imin,imax] = np.sort(np.array([aaa,bbb]))
     767      zeindex = np.array(range(imax-imin+1))+imin
     768      # because -180 and 180 are the same point in longitude,
     769      # we get rid of one for averaging purposes.
     770      if axis[imin] == -180 and axis[imax] == 180:
     771         zeindex = zeindex[0:len(zeindex)-1]
     772         print "whole longitude averaging asked, so last point is not taken into account."
     773  return zeindex
     774     
     775## Author: TN
     776def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
     777# Purpose of define_axis is to find x and y axis scales in a smart way
     778# x axis priority: 1/time 2/lon 3/lat 4/vertical
     779# To be improved !!!...
     780   from numpy import array,swapaxes
     781   x = None
     782   y = None
     783   count = 0
     784   what_I_plot = array(what_I_plot)
     785   shape = what_I_plot.shape
     786   if indextime is None:
     787      x = time
     788      count = count+1
     789   if indexlon is None:
     790      if count == 0: x = lon
     791      else: y = lon
     792      count = count+1
     793   if indexlat is None:
     794      if count == 0: x = lat
     795      else: y = lat
     796      count = count+1
     797   if indexvert is None and dim0 is 4:
     798      if vertmode == 0: # vertical axis is as is (GCM grid)
     799         if count == 0: x=range(len(vert))
     800         else: y=range(len(vert))
     801         count = count+1
     802      else: # vertical axis is in kms
     803         if count == 0: x = vert
     804         else: y = vert
     805         count = count+1
     806   x = array(x)
     807   y = array(y)
     808   if len(shape) == 1:
     809      if shape != len(x):
     810         print "WARNING HERE !!!"
     811         x = y
     812   elif len(shape) == 2:
     813      print shape[1], len(y), shape[0], len(x)
     814      if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
     815         what_I_plot = swapaxes(what_I_plot,0,1)
     816         print "swapaxes", what_I_plot.shape, shape
     817         #x0 = x
     818         #x = y
     819         #y = x0
     820   #print "define_axis", x, y
     821   return what_I_plot,x,y

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py

    r317 r345  
    8282    ### Open a figure and set subplots
    8383    fig = figure()
    84     sub = definesubplot( numplot, fig )
     84    subv,subh = definesubplot( numplot, fig )
    8585 
    8686    #################################
    … …  
    102102           if nplot > numplot: break
    103103           if numplot > 1:     
    104                if typefile not in ['geo']:  subplot(sub+nplot-1)
     104               if typefile not in ['geo']:  subplot(subv,subh,nplot)
    105105           found_lct = True
    106106       ### If only one local time is requested (numplot < 0)
Note: See TracChangeset for help on using the changeset viewer.