Changeset 763


Ignore:
Timestamp:
Aug 28, 2012, 3:08:35 PM (12 years ago)
Author:
acolaitis
Message:

###################################################
# PYTHON / PLANETPLOT #
###################################################

# ------------------- XY plots ------------------ #

# Added a new category of plot to unidim, contour, etc... called "xy"
# - xy plots are plots that do not use time,vert,lat or lon as axis
# - variables to be plotted are stored in plot_x and plot_y, which is done in
# select_getfield. (there is no "what_I_plot" var for "wy" plots)
# - plot_x and plot_y are also subject to reduce field. A None value indicates
# the plot is not 'xy'
# - "xy" plots are used for a specific subset of plots : histograms, fourier
# transforms, hodographs
# - added option --analysis to perform certain kinds of analysis on the data
# (corresponding to xy plots).
# - One selects the fields he wants to plot (e.g. -v UV --lon 0 --lat 20) and
# chooses a kind of analysis :

--analysis fft # in the particular case given above, this mode corresponds

# to the mean of the ampitude spectrum of the vertical spatial
# fast fourier transform taken at each time index
# note that for now, fft are always done along the vertical
# axis. This could be made more flexible.
# not that the minimum wavelength you can plot depends on the
# vertical step of your simulation. THIS STEP HAS TO BE
# CONSTANT, hence you MUST use API or ZRECAST with constant
# spacing.

--analysis histo # histogram on the flattened data. If the user asks for --lon 0,20,

# the average is done before computing the histogram (contrary to the fft).
# However, if given a 2D array (with only --lon and --lat on a
# 4D field for example), the data is flattened before computation
# so that the result is still a 1D histogram.

--analysis histodensity # histogram with a kernel density estimate to a

# gaussian distribution, also giving the mean, variance,
# skewness and kurtosis. Other distributions are available in
# the scipy.stats package and could be implemented.

# - added variables "-v hodograph" and "-v hodograph_2". First one is a
# regular hodograph with "u" and "v" as axis, with labels of the local times
# (use --axtime lt). This is a "xy" plot (you must specify a vertical level
# as well, usually --vert 0 with an interpolation at 5m (-i 4 -l 0.005)).
# Second one is the variation with local time (for exemple) of the wind
# rotation (arctan(v/u)). This is not a "xy" plot but "unidim".
# For --ope plots in "xy" cases, only the operation plot is displayed.

# ------------------- Operations ------------------ #

# - For operations --ope -, the histogram (fourth plot) has been removed. To get it
# back, call --ope -_histo.
# - _only has been added to "+","-","-%" operations (eg "-_only")
# - For operations "-","+","-%","-_only","+_only","-%_only","-_histo", it is now
# possible to call multiple files (sill one variable) and compare each of them
# with the unique given reference file. Ex: -f file1,file2 --fref file3 --ope - will give 6 plots:

file1 file3 file1-file3
file2 file3 file2-file3

# In the case of "xy" plots, the multiple operation plots are regrouped on a
# single plot (using multiple lines). the title of this plot is not 'fig(2)-fig(1)'
# (default) but the argument of the --title command. Labels have to be given
# as following : --labels "dummy","dummy","file1-file3 label","dummy","dummy","file2-file3 label" (dummy can be anything. this is to be improved)
# To be able to run on multiple files and easily introduce the correct number and order of plots, these operations have been moved outside of the main loop
# on namefiles and variables. Operation of the type "add_var" or "cat" have
# been left inside the loop and are unchanged.

# ------------------- Localtime ------------------ #

# Changed the way localtime is computed. Reasons:
# - it was assuming one timestep per hour in computations mixing indices and
# actual times
# - to determine the starting date, it was using the name of the file (for Meso), instead of using the netcdf
# attribute (START_DATE)
# - it was using the computed mean longitude of the domain, which is not correct for
# hemispheric domain. => better to use the netcdf attribute CEN_LON
# - it was using this mean longitude even for profile plots at given
# longitude => better to get the local time at this given lon, especially
# important for large domains
# => new localtime() is in myplot (old one is commented). Interv is obsolete (but not removed yet). Case "Geo" has not been looked at.
# new localtime in myplot correctly account for starting date of the file in
# all cases
# accounts for local longitude of the plot
# accounts for files that do not have per hour outputs, but per timestep.
# specific cases can be added in myplot in localtime()

# ------------------ Misc ------------------ #

# - added option --xlog to get x logarithmic axis (--ylog already existed)
# - added the possibility to use 2 files of different gridding (-f
# file1,file2) although of the same type (meso for exemple). For that purpose,
# lon, lat, alt and vert arrays are now indexed with 'index_f', as for
# all_var. => all_lon, all_lat, all_lat, all_vert
# - added function teta_to_tk in myplot so that a call to pp.py on standard
# LMD mesoscale file with "-v tk" can be done without the need to call API.
# The temperature is computed from T and PTOT, knowing P0, T0 and R_CP.
# - bug correction in determineplot() that was causing wrong plot number and
# slices number when not calling averaged lon, lat, vert or times. (which is
# often !)
# - added new options to redope: --redope edge_x1, edge_x2, edge_y1, edge_y2
# which plots the boundaries of the domain. This is different compared to
# asking for a fixed longitude, because the domain boundary might not be at constant
# longitude (hemispheric domain for exemple). x1 is the western boundary, x2
# the eastern, y1 the southern and y2 the northern. x1 and x2 reduce the
# dimension along --lon, y1 and y2 reduce the dimension along --lat.
# - added control in windamplitude() to determine whether winds are staggered or
# unstaggered. This is usefull when dealing with non LMD_MMM files.
# - corrected a bug in reduce_field where the mean was computed on the wrong
# axis !!! (pretty serious bug)

# Exemple of plots you can do with these new options can be found in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/Plots_MasterScript/bam.sh

# ------------------ API ------------------ #

# - changed maximum of levels from 299 to 1000 in API (interpolation on 1000 levels is
# usefull to get larger bandwidth in fourier transform)
# the following concerns users of MRAMS files.
# - API has not been modified for MRAMS files. Instead, a python script
# (ic.py) is run on MRAMS .ctl and .dat files, which automatically format those files
# to be API and pp.py compatible.
# - ic.py is in $YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/File_conversion

###################################################
# INTERCOMPARISON TOOLS: #
###################################################

