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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.