Ignore:
Timestamp:
Feb 11, 2013, 2:22:32 PM (12 years ago)
Author:
aslmd
Message:

UTIL PYTHON. Added a speedy mode for large file ! Not all options are available for the moment. Also better access to planetoplot by not preloading all functions (this is better for lisibility too).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/myplot.py

    r874 r876  
     1
     2import numpy as np
     3import netCDF4
     4
    15## Author: AS
    26def errormess(text,printvar=None):
    … …  
    812## Author: AS
    913def adjust_length (tab, zelen):
    10     from numpy import ones
    1114    if tab is None:
    12         outtab = ones(zelen) * -999999
     15        outtab = np.ones(zelen) * -999999
    1316    else:
    1417        if zelen != len(tab):
    1518            print "not enough or too much values... setting same values all variables"
    16             outtab = ones(zelen) * tab[0]
     19            outtab = np.ones(zelen) * tab[0]
    1720        else:
    1821            outtab = tab
    … …  
    3134## Author: AS + AC
    3235def localtime(time,lon,namefile): # lon is the mean longitude of the plot, not of the domain. central lon of domain is taken from cen_lon
    33     import numpy as np
    34     from netCDF4 import Dataset
    3536    ## THIS IS FOR MESOSCALE
    36     nc = Dataset(namefile)
     37    nc = netCDF4.Dataset(namefile)
    3738    ## get start date and intervals
    3839    dt_hour=1. ; start=0.
    … …  
    9394
    9495## Author: AS
     96def getfieldred (nc,var,indexlon,indexlat,indexvert,indextime):
     97    dimension = len(nc.variables[var].dimensions)
     98    ## this allows to get much faster and use much less memory esp. with large datasets
     99    if dimension == 2:    field = nc.variables[var][indextime,indexlon]
     100    elif dimension == 3:  field = nc.variables[var][indextime,indexlat,indexlon]
     101    elif dimension == 4:  field = nc.variables[var][indextime,indexvert,indexlat,indexlon]
     102    elif dimension == 1:  field = nc.variables[var][indextime]
     103    return field
     104
     105## Author: AS
    95106def getfield (nc,var):
     107    dimension = len(nc.variables[var].dimensions)
    96108    ## this allows to get much faster (than simply referring to nc.variables[var])
    97     import numpy as np
    98     dimension = len(nc.variables[var].dimensions)
    99     #print "   Opening variable with", dimension, "dimensions ..."
     109    print "   Opening variable",var," with", dimension, "dimensions ..."
    100110    if dimension == 2:    field = nc.variables[var][:,:]
    101111    elif dimension == 3:  field = nc.variables[var][:,:,:]
    … …  
    142152# The corresponding variable to call is UV or uvmet (to use api)
    143153def windamplitude (nc,mode):
    144     import numpy as np
    145154    varinfile = nc.variables.keys()
    146155    if "U" in varinfile: zu=getfield(nc,'U')
    … …  
    177186# this only requires co2col so that you can concat.nc at low cost
    178187def enrichment_factor(nc,lon,lat,time):
    179     import numpy as np
    180188    from myplot import reducefield
    181189    varinfile = nc.variables.keys()
    … …  
    213221# The corresponding variable to call is SLOPEXY
    214222def slopeamplitude (nc):
    215     import numpy as np
    216223    varinfile = nc.variables.keys()
    217224    if "slopex" in varinfile: zu=getfield(nc,'slopex')
    … …  
    233240# The corresponding variable to call is DELTAT
    234241def deltat0t1 (nc):
    235     import numpy as np
    236242    varinfile = nc.variables.keys()
    237243    if "tsurf" in varinfile: zu=getfield(nc,'tsurf')
    … …  
    251257    ### it would be actually better to name d4 d3 d2 d1 as t z y x
    252258    ### ... note, anomaly is only computed over d1 and d2 for the moment
    253     import numpy as np
    254259    from mymath import max,mean,min,sum,getmask
    255260    csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
    … …  
    444449    from mymath import max,mean
    445450    from scipy import integrate
    446     import numpy as np
    447451    if yint and vert is not None and indice is not None:
    448452       if type(input).__name__=='MaskedArray':
    … …  
    483487        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
    484488    elif numplot <= 6:
    485         subv = 2
    486         subh = 3
    487         #fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
    488         fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
     489        subv = 3#2
     490        subh = 2#3
     491        ##fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
     492        #fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
     493        fig.subplots_adjust(wspace = 0.3, hspace = 0.5)
    489494        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    490495    elif numplot <= 8:
    … …  
    531536## Author: AS
    532537def getlschar ( namefile, getaxis=False ):
    533     from netCDF4 import Dataset
    534538    from timestuff import sol2ls
    535     from numpy import array
    536539    from string import rstrip
    537540    import os as daos
    … …  
    542545      namefile = namefiletest
    543546      #### we assume that wrfout is next to wrfout_z and wrfout_zabg
    544       nc  = Dataset(namefile)
     547      nc  = netCDF4.Dataset(namefile)
    545548      zetime = None
    546549      days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
    … …  
    548551      if 'Times' in nc.variables:
    549552        zetime = nc.variables['Times'][0]
    550         shape = array(nc.variables['Times']).shape
     553        shape = np.array(nc.variables['Times']).shape
    551554        if shape[0] < 2: zetime = None
    552555    if zetime is not None \
    … …  
    621624## Author: AS
    622625def polarinterv (lon2d,lat2d):
    623     import numpy as np
    624626    wlon = [np.min(lon2d),np.max(lon2d)]
    625627    ind = np.array(lat2d).shape[0] / 2  ## to get a good boundlat and to get the pole
    … …  
    629631## Author: AS
    630632def simplinterv (lon2d,lat2d):
    631     import numpy as np
    632633    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
    633634
    … …  
    679680## Author: AS + AC
    680681def getcoorddef ( nc ):   
    681     import numpy as np
    682682    ## getcoord2d for predefined types
    683683    typefile = whatkindfile(nc)
    … …  
    712712## Author: AS
    713713def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
    714     import numpy as np
    715714    if nlon == "nothing" or nlat == "nothing":
    716715        print "NO LAT LON FIELDS. I AM TRYING MY BEST. I ASSUME GLOBAL FIELD."
    … …  
    740739## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
    741740def smooth1d(x,window_len=11,window='hanning'):
    742     import numpy
    743741    """smooth the data using a window with requested size.
    744742    This method is based on the convolution of a scaled window with the signal.
    … …  
    762760    TODO: the window parameter could be the window itself if an array instead of a string   
    763761    """
    764     x = numpy.array(x)
     762    x = np.array(x)
    765763    if x.ndim != 1:
    766764        raise ValueError, "smooth only accepts 1 dimension arrays."
    … …  
    771769    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
    772770        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
    773     s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
     771    s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
    774772    #print(len(s))
    775773    if window == 'flat': #moving average
    776         w=numpy.ones(window_len,'d')
     774        w=np.ones(window_len,'d')
    777775    else:
    778776        w=eval('numpy.'+window+'(window_len)')
    779     y=numpy.convolve(w/w.sum(),s,mode='valid')
     777    y=np.convolve(w/w.sum(),s,mode='valid')
    780778    return y
    781779
    … …  
    789787## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
    790788def gauss_kern(size, sizey=None):
    791         import numpy as np
    792789        # Returns a normalized 2D gauss kernel array for convolutions
    793790        size = int(size)
    … …  
    832829    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
    833830    import  matplotlib.pyplot               as plt
    834     import  numpy                           as np
    835831    #posx = np.min(x) - np.std(x) / 10.
    836832    #posy = np.min(y) - np.std(y) / 10.
    … …  
    874870def define_proj (char,wlon,wlat,back=None,blat=None,blon=None):
    875871    from    mpl_toolkits.basemap            import Basemap
    876     import  numpy                           as np
    877872    import  matplotlib                      as mpl
    878873    from mymath import max
    … …  
    978973## Author: AS
    979974def calculate_bounds(field,vmin=None,vmax=None):
    980     import numpy as np
    981975    from mymath import max,min,mean
    982976    ind = np.where(field < 9e+35)
    … …  
    10221016## Author : AC
    10231017def hole_bounds(what_I_plot,zevmin,zevmax):
    1024     import numpy as np
    10251018    zi=0
    10261019    for i in what_I_plot:
    … …  
    12641257## Author: TN
    12651258def separatenames (name):
    1266   from numpy import concatenate
    12671259  # look for comas in the input name to separate different names (files, variables,etc ..)
    12681260  if name is None:
    … …  
    12791271      else:
    12801272        name1 = currentname[0:indexvir]
    1281       names = concatenate((names,[name1]))
     1273      names = np.concatenate((names,[name1]))
    12821274      currentname = currentname[indexvir+1:len(currentname)]
    12831275  return names
    … …  
    12861278## Author: TN
    12871279def readslices(saxis):
    1288   from numpy import empty
    12891280  if saxis == None:
    12901281     zesaxis = None
    12911282  else:
    1292      zesaxis = empty((len(saxis),2))
     1283     zesaxis = np.empty((len(saxis),2))
    12931284     for i in range(len(saxis)):
    12941285        a = separatenames(saxis[i])
    … …  
    13041295def readdata(data,datatype,coord1,coord2):
    13051296## Read sparse data
    1306   from numpy import empty
    13071297  if datatype == 'txt':
    13081298     if len(data[coord1].shape) == 1:
    … …  
    13201310## Author: AS
    13211311def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
    1322    import numpy as np
    13231312   import matplotlib.pyplot as mpl
    13241313   if vlat is None:    array = (lon2d - vlon)**2
    … …  
    13421331# input  : all the desired slices and the good index
    13431332# output : all indexes to be taken into account for reducing field
    1344   import numpy as np
    13451333  if ( np.array(axis).ndim == 2):
    13461334      axis = axis[:,0]
    … …  
    13641352# x axis priority: 1/time 2/lon 3/lat 4/vertical
    13651353# To be improved !!!...
    1366    from numpy import array,swapaxes
    13671354   x = None
    13681355   y = None
    13691356   count = 0
    1370    what_I_plot = array(what_I_plot)
     1357   what_I_plot = np.array(what_I_plot)
    13711358   shape = what_I_plot.shape
    13721359   if indextime is None and len(time) > 1:
    … …  
    13941381         else: y = vert
    13951382         count = count+1
    1396    x = array(x)
    1397    y = array(y)
     1383   x = np.array(x)
     1384   y = np.array(y)
    13981385   print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
    13991386   if len(shape) == 1:
    … …  
    14021389   elif len(shape) == 2:
    14031390       if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
    1404            print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = swapaxes(what_I_plot,0,1)
     1391           print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = np.swapaxes(what_I_plot,0,1)
    14051392       else:
    14061393           if shape[0] != len(y):       print "WARNING: shape[0] != len(y). Correcting." ; what_I_plot = what_I_plot[0:len(y),:]
    … …  
    14451432               vecx=None,vecy=None,stride=2 ):
    14461433    ### an easy way to map a field over lat/lon grid
    1447     import numpy as np
    14481434    import matplotlib.pyplot as mpl
    14491435    from matplotlib.cm import get_cmap
    … …  
    14961482      specificname_gcm = ['enfact']
    14971483
     1484      zncvar = znc.variables
     1485
    14981486      ## Check for variable in file:
    14991487      if mode == 'check':     
    15001488          varname = zvarname
    1501           varinfile=znc.variables.keys()
    1502           logical_novarname = zvarname not in znc.variables
     1489          varinfile=zncvar.keys()
     1490          logical_novarname = zvarname not in zncvar
    15031491          logical_nospecificname_meso = not ((ztypefile in ['meso']) and (zvarname in specificname_meso))
    15041492          logical_nospecificname_gcm = not ((ztypefile in ['gcm']) and (zvarname in specificname_gcm))
    … …  
    15191507              all_var=get_tsat(yalt,tinput,zlon=ylon,zlat=ylat,zalt=yalt,ztime=ytime)
    15201508          ### ----------- 2. wind amplitude
    1521           elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
     1509          elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in zncvar)):
    15221510              all_var=windamplitude(znc,'amplitude')
    1523           elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
     1511          elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in zncvar)):
    15241512              plot_x, plot_y = windamplitude(znc,zvarname)
    15251513              if plot_x is not None: all_var=plot_x # dummy
    15261514              else: all_var=plot_y ; plot_x = None ; plot_y = None # Hodograph type 2 is not 'xy' mode
    1527           elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
     1515          elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in zncvar)):
    15281516              all_var=slopeamplitude(znc)
    15291517          ### ------------ 3. Near surface instability
    1530           elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
     1518          elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in zncvar)):
    15311519              all_var=deltat0t1(znc)
    15321520          ### ------------ 4. Enrichment factor
    … …  
    15341522              all_var=enrichment_factor(znc,ylon,ylat,ytime)
    15351523          ### ------------ 5. teta -> temp
    1536           elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in znc.variables.keys())):
     1524          elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in zncvar.keys())):
    15371525              all_var=teta_to_tk(znc)
    15381526          else:
    … …  
    15491537# (which is not efficient but it is still ok) and then, make the average (if the user wants to)
    15501538def spectrum(var,time,vert,lat,lon):
    1551     import numpy as np
    15521539    fft=np.fft.fft(var,axis=1)
    15531540    N=len(vert)
    … …  
    15651552# Computes temperature from potential temperature for mesoscale files, without the need to use API, i.e. using natural vertical grid
    15661553def teta_to_tk(nc):
    1567     import numpy as np
    15681554    varinfile = nc.variables.keys()
    15691555    p0=610.
    … …  
    15881574# 6/ in this slab, find the point at which the surface pressure is minimum
    15891575def find_devil(nc,indextime):
    1590     import numpy as np
    15911576    from scipy import ndimage
    15921577    from mymath import array2image,image2array
Note: See TracChangeset for help on using the changeset viewer.