#ic.py in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/File_conversion
#CDO installer with import_binary in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/CDO
#Plotting scripts in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/Plots_MasterScript
#1D sensibility tool in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/1D_sensibility_tool

# See README in each of these folders for details.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90

    r755 r763  
    4444      CHARACTER (LEN=20)                                 :: process
    4545      CHARACTER (LEN=2000)                               :: fields
    46       REAL, DIMENSION(299)                               :: interp_levels
     46      REAL, DIMENSION(1000)                               :: interp_levels
    4747      INTEGER                                            :: interp_method=1
    4848      INTEGER                                            :: extrapolate=0
    … …  
    154154      CHARACTER (LEN=20)                                 :: process, dummy
    155155      CHARACTER (LEN=2000)                               :: fields, process_these_fields
    156       REAL, DIMENSION(299)                               :: interp_levels
     156      REAL, DIMENSION(1000)                               :: interp_levels
    157157      REAL                                               :: rval
    158158      REAL                                               :: MISSING=1.e36

  • trunk/UTIL/PYTHON/myplot.py

    r760 r763  
    2929    return basename
    3030
    31 ## Author: AS
    32 def localtime(utc,lon):
    33     ltst = utc + lon / 15.
     31## Author: AS + AC
     32def localtime(time,lon,nc): # 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    dt_hour=1.
     35    start=0.
     36    if lon is not None:
     37       if lon[0,1]!=lon[0,0]: mean_lon_plot = 0.5*(lon[0,1]-lon[0,0])
     38       else: mean_lon_plot=lon[0,0]
     39    elif hasattr(nc, 'CEN_LON'): mean_lon_plot=getattr(nc, 'CEN_LON')
     40    else: mean_lon_plot=0.
     41    if hasattr(nc,'TITLE'):
     42       title=getattr(nc, 'TITLE')
     43       if hasattr(nc,'DT') and 'MRAMS' in title: dt_hour=getattr(nc, 'DT')/60. #LMD MMM is 1 output/hour (and not 1 output/timestep)
     44                        # MRAMS is 1 output/timestep, unless stride is added in ic.py
     45       if hasattr(nc,'START_DATE'):
     46          start_date=getattr(nc, 'START_DATE')
     47          start_hour=np.float(start_date[11:13])
     48          start_minute=np.float(start_date[14:16])/60.
     49          start=start_hour+start_minute # start is the local time of simu at longitude 0
     50          ltst = start + time*dt_hour + mean_lon_plot / 15.
     51       else: ltst = time*dt_hour + mean_lon_plot / 15.
     52    else: ltst = time*dt_hour + mean_lon_plot / 15.
    3453    ltst = int (ltst * 10) / 10.
    3554    ltst = ltst % 24
    … …  
    3958def check_localtime(time):
    4059    a=-1
     60    print time
    4161    for i in range(len(time)-1):
    4262       if (time[i] > time[i+1]): a=i
    … …  
    102122#          if (True in np.array(mask)):out=masked
    103123#          else:out=field
    104        else:out=field
     124       else:
     125#       # missing values from MRAMS files are 0.100E+32
     126          masked=np.ma.masked_where(field > 1e30,field)
     127          masked.set_fill_value([np.NaN])
     128          mask = np.ma.getmask(masked)
     129          if (True in np.array(mask)):out=masked
     130          else:out=field
     131#       else:out=field
    105132    return out
    106133
    107134## Author: AC
    108 # Compute the norm of the winds
     135# Compute the norm of the winds or return an hodograph
    109136# The corresponding variable to call is UV or uvmet (to use api)
    110 def windamplitude (nc):
     137def windamplitude (nc,mode):
    111138    import numpy as np
    112139    varinfile = nc.variables.keys()
    113140    if "U" in varinfile: zu=getfield(nc,'U')
    114141    elif "Um" in varinfile: zu=getfield(nc,'Um')
     142    else: errormess("you need slopex or U or Um in your file.")
    115143    if "V" in varinfile: zv=getfield(nc,'V')
    116144    elif "Vm" in varinfile: zv=getfield(nc,'Vm')
     145    else: errormess("you need V or Vm in your file.")
    117146    znt,znz,zny,znx = np.array(zu).shape
    118     if "U" in varinfile:znx=znx-1
     147    if hasattr(nc,'WEST-EAST_PATCH_END_UNSTAG'):znx=getattr(nc, 'WEST-EAST_PATCH_END_UNSTAG')
    119148    zuint = np.zeros([znt,znz,zny,znx])
    120149    zvint = np.zeros([znt,znz,zny,znx])
    121150    if "U" in varinfile:
    122        for xx in np.arange(znx):      zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
    123        for yy in np.arange(zny):      zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
     151       if hasattr(nc,'SOUTH-NORTH_PATCH_END_STAG'): zny_stag=getattr(nc, 'SOUTH-NORTH_PATCH_END_STAG')
     152       if hasattr(nc,'WEST-EAST_PATCH_END_STAG'): znx_stag=getattr(nc, 'WEST-EAST_PATCH_END_STAG')
     153       if zny_stag == zny: zvint=zv
     154       else:
     155          for yy in np.arange(zny):      zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
     156       if znx_stag == znx: zuint=zu
     157       else:
     158          for xx in np.arange(znx):      zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
    124159    else:
    125160       zuint=zu
    126161       zvint=zv
    127     return np.sqrt(zuint**2 + zvint**2)
     162    if mode=='amplitude': return np.sqrt(zuint**2 + zvint**2)
     163    if mode=='hodograph': return zuint,zvint
     164    if mode=='hodograph_2': return None, 360.*np.arctan(zvint/zuint)/(2.*np.pi)
    128165
    129166## Author: AC
    … …  
    214251       if   redope == "mint":     input = min(input,axis=0) ; d1 = None
    215252       elif redope == "maxt":     input = max(input,axis=0) ; d1 = None
     253       elif redope == "edge_y1":  input = input[:,:,0,:]    ; d2 = None
     254       elif redope == "edge_y2":  input = input[:,:,-1,:]   ; d2 = None
     255       elif redope == "edge_x1":  input = input[:,:,:,0]    ; d1 = None
     256       elif redope == "edge_x2":  input = input[:,:,:,-1]   ; d1 = None
    216257       else:                      errormess("not supported. but try lines in reducefield beforehand.")
    217258       #elif redope == "minz":     input = min(input,axis=1) ; d2 = None
    … …  
    372413               totalarea = mean(totalarea[:,:,d2,:],axis=2); totalarea = mean(totalarea[:,:,d1],axis=1)
    373414             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
    374              output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=1)/totalarea
     415             output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)/totalarea
    375416        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
    376417        elif d2 is not None:
    … …  
    12721313     
    12731314## Author: TN
    1274 def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
     1315def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode,redope):
    12751316# Purpose of define_axis is to find x and y axis scales in a smart way
    12761317# x axis priority: 1/time 2/lon 3/lat 4/vertical
    … …  
    12861327      x = time
    12871328      count = count+1
    1288    if indexlon is None and len(lon) > 1:
     1329   if indexlon is None and len(lon) > 1 and redope not in ['edge_x1','edge_x2']:
    12891330      print "AXIS is lon"
    12901331      if count == 0: x = lon
    12911332      else: y = lon
    12921333      count = count+1
    1293    if indexlat is None and len(lat) > 1:
     1334   if indexlat is None and len(lat) > 1 and redope not in ['edge_y1','edge_y2']:
    12941335      print "AXIS is lat"
    12951336      if count == 0: x = lat
    … …  
    13221363   return what_I_plot,x,y
    13231364
    1324 # Author: TN + AS
    1325 def determineplot(slon, slat, svert, stime):
     1365# Author: TN + AS + AC
     1366def determineplot(slon, slat, svert, stime, redope):
    13261367    nlon = 1 # number of longitudinal slices -- 1 is None
    13271368    nlat = 1
    … …  
    13301371    nslices = 1
    13311372    if slon is not None:
    1332         nslices = nslices*len(slon)
     1373        if slon[0,0]!=slon[0,1]: length=len(slon)
     1374        else: length=1
     1375        nslices = nslices*length
    13331376        nlon = len(slon)
    13341377    if slat is not None:
    1335         nslices = nslices*len(slat)
     1378        if slat[0,0]!=slat[0,1]: length=len(slat)
     1379        else: length=1
     1380        nslices = nslices*length
    13361381        nlat = len(slat)
    13371382    if svert is not None:
    1338         nslices = nslices*len(svert)
     1383        if svert[0,0]!=svert[0,1]: lenth=len(svert)
     1384        else: length=1
     1385        nslices = nslices*length
    13391386        nvert = len(svert)
    13401387    if stime is not None:
    1341         nslices = nslices*len(stime)
     1388        if stime[0,0]!=stime[0,1]: length=len(stime)
     1389        else: length=1
     1390        nslices = nslices*length
    13421391        ntime = len(stime)
    13431392    #else:
    13441393    #    nslices = 2 
    13451394    mapmode = 0
    1346     if slon is None and slat is None:
     1395    if slon is None and slat is None and redope not in ['edge_x1','edge_x2','edge_y1','edge_y2']:
    13471396       mapmode = 1 # in this case we plot a map, with the given projection
    1348 
    13491397    return nlon, nlat, nvert, ntime, mapmode, nslices
    1350 
    1351 ## Author: AC
    1352 ## Reduce complexity of main script by moving the contour part here. Also allow to call it from elsewhere
    1353 ## which can be usefull
    1354 #
    1355 #def call_contour(what_I_plot,error,x,y,m,lon,lat,vert,time,vertmode,ze_var2,indextime,indexlon,indexlat,indexvert,yintegral,mapmode,typefile,var2,ticks):
    1356 #      from matplotlib.pyplot import contour, plot, clabel
    1357 #      import numpy as np
    1358 #      #what_I_plot = what_I_plot*mult
    1359 #      if not error:
    1360 #         if mapmode == 1:
    1361 #             if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    1362 #         zevmin, zevmax = calculate_bounds(what_I_plot)
    1363 #         zelevels = np.linspace(zevmin,zevmax,ticks) #20)
    1364 #         if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
    1365 #         if mapmode == 0:
    1366 #             #if typefile in ['mesoideal']:    what_I_plot = dumpbdy(what_I_plot,0,stag='W')
    1367 #             itime=indextime
    1368 #             if len(what_I_plot.shape) is 3: itime=[0]
    1369 #             what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    1370 #                   itime,what_I_plot, len(ze_var2.shape),vertmode)
    1371 #         ### If we plot a 2-D field
    1372 #         if len(what_I_plot.shape) is 2:
    1373 #             #zelevels=[1.]
    1374 #             if mapmode == 0:cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
    1375 #             elif mapmode == 1:cs = m.contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
    1376 #             #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7)
    1377 #         ### If we plot a 1-D field
    1378 #         elif len(what_I_plot.shape) is 1:
    1379 #             plot(what_I_plot,x)
    1380 #      else:
    1381 #         errormess("There is an error in reducing field !")
    1382 #      return error
    13831398
    13841399## Author : AS
    … …  
    14281443## Author : AC
    14291444## Handles calls to specific computations (e.g. wind norm, enrichment factor...)
    1430 def select_getfield(zvarname=None,znc=None,ztypefile=None,mode=None,ztsat=None,ylon=None,ylat=None,yalt=None,ytime=None):
     1445def select_getfield(zvarname=None,znc=None,ztypefile=None,mode=None,ztsat=None,ylon=None,ylat=None,yalt=None,ytime=None,analysis=None):
    14311446      from mymath import get_tsat
    14321447 
    14331448      ## Specific variables are described here:
    14341449      # for the mesoscale:
    1435       specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']
     1450      specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT','hodograph','tk','hodograph_2']
    14361451      # for the gcm:
    14371452      specificname_gcm = ['enfact']
    … …  
    14521467      ## Get the corresponding variable:
    14531468      if mode == 'getvar':
     1469          plot_x = None ; plot_y = None ;
    14541470          ### ----------- 1. saturation temperature
    14551471          if zvarname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and ztsat:
    … …  
    14601476          ### ----------- 2. wind amplitude
    14611477          elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
    1462               all_var=windamplitude(znc)
     1478              all_var=windamplitude(znc,'amplitude')
     1479          elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
     1480              plot_x, plot_y = windamplitude(znc,zvarname)
     1481              if plot_x is not None: all_var=plot_x # dummy
     1482              else: all_var=plot_y ; plot_x = None ; plot_y = None # Hodograph type 2 is not 'xy' mode
    14631483          elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
    14641484              all_var=slopeamplitude(znc)
    … …  
    14691489          elif ((ztypefile in ['gcm']) and (zvarname in ['enfact'])):
    14701490              all_var=enrichment_factor(znc,ylon,ylat,ytime)
     1491          ### ------------ 5. teta -> temp
     1492          elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in znc.variables.keys())):
     1493              all_var=teta_to_tk(znc)
    14711494          else:
    14721495          ### -----------  999. Normal case
    14731496              all_var = getfield(znc,zvarname)
    1474           return all_var
    1475        
     1497          if analysis is not None:
     1498             if analysis in ['histo','density','histodensity']: plot_y=all_var ; plot_x = plot_y
     1499             elif analysis == 'fft': plot_y, plot_x = spectrum(all_var,ytime,yalt,ylat,ylon) ; all_var = plot_y
     1500          return all_var, plot_x, plot_y
     1501
     1502# Author : A.C
     1503# FFT is computed before reducefield voluntarily, because we dont want to compute
     1504# ffts on averaged fields (which would kill all waves). Instead, we take the fft everywhere
     1505# (which is not efficient but it is still ok) and then, make the average (if the user wants to)
     1506def spectrum(var,time,vert,lat,lon):
     1507    import numpy as np
     1508    fft=np.fft.fft(var,axis=1)
     1509    N=len(vert)
     1510    step=(vert[1]-vert[0])*1000.
     1511    print "step is: ",step
     1512    fftfreq=np.fft.fftfreq(N,d=step)
     1513    fftfreq=np.fft.fftshift(fftfreq) # spatial FFT => this is the wavenumber
     1514    fft=np.fft.fftshift(fft)
     1515    fftfreq = 1./fftfreq # => wavelength (div by 0 expected, don't panic)
     1516    fft=np.abs(fft) # => amplitude spectrum
     1517#    fft=np.abs(fft)**2 # => power spectrum
     1518    return fft,fftfreq
     1519
     1520# Author : A.C.
     1521# Computes temperature from potential temperature for mesoscale files, without the need to use API, i.e. using natural vertical grid
     1522def teta_to_tk(nc):
     1523    import numpy as np
     1524    varinfile = nc.variables.keys()
     1525    p0=610.
     1526    t0=220.
     1527    r_cp=1./3.89419
     1528    if "T" in varinfile: zteta=getfield(nc,'T')
     1529    else: errormess("you need T in your file.")
     1530    if "PTOT" in varinfile: zptot=getfield(nc,'PTOT')
     1531    else: errormess("you need PTOT in your file.")
     1532    zt=(zteta+220.)*(zptot/p0)**(r_cp)
     1533    return zt

  • trunk/UTIL/PYTHON/myscript.py

    r760 r763  
    6464    parser.add_option('--ymin',         action='store',dest='ymin',    type="float",   default=None, help='min value for y-axis in contour-plots [min(yaxis)]')
    6565    parser.add_option('--inverty',      action='store_true',dest='inverty',            default=False,help='force decreasing values along y-axis (e.g. p-levels)')
    66     parser.add_option('--logy',         action='store_true',dest='logy',               default=False,help='set y-axis to logarithmic')   
     66    parser.add_option('--logx',         action='store_true',dest='logx',               default=False,help='set x-axis to logarithmic')
     67    parser.add_option('--logy',         action='store_true',dest='logy',               default=False,help='set y-axis to logarithmic')
    6768    parser.add_option('--axtime',       action='store',dest='axtime',  type="string",  default=None, help='choose "ls","sol","lt" for time ref (1D or --time)')
    6869
    6970    ### OPERATIONS BETWEEN FILES
    70     parser.add_option('--operation',    action='store',dest='operat',  type="string",  default=None,  help='operation to perform on input files given through -f. "+" or "-" acts on each input file by adding or substracting the ref file specified through --fref. "cat" acts on all input files in-a-row. "add_var" "sub_var" "mul_var" "div_var" acts on two variables (add _only to get only operation plot).')
     71    parser.add_option('--operation',    action='store',dest='operat',  type="string",  default=None,  help='operation to perform on input files given through -f. "+" or "-" acts on each input file by adding or substracting the ref file specified through --fref. "cat" acts on all input files in-a-row. "add_var" "sub_var" "mul_var" "div_var" acts on two variables (add _only to get only operation plot). "-_histo" will add an histogram plot for the "-" operation.')
    7172    parser.add_option('--fref',         action='store',dest='fref',    type="string",  default=None,  help='reference namefile for the --operation option.')
    7273    parser.add_option('--mope',         action='store',dest='vminope', type="float",   default=0.,  help='bounding minimum value for inter-file operation')
    … …  
    7778    parser.add_option('--tsat',         action='store_true',dest='tsat',               default=False,help='convert temperature field T in Tsat-T using pressure')
    7879    parser.add_option('--stream',       action='store_true',dest='stream',             default=False,help='plot streamlines from streamfunction.e for vert/lat slices.')
     80    parser.add_option('--analysis',     action='store'     ,dest='analysis',           default=None ,help='Analysis operation. histo, density (kernel distribution estimate, with gaussian kernel only for the moment (many other distributions are available and can be added)), histodensity (overplot of both density and histo), fft. (currently fft works only on the z-axis, i.e. spatial Fourier Transform, and yields spectrum amplitude. To get enough bandwith, use API with a step of 10m on your data. Note that if you use --time A,B, the result will be the mean of FT at each timestep, and not the FT of the mean. The same apply to --lon and --lat, but does not apply to histo and density (for which arrays are flattened -> no mean).) [None]')
    7981
    8082    return parser

  • trunk/UTIL/PYTHON/planetoplot.py

    r760 r763  
    2626           winds=False,\
    2727           addchar=None,\
    28            interv=[0,1],\
    2928           vmin=None,\
    3029           vmax=None,\
    … …  
    5453           yaxis=[None,None],\
    5554           ylog=False,\
     55           xlog=False,\
    5656           yintegral=False,\
    5757           blat=None,\
    … …  
    7373           facwind=1.,\
    7474           trycol=False,\
    75            streamflag=False):
     75           streamflag=False,\
     76           analysis=None):
    7677
    7778    ####################################################################################################################
    … …  
    9596    import numpy as np
    9697    from numpy.core.defchararray import find
     98    from scipy.stats import gaussian_kde,describe
    9799    from videosink import VideoSink
    98100    from times import sol2ls
    … …  
    126128### we prepare number of plot fields [zelen] and number of plot instances [numplot] according to user choices
    127129### --> we support multifile and multivar plots : files + vars separated by commas are on the same figure
    128     nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
     130    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime, redope)
    129131    zelen = len(namefiles)*len(var)
     132### we correct number of plot fields for possible operation (substract, etc...)
     133    if ope is not None:
     134        if fileref is not None:       zelen = 3*len(var)*len(namefiles)
     135        elif "var" in ope:            zelen = zelen + 1
    130136    numplot = zelen*nslices
    131137    print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
    132138    print "********** MAPMODE: ", mapmode
    133 ### we correct number of plot fields for possible operation (substract, etc...)
    134     if ope is not None:
    135         if fileref is not None:       zelen = zelen + 2
    136         elif "var" in ope:            zelen = zelen + 1
    137139### we define the arrays for plot fields -- which will be placed within multiplot loops
    138     all_var  = [[]]*zelen ; all_var2  = [[]]*zelen
    139     all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen ; all_colorb = [[]]*zelen
     140    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen
     141    all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_colorb = [[]]*zelen
     142    all_time = [[]]*zelen ; all_vert = [[]]*zelen ; all_lat = [[]]*zelen ; all_lon = [[]]*zelen
    140143    all_windu = [[]]*zelen ; all_windv = [[]]*zelen
    141 
     144    plot_x = [[]]*zelen ; plot_y = [[]]*zelen ; multiplot = [[]]*zelen
     145### tool (should be move to mymath)
     146    getVar = lambda searchList, ind: [searchList[i] for i in ind]
    142147#############################
    143148### LOOP OVER PLOT FIELDS ###
    144149#############################
     150   
    145151    for nnn in range(len(namefiles)):
    146152     for vvv in range(len(var)):
    … …  
    290296              ftime = np.zeros(len(time))
    291297              for i in range(len(time)):
    292                  ftime[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
     298                 ftime[i] = localtime ( time[i], slon , nc )
    293299              time=ftime
    294300              time=check_localtime(time)
    … …  
    312318      all_namefile[k] = namefile
    313319      all_time[k] = time
     320      all_vert[k] = vert
     321      all_lat[k] = lat
     322      all_lon[k] =  lon
    314323      all_colorb[k] = clb[vvv]
    315       if var2:  all_var2[k] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
     324      if var2:  all_var2[k], plot_x[k], plot_y[k] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
    316325      if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
    317326    ### we fill the arrays of fields to be plotted at the current step considered
    318       all_var[k] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
    319     ### we inform the user about the loop then increment the loop. this is the last line of "for namefile in namefiles"
     327      all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
     328
     329      # last line of for namefile in namefiles
    320330      print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
    321331
    … …  
    323333
    324334    ##################################
    325     ### Operation on files
    326     if ope is not None:
    327         print "********** OPERATION: ",ope
    328         if "var" not in ope:
    329              if len(var) > 1: errormess("for this operation... please set only one var !")
    330              if ope in ["-","+","-%"]:
    331                 if fileref is not None:   
    332                    ncref = Dataset(fileref)
    333                    all_var[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
    334                    all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1]
    335                    if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
    336                 else: errormess("fileref is missing!")
    337                 if ope == "-":     all_var[k+1]= all_var[k-1] - all_var[k]
    338                 elif ope == "+":   all_var[k+1]= all_var[k-1] + all_var[k]
    339                 elif ope == "-%":
    340                     masked = np.ma.masked_where(all_var[k] == 0,all_var[k])
    341                     masked.set_fill_value([np.NaN])
    342                     all_var[k+1]= 100.*(all_var[k-1] - masked)/masked
    343                 all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; all_var2[k+1] = all_var2[k] ; numplot = numplot+2
    344                 if len(clb) >= zelen:  all_colorb[k+1] = clb[-1]   # last additional user-defined color is for operation plot
    345                 else:                  all_colorb[k+1] = "RdBu_r"  # if no additional user-defined color... set a good default one
    346                 if winds: all_windu[k+1] = all_windu[k-1]-all_windu[k] ; all_windv[k+1] = all_windv[k-1] - all_windv[k]
    347              elif ope in ["cat"]:
    348                 tabtime = all_time[0];tab = all_var[0];k = 1
    349                 if var2: tab2 = all_var2[0]
    350                 while k != len(namefiles) and len(all_time[k]) != 0:
    351                     if var2: tab2 = np.append(tab2,all_var2[k],axis=0)
    352                     tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
    353                 all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
    354                 if var2: all_var2[0] = np.array(tab2)
    355              else: errormess(ope+" : non-implemented operation. Check pp.py --help")
    356         else:
    357              if len(namefiles) > 1: errormess("for this operation... please set only one file !") 
    358              if len(var) > 2:       errormess("not sure this works for more than 2 vars... please check.")
    359              if   "div_var" in ope: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
    360              elif "mul_var" in ope: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
    361              elif "add_var" in ope: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
    362              elif "sub_var" in ope: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
    363              else:                    errormess(ope+" : non-implemented operation. Check pp.py --help")
    364              numplot = numplot + 1 ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
    365              all_varname[k] = all_varname[k-2] + insert + all_varname[k-1]
    366              if len(clb) >= zelen: all_colorb[k] = clb[-1]   # last additional user-defined color is for operation plot
    367              else:                 all_colorb[k] = all_colorb[k-1]  # if no additional user-defined color... set same as var
    368              ### only the operation plot. do not mention colorb so that it is user-defined?
    369              if "only" in ope:
    370                  numplot = 1 ; all_var[0] = all_var[k]
    371                  all_time[0] = all_time[k] ; all_namefile[0] = all_namefile[k]
    372                  all_varname[0] = all_varname[k-2] + insert + all_varname[k-1]
    373 
     335    ### Operation on files (I) with _var
     336    if ope is not None and "var" in ope:
     337         print "********** OPERATION: ",ope
     338         if len(namefiles) > 1: errormess("for this operation... please set only one file !")
     339         if len(var) > 2:       errormess("not sure this works for more than 2 vars... please check.")
     340         if   "div_var" in ope: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
     341         elif "mul_var" in ope: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
     342         elif "add_var" in ope: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
     343         elif "sub_var" in ope: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
     344         else:                    errormess(ope+" : non-implemented operation. Check pp.py --help")
     345         plot_x[k] = None ; plot_y[k] = None
     346         all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1]
     347         all_varname[k] = all_varname[k-2] + insert + all_varname[k-1]
     348         if len(clb) >= zelen: all_colorb[k] = clb[-1]   # last additional user-defined color is for operation plot
     349         else:                 all_colorb[k] = all_colorb[k-1]  # if no additional user-defined color... set same as var
     350         ### only the operation plot. do not mention colorb so that it is user-defined?
     351         if "only" in ope:
     352             numplot = 1 ; all_var[0] = all_var[k]
     353             all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k]
     354             all_varname[0] = all_varname[k-2] + insert + all_varname[k-1]
     355
     356
     357    ##################################
     358    ### Operation on files (II) without _var
     359    # we re-iterate on the plots to set operation subplots to make it compatible with multifile (using the same ref file)
     360    # (k+1)%3==0 is the index of operation plots
     361    # (k+2)%3==0 is the index of reference plots
     362    # (k+3)%3==0 is the index of first plots
     363    opefirstpass=True
     364    if ope is not None and "var" not in ope:
     365       print "********** OPERATION: ",ope
     366       for k in np.arange(zelen):
     367               if len(var) > 1: errormess("for this operation... please set only one var !")
     368               if ope in ["-","+","-%","-_only","+_only","-%_only","-_histo"]:
     369                  if fileref is None: errormess("fileref is missing!")
     370                  else:ncref = Dataset(fileref)
     371
     372                  if opefirstpass: ## first plots
     373                     for ll in np.arange(len(namefiles)):
     374                        print "SETTING FIRST PLOT"
     375                        all_varname[3*ll] = all_varname[ll] ; all_time[3*ll] = all_time[ll] ; all_vert[3*ll] = all_vert[ll] ; all_lat[3*ll] = all_lat[ll] ; all_lon[3*ll] = all_lon[ll] ; all_namefile[3*ll] = all_namefile[ll] ; all_var2[3*ll] = all_var2[ll] ; all_colorb[3*ll] = all_colorb[ll] ; all_var[3*ll] = all_var[ll]
     376                        if plot_y[ll] is not None: plot_y[3*ll] = plot_y[ll] ; plot_x[3*ll] = plot_x[ll]
     377                        else: plot_y[3*ll] = None ; plot_x[3*ll] = None
     378                        opefirstpass=False
     379
     380                  if (k+2)%3==0: ## reference plots
     381                        print "SETTING REFERENCE PLOT"
     382                        all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1]
     383                        all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
     384                        if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
     385
     386                  if (k+1)%3==0: ## operation plots
     387                     print "SETTING OPERATION PLOT"
     388                     all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1]
     389                     if ope in ["-","-_only","-_histo"]:
     390                         all_var[k]= all_var[k-2] - all_var[k-1]
     391                         if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] - plot_y[k-1]
     392                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
     393                     elif ope in ["+","+_only"]:   
     394                         all_var[k]= all_var[k-2] + all_var[k-1]
     395                         if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] + plot_y[k-1]
     396                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
     397                     elif ope in ["-%","-%_only"]:
     398                         masked = np.ma.masked_where(all_var[k-1] == 0,all_var[k-1])
     399                         masked.set_fill_value([np.NaN])
     400                         all_var[k]= 100.*(all_var[k-2] - masked)/masked
     401                         if plot_y[k-1] is not None and plot_y[k-2] is not None:
     402                            masked = np.ma.masked_where(plot_y[k-1] == 0,plot_y[k-1])
     403                            masked.set_fill_value([np.NaN])
     404                            plot_y[k]= 100.*(plot_y[k-2] - masked)/masked
     405                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
     406                     if len(clb) >= zelen: all_colorb[k] = clb[-1]
     407                     else: all_colorb[k] = "RdBu_r" # if no additional user-defined color... set a good default one
     408                     if winds: all_windu[k] = all_windu[k-2]-all_windu[k-1] ; all_windv[k] = all_windv[k-2] - all_windv[k-1]
     409
     410               elif ope in ["cat"]:
     411                  tabtime = all_time[0];tab = all_var[0];k = 1
     412                  if var2: tab2 = all_var2[0]
     413                  while k != len(namefiles) and len(all_time[k]) != 0:
     414                      if var2: tab2 = np.append(tab2,all_var2[k],axis=0)
     415                      tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
     416                  all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
     417                  if var2: all_var2[0] = np.array(tab2)
     418               else: errormess(ope+" : non-implemented operation. Check pp.py --help")
     419       if "only" in ope:
     420           numplot = 1 ; all_var[0] = all_var[k]
     421           all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k] ; plot_x[0]=plot_x[k] ; plot_y[0]=plot_y[k]
     422           all_varname[0] = all_varname[k]
    374423
    375424    ##################################
    … …  
    377426    fig = figure()
    378427    subv,subh = definesubplot( numplot, fig )
    379     if ope in ['-','-%']: subv,subh = 2,2
     428    if ope in ['-','-%','-_histo'] and len(namefiles) ==1 : subv,subh = 2,2
     429    elif ope in ['-','-%'] and len(namefiles)>1 : subv, subh = len(namefiles), 3
    380430 
    381431    #################################
    382432    ### Time loop for plotting device
    383     nplot = 1 ; error = False
     433    nplot = 1 ; error = False ; firstpass = True
    384434    if lstyle is not None: linecycler = cycle(lstyle)
    385435    else: linecycler = cycle(["-","--","-.",":"])
    … …  
    387437    while error is False:
    388438     
    389        print "********** NPLOT", nplot
     439       print "********** PLOT", nplot, " OF ",numplot
    390440       if nplot > numplot: break
    391441
    … …  
    393443       ## get all indexes to be taken into account for this subplot and then reduce field
    394444       ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
     445       if ope is not None:
     446           if fileref is not None:      index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(3*len(namefiles))  ## OK only 1 var,  see test in the beginning
     447           elif "var" in ope:           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
     448           elif "cat" in ope:           index_f = 0
     449       elif not firstpass:
     450          if len(namefiles) > 1 and len(var) == 1 and which == "unidim": pass
     451          else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
     452       else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
     453       time = all_time[index_f] ; vert = all_vert[index_f] ; lat = all_lat[index_f] ; lon = all_lon[index_f]
    395454       indexlon  = getsindex(sslon,(nplot-1)%nlon,lon)
    396455       indexlat  = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
    397        indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
    398        if ope is not None:
    399            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
    400            elif "var" in ope:           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
    401            elif "cat" in ope:           index_f = 0
    402        else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
    403        time = all_time[index_f]
     456       indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
     457       plotx=plot_x[index_f] ; ploty=plot_y[index_f]
    404458       if mrate is not None:                 indextime = None
    405459       else:                                 indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
    406460       ltst = None
    407        if typefile in ['meso'] and indextime is not None:  ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) ) 
     461       if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = localtime ( indextime, slon , nc)
    408462       print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
    409463       ##var = nc.variables["phisinit"][:,:]
    … …  
    419473       error = False
    420474       varname = all_varname[index_f]
     475       which=''
    421476       if varname:   ### what is shaded.
    422477           what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
    … …  
    432487           vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
    433488           if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind"
    434 
     489 
     490       if plotx is not None:
     491          plotx, error = reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
     492                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
     493          ploty, error = reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
     494                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
     495          which='xy'
    435496       #####################################################################
    436497       #if truc:
    … …  
    452513       ####################################################################
    453514       ### General plot settings
    454        changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) ## default for 1D plots: superimposed. to be reworked for better flexibility.
     515       changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) and (which != "xy") ## default for 1D plots: superimposed. to be reworked for better flexibility.
    455516       if changesubplot: subplot(subv,subh,nplot) #; subplots_adjust(wspace=0,hspace=0)
    456517       colorb = all_colorb[index_f]
    … …  
    474535                       if slat is not None: lonyeah = lon2d[0,:]
    475536                   what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
    476                          itime,what_I_plot, len(all_var[index_f].shape),vertmode)
     537                         itime,what_I_plot, len(all_var[index_f].shape),vertmode,redope)
    477538               ###
    478                if (fileref is not None) and (index_f == numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
     539               if (fileref is not None) and ((index_f+1)%3 == 0):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
    479540               else:                                                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    480541               #if (fileref is not None) and (index_f == numplot-1):    colorb = "RdBu_r"
    … …  
    493554                 if mrate is not None: iend=len(time)-1
    494555                 else:                 iend=0
    495                  imov = 0
    496                  if len(what_I_plot.shape) == 3:
    497                     if var2:               which = "contour" ## have to start with contours rather than shading
    498                     else:                  which = "regular"
     556                 imov = 0
     557                 if analysis in ['density','histo','fft','histodensity']: which="xy"
     558                 elif len(what_I_plot.shape) == 3:
     559                    if var2 and which == '':               which = "contour" ## have to start with contours rather than shading
     560                    elif which == '':                      which = "regular"
    499561                    if mrate is None:      errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
    500562                 elif len(what_I_plot.shape) == 2:
    501                     if var2:               which = "contour" ## have to start with contours rather than shading
    502                     else:                  which = "regular"
    503                     if mrate is not None:  which = "unidim"
    504                  elif len(what_I_plot.shape) == 1:
     563                    if var2 and which == '':               which = "contour" ## have to start with contours rather than shading
     564                    elif which == '':                      which = "regular"
     565                    if mrate is not None and which == '':  which = "unidim"
     566                 elif len(what_I_plot.shape) == 1 and which == '' :
    505567                    which = "unidim"
    506568                    if what_I_plot.shape[-1] == 1:      print "VALUE VALUE VALUE VALUE ::: ", what_I_plot[0] ; save = 'donothing'
     569##                 if which == "unidim" and len(namefiles) > 1: numplot = 1 # this case is similar to several vars from one file
    507570                 ##### IMOV LOOP #### IMOV LOOP
    508571                 while imov <= iend:
    … …  
    522585                    #if mrate is not None:     
    523586                    if mapmode == 1:
    524                             m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
    525                             x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
     587                        m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
     588                        x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
    526589                    if typefile in ['meso'] and mapmode == 1:   what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
    527590#                   if typefile in ['mesoideal']:    what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
    528591
    529592                    if which == "unidim":
    530                         if lbls is not None: lbl=lbls[nplot-1]
     593#                        if lbls is not None: lbl=lbls[nplot-1]
     594                        if lbls is not None: lbl=lbls[index_f]
    531595                        else:
    532596                           lbl = ""
    … …  
    549613                        if indextime is None and axtime is not None and xlab is None:    xlabel(axtime.upper()) ## define the right label
    550614                        if save == 'txt':  writeascii(np.transpose([x,np.transpose(what_I_plot)]),'profile'+str(nplot*1000+imov)+'.txt')
     615
     616                    elif which == "xy":
     617                        if lbls is not None: lbl=lbls[index_f]
     618                        else: lbl=None
     619                        if analysis is not None:
     620                           if analysis == 'histo':
     621                               if zelen == 1:
     622                                  mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.67, facecolor = 'green', label = lbls)
     623                                  mpl.pyplot.legend()
     624                                  mpl.pyplot.grid(True)
     625                                  mpl.pyplot.title(zetitle)
     626                               else:
     627                                  multiplot[index_f]=ploty.flatten()
     628                                  if index_f == zelen-1:
     629                                     if ope is not None: multiplot = getVar(multiplot,3*np.arange((index_f+1)/3)+2) ## we only compute histograms for the operation plots
     630                                     mpl.pyplot.hist(multiplot,bins=ndiv,normed=True, alpha=0.75, label = lbls)
     631                                     mpl.pyplot.legend()
     632                                     mpl.pyplot.grid(True)
     633                                     mpl.pyplot.title(zetitle)
     634                           elif analysis in ['density','histodensity']:
     635                                  if ope is not None and (index_f+1)%3 !=0: pass
     636                                  else:
     637                                     plotx = np.linspace(min(ploty.flatten()),max(ploty.flatten()),1000)
     638                                     density = gaussian_kde(ploty.flatten())
     639   #                                  density.covariance_factor = lambda : .25  # adjust the covariance factor to change the bandwidth if needed
     640   #                                  density._compute_covariance()
     641                                        # display the mean and variance of the kde:
     642                                     sample = density.resample(size=20000)
     643                                     n, (smin, smax), sm, sv, ss, sk = describe(sample[0])
     644                                     mpl.pyplot.plot(plotx,density(plotx), label = lbl+'\nmean: '+str(sm)[0:5]+'   std: '+str(np.sqrt(sv))[0:5]+'\nskewness: '+str(ss)[0:5]+'   kurtosis: '+str(sk)[0:5])
     645                                     if analysis == 'histodensity':  # plot both histo and density (to assess the rightness of the kernel density estimate for exemple) and display the estimated variance
     646                                        mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.30, label = lbl)
     647                                     if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
     648                           else:
     649                              plot(plotx,ploty,label = lbl)
     650                              if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
     651                              mpl.pyplot.grid(True)
     652                        else:
     653                           plot(plotx,ploty,label = lbl)
     654                           if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
     655                        if varname == 'hodograph':
     656                            a=0
     657                            for ii in np.arange(len(time)):
     658                               if a%6 == 0: mpl.pyplot.text(plotx[ii],ploty[ii],time[ii])
     659                               a=a+1
     660                            mpl.pyplot.grid(True)
    551661
    552662                    elif which == "regular":
    … …  
    590700
    591701                        if colorb not in ['nobar','onebar']:       
    592                             if (fileref is not None) and (index_f == numplot-1):   daformat = "%.3f"
     702                            if (fileref is not None) and ((index_f+1)%3 == 0):   daformat = "%.3f"
    593703                            elif mult != 1:                                        daformat = "%.1f"
    594704                            else:                                                  daformat = fmtvar(fvar.upper())
    … …  
    614724                                                              #200.         ## or csmooth=stride
    615725                        ### THIS IS A QUITE SPECIFIC PIECE (does not work for mesoscale files)
    616                         if ope == '-' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
     726                        if ope == '-_histo' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
    617727                            subplot(subv,subh,nplot+1)
    618728                            rcParams["legend.fontsize"] = 'xx-large'
    … …  
    632742                            plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
    633743                            legend(loc=0,frameon=False)
    634                             title("Histogram fig(1) "+ope+" fig(2)")
    635744                            subplot(subv,subh,nplot) # go back to last plot for title of contour difference
     745                        if ope is not None and "only" not in ope: title("fig(1) "+ope+" fig(2)")
     746                        elif ope is not None and "only" in ope: title("fig(1) "+ope[0]+" fig(2)")
    636747                           
    637748                    elif which == "contour":
    … …  
    647758                        if mapmode == 0:   
    648759                            what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
    649                                                               itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode )
     760                                                              itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode,redope)
    650761                            ## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape)
    651762                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
    … …  
    655766                            cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
    656767                            #clabel(cs, inline=0, fontsize = rcParams['font.size'], fmt="%.0f") #fmtvar(var2.upper()))
    657                     if which in ["regular","unidim"]:
    658 
    659                         if nplot > 1 and which == "unidim":
     768                    if which in ["regular","unidim","xy"]:
     769
     770                        if nplot > 1 and which in ["unidim","xy"]:
    660771                           pass  ## because we superimpose nplot instances
    661772                        else:
    … …  
    666777                           if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
    667778                           if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
    668                            if ylog:      mpl.pyplot.semilogy()
     779                           if ylog and not xlog:      mpl.pyplot.semilogy()
     780                           if xlog and not ylog:      mpl.pyplot.semilogx()
     781                           if xlog and ylog:
     782                                mpl.pyplot.xscale('log')
     783                                mpl.pyplot.yscale('log')
    669784                           if invert_y:  ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
    670785                           if xlab is not None: xlabel(xlab)
    … …  
    715830       else:
    716831            if fileref is not None:
    717                 if index_f is numplot-1:     plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
    718                 elif index_f is numplot-2:   plottitle = basename+' '+fileref
    719                 else:                        plottitle = basename+' '+all_namefile[index_f]
     832               if (index_f+1)%3 == 0:     plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
     833               elif (index_f+2)%3 == 0:   plottitle = basename+' '+fileref
     834               else:                        plottitle = basename+' '+all_namefile[index_f]
    720835            else:                            plottitle = basename+' '+all_namefile[index_f]
    721836       if mult != 1:                         plottitle = '{:.0e}'.format(mult) + "*" + plottitle
    … …  
    724839          else: plottitle = zetitle[0]
    725840          if titleref is "fill":             titleref=zetitle[index_f]
    726           if fileref is not None:
    727              if index_f is numplot-2:        plottitle = titleref
    728              if index_f is numplot-1:        plottitle = "fig(1) "+ope+" fig(2)"
     841          if fileref is not None and which != "xy":
     842             if (index_f+2)%3 == 0:        plottitle = titleref
     843             if (index_f+1)%3 == 0:        plottitle = "fig(1) "+ope+" fig(2)"
    729844#       if indexlon is not None:      plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
    730845#       if indexlat is not None:      plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
    … …  
    734849       if nplot >= numplot: error = True
    735850       nplot += 1
     851
     852       if len(namefiles) > 1 and len(var) == 1 and which == "unidim": index_f=index_f+1
     853       firstpass=False
    736854
    737855    if colorb == "onebar":

  • trunk/UTIL/PYTHON/pp.py

    r760 r763  
    5757          if opt.var is None:  opt.var = ["HGT"] #; opt.clb = "nobar" # otherwise breaks the automatic one-var file capability
    5858      elif typefile in ["geo"]:
    59           [lschar,zehour,zehourin] = ["",0,0]
     59          lschar=""
    6060          if opt.var is None:  opt.var = ["HGT_M"] ; opt.clb = "nobar"
    61       else:                                       
    62           [lschar,zehour,zehourin] = ["",0,0]
     61      else:                                     
     62          lschar=""
    6363          if opt.var is None: 
    6464             opt.var = ["phisinit"] ; opt.clb = "nobar"
    … …  
    106106            if opt.winds            :  zefields = 'uvmet'
    107107            elif opt.var[j] in ['deltat','DELTAT'] : zefields = 'tk,TSURF'
     108            elif opt.var[j] in ['uv','UV','hodograph','hodograph_2'] : zefields = 'U,V'
    108109            else                    :  zefields = ''
    109110            ### var or no var
    … …  
    163164                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevars,\
    164165                clb=separatenames(opt.clb),winds=opt.winds,\
    165                 addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
     166                addchar=lschar,vmin=zevmin,vmax=zevmax,\
    166167                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
    167168                hole=opt.hole,save=opt.save,\
    … …  
    171172                outputname=opt.out,resolution=opt.res,\
    172173                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
    173                 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
     174                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,xlog=opt.logx,yintegral=opt.column,\
    174175                blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\
    175176                mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\
    176177                redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels),\
    177178                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),facwind=opt.facwind,\
    178                 trycol=opt.trycol,streamflag=opt.stream)
     179                trycol=opt.trycol,streamflag=opt.stream,analysis=opt.analysis)
    179180        print 'DONE: '+name
    180181        system("rm -f to_be_erased")
Note: See TracChangeset for help on using the changeset viewer.