FLEXPART Testing Environment CTBTO WO8
 All Classes Namespaces Files Functions Variables Pages
Classes | Functions | Variables
flextest.flexread.pflexible Namespace Reference

Classes

class  BinaryFile
 HELPER FUNCTIONS ##########. More...
 
class  Structure
 
class  Header
 

Functions

def read_command
 Functions for reading FLEXPART output #####. More...
 
def read_releases
 
def read_trajectories
 
def curtain_for_flightrack
 
def curtain_for_line
 
def groundlevel_for_line
 
def curtain_agltoasl
 
def read_agespectrum
 
def save_spectrum
 
def gridarea
 
def _get_header_version
 
def read_header
 
def _readV6
 
def _readgrid_noFF
 Grid Reading Routines ###############. More...
 
def _readgridBF
 
def _read_headerFF
 
def read_grid
 
def readgridV8
 
def readgridV6
 
def monthly_footprints
 
def fill_grids
 
def read_emissions
 
def get_slabs
 
def _plot_dropm
 FLEXPART Plotting Functions ########. More...
 
def plot_trajchar
 
def plot_releases
 
def plot_spectra
 
def plot_agespectra
 
def plot_clusters
 
def plot_trajectory_ellipses
 
def plot_markers
 
def plot_trajectory
 
def plot_at_level
 
def plot_totalcolumn
 
def plot_sourcecontribution
 
def plot_sensitivity
 
def plot_curtain
 
def plot_METDATA
 
def closest
 
def data_range
 
def _normdatetime
 
def _log_clevs
 
def add_nilu_logo
 
def _cum_spec
 
def _getfile_lines
 
def _gen_MapPar_fromHeader
 
def _gen_altitude_color
 
def _gen_daylabels
 
def _datarange
 
def _genEllipse
 
def _shout
 
def _gen_flexpart_colormap
 
def main
 Below for use from command line ##########################################. More...
 

Variables

string __version__ = '0.9.5'
 
tuple __path__ = os.path.abspath(os.curdir)
 
 readheader = read_header
 
 readheaderV8 = read_header
 
 readheaderV6 = read_header
 
 readgridBF = _readgrid_noFF
 
 plot_footprint = plot_at_level
 
tuple start_time = time.time()
 
tuple parser
 
string default = 'verbose output'
 
tuple exit_code = main()
 LateX #####. More...
 

Function Documentation

def flextest.flexread.pflexible._cum_spec (   inspectra,
  cum = True 
)
private 
helper function for the plot_spectra routines

Definition at line 4922 of file pflexible.py.

4922 
4923 def _cum_spec(inspectra, cum=True):
4924  """ helper function for the plot_spectra routines """
4925 
4926  if cum == 'norm':
4927  # # Normalize the data so it fills to 100%
4928  # spectra = np.zeros(inspectra.shape)
4929  spectra = (inspectra.transpose() / np.sum(inspectra, axis=1)).transpose()
4930  spectra = np.cumsum(spectra[:, :], axis=1)
4931  # sums = np.sum(inspectra,axis=1)
4932  # for i,elem in enumerate(inspectra):
4933  # spectra[i,:] = elem/sums[i]
4934  elif cum is True and cum != 'norm':
4935  spectra = np.cumsum(inspectra[:, :], axis=1)
4936  else:
4937  spectra = inspectra
4938 
4939  return spectra
def flextest.flexread.pflexible._datarange (   H,
  G,
  index = None 
)
private 
returns max and min for footprint and total column
    for a given range of releases [default is all]

Definition at line 4992 of file pflexible.py.

4992 
4993 def _datarange(H, G, index=None):
4994  """ returns max and min for footprint and total column
4995  for a given range of releases [default is all] """
4996  Heightnn = H['Heightnn']
4997  fpmax = -999.
4998  if index == None:
4999  seek = range(G.shape[-1])
5000  else:
5001  seek = [index]
5002 
5003  for i in seek:
5004  zpmax = np.max(G[:, :, 0, i] / Heightnn[:, :, 0])
5005  if zpmax > fpmax:
5006  fpmax = zpmax
5007  # print fpmax
5008  tcmax = np.max(G)
5009  return ((0, fpmax), (0, tfcmax))
5010 
5011 

Here is the call graph for this function:

def flextest.flexread.pflexible._gen_altitude_color (   p,
  altlim = (0, 5000,
  numcon = 1,
  cmap_id = 'jet' 
)
private 
Generates color based on altlim = (min,max)

Definition at line 4966 of file pflexible.py.

4966 
4967 def _gen_altitude_color(p, altlim=(0, 5000), numcon=1, cmap_id='jet'):
4968  """
4969  Generates color based on altlim = (min,max)
4970  """
4971  norm = mpl.colors.Normalize(vmin=altlim[0], vmax=altlim[1])
4972  cmap = plt.cm.get_cmap(cmap_id)
4973  # p = p/altmax
4974  cm = cmap(norm(p))
4975  # print p, cm
4976  return cm
4977 

Here is the caller graph for this function:

def flextest.flexread.pflexible._gen_daylabels (   P,
  H = None,
  dt = 86400 
)
private 
Uses H.loutstep to calculate 'days back' for plotting on clusters and trajectories.

Definition at line 4978 of file pflexible.py.

4978 
4979 def _gen_daylabels(P, H=None, dt=86400):
4980  """
4981  Uses H.loutstep to calculate 'days back' for plotting on clusters and trajectories.
4982  """
4983  if isinstance(P, int):
4984 
4985  if H:
4986  dt = abs(H.loutstep)
4987  return str(1 + int(abs(P)) / dt)
4988  else:
4989  return [str(1 + int(abs(p)) / dt) for p in P]
4990 
4991 

Here is the caller graph for this function:

def flextest.flexread.pflexible._gen_flexpart_colormap (   ctbfile = None,
  colors = None 
)
private 
Generate the ast colormap for FLEXPART

Definition at line 5068 of file pflexible.py.

5068 
5069 def _gen_flexpart_colormap(ctbfile=None, colors=None):
5070  """
5071  Generate the ast colormap for FLEXPART
5072  """
5073  from matplotlib.colors import ListedColormap
5074  if ctbfile:
5075  try:
5076  colors = np.loadtxt(ctbfile)
5077  except:
5078  print "WARNING: cannot load ctbfile. using colors"
5079  if colors:
5080  name = 'user_colormap'
5081  if not colors:
5082  # # AST Colorset for FLEXPART
5083  colors = [1.0000000e+00, 1.0000000e+00, 1.0000000e+00
5084  , 9.9607843e-01, 9.1372549e-01, 1.0000000e+00
5085  , 9.8431373e-01, 8.2352941e-01, 1.0000000e+00
5086  , 9.6470588e-01, 7.1764706e-01, 1.0000000e+00
5087  , 9.3333333e-01, 6.0000000e-01, 1.0000000e+00
5088  , 8.9019608e-01, 4.4705882e-01, 1.0000000e+00
5089  , 8.3137255e-01, 2.0000000e-01, 1.0000000e+00
5090  , 7.5686275e-01, 0.0000000e+00, 1.0000000e+00
5091  , 6.6274510e-01, 0.0000000e+00, 1.0000000e+00
5092  , 5.4901961e-01, 0.0000000e+00, 1.0000000e+00
5093  , 4.0784314e-01, 0.0000000e+00, 1.0000000e+00
5094  , 2.4705882e-01, 0.0000000e+00, 1.0000000e+00
5095  , 7.4509804e-02, 0.0000000e+00, 1.0000000e+00
5096  , 0.0000000e+00, 2.8235294e-01, 1.0000000e+00
5097  , 0.0000000e+00, 4.8627451e-01, 1.0000000e+00
5098  , 0.0000000e+00, 6.3137255e-01, 1.0000000e+00
5099  , 0.0000000e+00, 7.4509804e-01, 1.0000000e+00
5100  , 0.0000000e+00, 8.4705882e-01, 1.0000000e+00
5101  , 0.0000000e+00, 9.3725490e-01, 1.0000000e+00
5102  , 0.0000000e+00, 1.0000000e+00, 9.7647059e-01
5103  , 0.0000000e+00, 1.0000000e+00, 8.9411765e-01
5104  , 0.0000000e+00, 1.0000000e+00, 8.0000000e-01
5105  , 0.0000000e+00, 1.0000000e+00, 6.9019608e-01
5106  , 0.0000000e+00, 1.0000000e+00, 5.6470588e-01
5107  , 0.0000000e+00, 1.0000000e+00, 4.0000000e-01
5108  , 0.0000000e+00, 1.0000000e+00, 0.0000000e+00
5109  , 3.9607843e-01, 1.0000000e+00, 0.0000000e+00
5110  , 5.6470588e-01, 1.0000000e+00, 0.0000000e+00
5111  , 6.9019608e-01, 1.0000000e+00, 0.0000000e+00
5112  , 7.9607843e-01, 1.0000000e+00, 0.0000000e+00
5113  , 8.9411765e-01, 1.0000000e+00, 0.0000000e+00
5114  , 9.7647059e-01, 1.0000000e+00, 0.0000000e+00
5115  , 1.0000000e+00, 9.4509804e-01, 0.0000000e+00
5116  , 1.0000000e+00, 8.7450980e-01, 0.0000000e+00
5117  , 1.0000000e+00, 7.9215686e-01, 0.0000000e+00
5118  , 1.0000000e+00, 7.0588235e-01, 0.0000000e+00
5119  , 1.0000000e+00, 6.0392157e-01, 0.0000000e+00
5120  , 1.0000000e+00, 4.8235294e-01, 0.0000000e+00
5121  , 1.0000000e+00, 3.1372549e-01, 0.0000000e+00
5122  , 1.0000000e+00, 0.0000000e+00, 1.4901961e-01
5123  , 1.0000000e+00, 0.0000000e+00, 3.3333333e-01
5124  , 1.0000000e+00, 0.0000000e+00, 4.4705882e-01
5125  , 1.0000000e+00, 0.0000000e+00, 5.3725490e-01
5126  , 1.0000000e+00, 0.0000000e+00, 6.1176471e-01
5127  , 9.7647059e-01, 0.0000000e+00, 6.6666667e-01
5128  , 8.9411765e-01, 0.0000000e+00, 6.6666667e-01
5129  , 7.9607843e-01, 0.0000000e+00, 6.3921569e-01
5130  , 6.9019608e-01, 0.0000000e+00, 5.9215686e-01
5131  , 5.6470588e-01, 0.0000000e+00, 5.0980392e-01
5132  , 3.9607843e-01, 0.0000000e+00, 3.8039216e-01]
5133  colors = np.reshape(colors, (-1, 3))
5134  name = 'flexpart_cmap'
5135  cmap = ListedColormap(colors, name)
5136  return cmap
5137 
5138 
def flextest.flexread.pflexible._gen_MapPar_fromHeader (   H)
private 
Define some default map parameters from the Header File.

Definition at line 4950 of file pflexible.py.

4950 
4951 def _gen_MapPar_fromHeader(H):
4952  """
4953  Define some default map parameters from the Header File.
4954  """
4955 
4956  MapPar = Structure()
4957  MapPar.llcrnrlat = H.outlat0
4958  MapPar.llcrnrlon = H.outlon0
4959  MapPar.urcrnrlat = H.outlat0 + (H.dyout * H.numygrid)
4960  MapPar.urcrnrlon = H.outlon0 + (H.dxout * H.numxgrid - 1)
4961  for k in MapPar.keys():
4962  print k, MapPar[k]
4963 
4964  return MapPar
4965 
def flextest.flexread.pflexible._genEllipse (   data,
  m,
  sizescale = 20000,
  altlim = None 
)
private 
Generates ellipses based on input array 'data'. Requires basemap instance, 'm'. NOTE::

    data = [itime,x,y,z,[size]]

    r,c = data.shape

    if c == 5:
        size function of data[:,4]
    else:
        size = 1*sizescale

uses functions:

    * _gen_altitude_color
    * _gen_daylabels

for labeling/color of ellipses

Definition at line 5013 of file pflexible.py.

5013 
5014  altlim=None):
5015  """
5016 
5017  Generates ellipses based on input array 'data'. Requires basemap instance, 'm'. NOTE::
5018 
5019  data = [itime,x,y,z,[size]]
5020 
5021  r,c = data.shape
5022 
5023  if c == 5:
5024  size function of data[:,4]
5025  else:
5026  size = 1*sizescale
5027 
5028  uses functions:
5029 
5030  * _gen_altitude_color
5031  * _gen_daylabels
5032 
5033  for labeling/color of ellipses
5034 
5035  """
5036  r, c = data.shape
5037  if altlim == None:
5038  altlim = (np.min(data[:, 3]), np.max(data[:, 3]))
5039 
5040  if c == 5:
5041 
5042  ell = [(Ellipse(xy=np.array(m(p[1], p[2])),
5043  width=p[4] * sizescale, height=p[4] * sizescale,
5044  angle=360,
5045  facecolor=_gen_altitude_color(p[3], altlim=altlim, cmap_id='gray'), \
5046  label=_gen_daylabels(p[0])), \
5047  np.array(m(p[1], p[2]))) for p in data]
5048  # np.array( m(p[1],p[2])) ) for p in data]
5049  else:
5050  ell = [(Ellipse(xy=np.array(m(p[1], p[2])),
5051  width=1e4 * sizescale, height=1e4 * sizescale,
5052  angle=360, \
5053  facecolor=_gen_altitude_color(p[3], altlim=altlim, cmap_id='gray'), \
5054  label=_gen_daylabels(p[0])), \
5055  np.array(m(p[1], p[2]))) for p in data]
5056 
5057  return ell
5058 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible._get_header_version (   bf)
private 
Open and read the binary file (bf) header only to the point of 
the Flexpart version string, return the string, seek to the
start of the file.

Definition at line 794 of file pflexible.py.

794 
795 def _get_header_version(bf):
796  """
797  Open and read the binary file (bf) header only to the point of
798  the Flexpart version string, return the string, seek to the
799  start of the file.
800 
801  """
802  try:
803  bf = BinaryFile(bf)
804  except:
805  bf = bf
806 
807  ret = bf.tell()
808  bf.seek(12) # start of version string
809  version = bf.read('13S')
810  bf.seek(ret)
811 
812  return version
flextest.flexread.pflexible.BinaryFile
HELPER FUNCTIONS ##########.
Definition: pflexible.py:4484

Here is the caller graph for this function:

def flextest.flexread.pflexible._getfile_lines (   infile)
private 
returns all lines from a file or file string
reverts to beginning of file.

Definition at line 4940 of file pflexible.py.

4940 
4941 def _getfile_lines(infile):
4942  """ returns all lines from a file or file string
4943  reverts to beginning of file."""
4944 
4945  if isinstance(infile, str):
4946  return file(infile, 'r').readlines()
4947  else:
4948  infile.seek(0)
4949  return infile.readlines()

Here is the caller graph for this function:

def flextest.flexread.pflexible._log_clevs (   dat_min,
  dat_max 
)
private 
## create logorithmic color scale
## method uses strings to get the magnitude of variable
#ss = "%1.2e" % (dat_max)
#dmx = int('%s%s' % (ss[-3],int(ss[-2:])))
#dmx=float(len(str(int(np.ceil(dat_max)))))
#if data_range is None:
#     dmn = int(np.round(np.log((1e-10*(dat_max-dat_min)))))
#    ss = "%1.2e" % (1e-10*(dat_max-dat_min))
#else:
#    ss = "%1.2e" % (dat_min)
#dmn = int('%s%s' % (ss[-3],int(ss[-2:])))

Definition at line 4875 of file pflexible.py.

4875 
4876 def _log_clevs(dat_min, dat_max):
4877  """
4878  ## create logorithmic color scale
4879  ## method uses strings to get the magnitude of variable
4880  #ss = "%1.2e" % (dat_max)
4881  #dmx = int('%s%s' % (ss[-3],int(ss[-2:])))
4882  #dmx=float(len(str(int(np.ceil(dat_max)))))
4883  #if data_range is None:
4884  # dmn = int(np.round(np.log((1e-10*(dat_max-dat_min)))))
4885  # ss = "%1.2e" % (1e-10*(dat_max-dat_min))
4886  #else:
4887  # ss = "%1.2e" % (dat_min)
4888  #dmn = int('%s%s' % (ss[-3],int(ss[-2:])))
4889  """
4890 
4891  if dat_max > 0:
4892  dmx = int(np.round(np.log10(dat_max))) + 1
4893  else:
4894  print 'dat_max not positive'
4895  print dat_max
4896  dmx = 1
4897 
4898  if dat_min > 0:
4899  dmn = int(np.round(np.log10(dat_min)))
4900  elif dat_min == 0. or np.isnan(dat_min):
4901  print 'dat_min not positive'
4902  print dat_min
4903  dmn = dmx - 3
4904 
4905 
4906  # # create equally spaced range
4907  if dmx == dmn:
4908  dmx = dmn + 1
4909  clevs = np.logspace(dmn, dmx, 100)
4910 
4911  return clevs
4912 
4913 
4914 
def flextest.flexread.pflexible._normdatetime (   Y,
  M,
  D,
  h,
  m,
  s 
)
private

Definition at line 4872 of file pflexible.py.

4872 
4873 def _normdatetime(Y, M, D, h, m, s):
4874  return datetime.datetime(Y, M, D) + datetime.timedelta(hours=h, minutes=m, seconds=s)
def flextest.flexread.pflexible._plot_dropm ( )
private

FLEXPART Plotting Functions ########.

Definition at line 2622 of file pflexible.py.

2622 
2623 def _plot_dropm():
2624  try:
2625  del FIGURE.m
2626  plt.close('all')
2627  del m
2628  "m dropped"
2629  return 1
2630  except:
2631  print 'Cannot delete m, does not exist.'
2632  return 0
2633 
2634 
2635 
flextest.flexread.pflexible._plot_dropm
def _plot_dropm
FLEXPART Plotting Functions ########.
Definition: pflexible.py:2622

Here is the call graph for this function:

def flextest.flexread.pflexible._read_headerFF (   pathname,
  h = None,
  maxpoint = 800000,
  maxspec = 4,
  maxageclass = 20,
  nxmax = 722,
  nymax = 362,
  nzmax = 40,
  verbose = True 
)
private 
DEPRECATED

Called from read_header if readp_ff is True, uses FortFlex.readheader 
to read all releasepoints.

This function is dependant on the FortFlex.so module
see FortFlex.f and the f2py directory

Definition at line 1758 of file pflexible.py.

1758 
1759  nxmax=722, nymax=362, nzmax=40, verbose=True):
1760  """
1761  DEPRECATED
1762 
1763  Called from read_header if readp_ff is True, uses FortFlex.readheader
1764  to read all releasepoints.
1765 
1766  This function is dependant on the FortFlex.so module
1767  see FortFlex.f and the f2py directory
1768  """
1769  try:
1770  from FortFlex import readheader
1771  except:
1772  print "Error with FortFlex.readheader, use read_header"
1773 
1774  headervars = ['numxgrid', 'numygrid', 'numzgrid', 'outlon0', 'outlat0', 'compoint', \
1775  'dxout', 'dyout', 'outheight', 'ibdate', 'ibtime', 'loutstep', \
1776  'nspec', 'nageclass', 'lage', 'ireleasestart', 'ireleaseend', \
1777  'numpoint', 'xpoint', 'ypoint', 'zpoint1', 'zpoint2', 'heightnn', 'area', \
1778  'maxpoint', 'maxspec', 'maxageclass', 'nxmax', 'nymax', 'nzmax', \
1779  'npart', 'kind', 'lage', 'loutaver', 'loutsample', 'yyyymmdd', 'hhmmss', 'method']
1780  if h == None:
1781  h = Structure()
1782 
1783  if verbose:
1784  print """Reading Header with:
1785  maxpoint : %s
1786  maxspec : %s
1787  maxageclass : %s
1788  nxmax : %s
1789  nymax : %s
1790  nzmax : %s
1791  """ % (maxpoint, maxspec, maxageclass, nxmax, nymax, nzmax)
1792 
1793 
1794  numxgrid, numygrid, numzgrid, outlon0, outlat0, dxout, dyout, outheight, \
1795  ibdate, ibtime, loutstep, nspec, nageclass, lage, ireleasestart, ireleaseend, \
1796  numpoint, xpoint, ypoint, zpoint1, zpoint2, heightnn, area, compoint, species_f, \
1797  npart, kind, loutaver, loutsample, yyyymmdd, hhmmss, method = \
1798  readheader(pathname, maxpoint, maxspec, maxageclass, nxmax, nymax, nzmax)
1799 
1800  for v in headervars:
1801  if v not in h.keys():
1802  exec("h.%s = %s" % (v, v))
1803  return h

Here is the caller graph for this function:

def flextest.flexread.pflexible._readgrid_noFF (   H,
  kwargs 
)
private

Grid Reading Routines ###############. 

accepts a header dictionary as input, returns dictionary of Grid values
%===========================================
%
%-------------------------------------------
% input
%   - H: required header dictionary
%
% optional (most is obtained from header dict)
%   - date: which yyyymmddhhmmss from dates
%   - unit: 'conc', 'pptv', ['time'], 'footprint'
%   - nspec_ret: numspecies 
%   - pspec_ret: 
%   - age_ret:
%   - time_ret: 
%   - nested: nested = True
%
% 
% output
%   - grid dictionary object
%   - fail indicator (-1=fail, 0=success)
%-------------------------------------------
% FLEXPART python import routines
%-------------------------------------------
% last changes: JFB, 10.10.2008
%===========================================

Definition at line 1336 of file pflexible.py.

1336 
1337 def _readgrid_noFF(H, **kwargs):
1338  """ accepts a header dictionary as input, returns dictionary of Grid values
1339  %===========================================
1340  %
1341  %-------------------------------------------
1342  % input
1343  % - H: required header dictionary
1344  %
1345  % optional (most is obtained from header dict)
1346  % - date: which yyyymmddhhmmss from dates
1347  % - unit: 'conc', 'pptv', ['time'], 'footprint'
1348  % - nspec_ret: numspecies
1349  % - pspec_ret:
1350  % - age_ret:
1351  % - time_ret:
1352  % - nested: nested = True
1353  %
1354  %
1355  % output
1356  % - grid dictionary object
1357  % - fail indicator (-1=fail, 0=success)
1358  %-------------------------------------------
1359  % FLEXPART python import routines
1360  %-------------------------------------------
1361  % last changes: JFB, 10.10.2008
1362  %===========================================
1363  """
1364  # if os.sys.platform != 'win32':
1365  # try:
1366  # from FortFlex import readgrid, sumgrid
1367  # useFortFlex = 1
1368  # print 'using FortFlex'
1369  # except:
1370  # useFortFlex = 0
1371  # print 'Cannot find FortFlex.so'
1372  useFortFlex = 0
1373 
1374  if 'date' in kwargs.keys():
1375  date = kwargs['date']
1376  # else: date = H['ibtime']
1377  else: date = None
1378 
1379  if 'unit' in kwargs.keys():
1380  unitname = kwargs['unit']
1381  else: unitname = 'time'
1382 
1383  units = ['conc', 'pptv', 'time', 'footprint']
1384  unit = units.index(unitname)
1385 
1386  if 'nspec_ret' in kwargs.keys():
1387  nspec_ret = kwargs['nspec_ret']
1388  else: nspec_ret = 1
1389 
1390  if 'pspec_ret' in kwargs.keys():
1391  pspec_ret = kwargs['ppsec_ret']
1392  else: pspec_ret = 1
1393 
1394  if 'age_ret' in kwargs.keys():
1395  age_ret = kwargs['age_ret']
1396  else: age_ret = 1
1397 
1398  if 'nested' in kwargs.keys():
1399  nested = kwargs['nested']
1400  else: nested = False
1401 
1402  if 'time_ret' in kwargs.keys():
1403  time_ret = kwargs['time_ret']
1404  else:
1405  time_ret = range(len(H['available_dates']))
1406 
1407  if 'scaledepo' in kwargs.keys():
1408  scaledepo = kwargs['scaledepo']
1409  else:
1410  scaledepo = 1.0
1411 
1412  if 'scaleconc' in kwargs.keys():
1413  scaleconc = kwargs['scaleconc']
1414  else:
1415  scaleconc = 1.0
1416 
1417  if 'decaycons' in kwargs.keys():
1418  decaycons = kwargs['decaycons']
1419  else:
1420  decaycons = 99999999990
1421  fail = -1
1422 
1423  # set filenames
1424  prefix = ['grid_conc_', 'grid_pptv_', \
1425  'grid_time_', 'footprint_total', \
1426  'grid_conc_nest_', 'grid_pptv_nest_', \
1427  'grid_time_nest_', 'footprint_total_nest']
1428 
1429  # local functions for reading binary data and creating grid dump
1430 # skip = lambda n=8 : f.read(n)
1431 # getbin = lambda fmt,n=1 : struct.unpack(fmt*n,f.read(struct.calcsize(fmt)))[0]
1432 
1433  # Utility functions
1434  skip = lambda n = 8 : f2.seek(n, 1)
1435  getbin = lambda dtype, n = 1 : f2.read(dtype, (n,))
1436 
1437 
1438  def getdump(n, fmt='f'):
1439  """ function to get the dump values for the sparse format """
1440  # print n
1441  skip()
1442  # Dfmt=[fmt]*n
1443 # a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt]
1444  a = f2.read(fmt, n)
1445 # dumplist=[a[j][0] for j in range(len(a))]
1446  # dumplist=[a[j] for j in range(len(a))]
1447  return a # dumplist
1448 
1449  def key2var(D, key):
1450  cmd = "global %s; %s = D['%s'];" % (key, key, key)
1451  exec(cmd)
1452 
1453  def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage):
1454  """ function to dump sparse elements into grid """
1455  ii = 0
1456  fact = 1
1457  for ir in range(cnt_r):
1458  if dmp_r[ir] * fact > 0:
1459  n = dmp_i[ii]; ii = ii + 1; fact = fact * -1.
1460  else:
1461  n = n + 1
1462 
1463  kz = n / (numxgrid * numygrid)
1464  jy = (n - kz * numxgrid * numygrid) / numxgrid
1465  ix = n - numxgrid * numygrid * kz - numxgrid * jy
1466  # print "n ==> ix,jy,kz,k,nage"
1467  # print "%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage)
1468  # print grd.shape
1469  # print grd[0,0,0,0,0]
1470  grd[ix, jy, kz - 1, k, nage] = abs(dmp_r[ir])
1471  return grd # flipud(grd.transpose())
1472 
1473  G = {}
1474  # get values from header file
1475  H['nageclass'] = H['numageclasses']
1476  headervars = ['nested', 'nspec', 'numxgrid', 'numygrid', 'numzgrid', 'nageclass', \
1477  'available_dates', 'pathname', 'decayconstant', 'numpoint', 'Area', 'Heightnn']
1478 
1479  for k in headervars:
1480  key2var(H, k)
1481 
1482  # create zero arrays for datagrid and zplot
1483  datagrid = np.zeros((numxgrid, numygrid, numzgrid, numpoint), dtype=np.float32)
1484  zplot = np.empty((numxgrid, numygrid, numzgrid, numpoint))
1485 
1486  #--------------------------------------------------
1487  # Loop over all times, given in field H['dates']
1488  #--------------------------------------------------
1489  for ks in range(nspec_ret, nspec_ret + 1):
1490  # print 'processing: ' + str(H['dates'][date_i]) + ' ' + str(ks)
1491  # print 'processing: ' + str(date) + ' ' + str(ks)
1492  specs = '_' + str(ks).zfill(3)
1493  fpmax = -999.
1494  if date == None and time_ret == None:
1495  get_dates = [available_dates[0]]
1496  elif time_ret == None:
1497  get_dates = []
1498  date = date.strip().split(',')
1499  for d in date:
1500  try:
1501  get_dates.append(available_dates.index(d))
1502  except:
1503  shout("Cannot find date: %s in H['dates']\n" % d)
1504  else:
1505  get_dates = available_dates
1506 
1507  #print 'getting grid for: ', get_dates
1508 
1509  for date_i in range(len(get_dates)):
1510 
1511  if unit != 4:
1512  datestring = H['available_dates'][date_i]
1513  # datestring = date
1514  filename = os.path.join(H['pathname'], \
1515  prefix[(unit) + (nested * 4)] + datestring + specs)
1516  else:
1517  filename = os.path.join(H['pathname'], \
1518  prefix[(unit) + (nested * 4)])
1519 
1520  # print 'reading: ' + filename
1521 
1522  if os.path.exists(filename):
1523  # print nspec,numxgrid,numygrid,numzgrid,nageclass,scaledepo,scaleconc,decaycons
1524  # print ks
1525  # print date
1526 
1527  if useFortFlex == 1:
1528  # print 'Using FortFLEX'
1529  #####################################################################################
1530  # # FORTRAN WRAPPER CODE - only works on linux
1531  # #
1532  # # USAGE: grid = readgrid(filegrid,numxgrid,numygrid,numzgrid,\
1533  # # nspec,nageclass,scaleconc,decayconstant)
1534  # #
1535  # # zplot = sumgrid(zplot,grid,\
1536  # # area,heightnn,\
1537  # # [numxgrid,numygrid,numzgrid,numpoint,nageclass])
1538  # #
1539  # #
1540  # # NOTE: numpoint = number of releases, ageclass always 1 for backward
1541  # #
1542  # # RETURN: grid(numxgrid,numygrid,numzgrid,numpoint,nageclass)
1543  #####################################################################################
1544 
1545  concgrid = readgrid(filename, numxgrid, numygrid, numzgrid, \
1546  numpoint, nageclass, \
1547  scaleconc, decayconstant)
1548 
1549  # contribution[:,:,:] = concgrid[:,:,:,:,0]
1550  print np.min(concgrid)
1551  print np.max(concgrid)
1552 
1553  # altitude = 50000
1554  zplot = sumgrid(zplot, concgrid, \
1555  H.area, H.Heightnn)
1556 
1557 
1558  else:
1559  dat_cnt = 0
1560  nage = 0
1561  # read data:
1562  # #datagrid=np.zeros((numxgrid,numygrid,numzgrid[nspec-1],nspec,nageclass),np.float)
1563  datagrid = np.zeros((numxgrid, numygrid, numzgrid, 1, 1), np.float)
1564  # f = file(filename, 'rb')
1565  # print filename
1566  f2 = BinaryFile(filename, order='fortran')
1567  skip(4)
1568  G['itime'] = getbin('i');
1569  print H['available_dates'][date_i]
1570 
1571  # Read Wet Depostion
1572  skip()
1573  cnt_i = getbin('i')
1574  dmp_i = getdump(cnt_i, 'i')
1575  skip()
1576  cnt_r = getbin('i')
1577  dmp_r = getdump(cnt_r)
1578  # wet=dumpgrid(dmp_i, cnt_r, dmp_r, datagrid, ks-1, nage)
1579  # Read Dry Deposition
1580  skip()
1581  cnt_i = getbin('i')
1582  dmp_i = getdump(cnt_i, 'i')
1583  skip()
1584  cnt_r = getbin('i')
1585  dmp_r = getdump(cnt_r)
1586  # dry=dumpgrid(dmp_i, cnt_r, dmp_r, datagrid, ks-1, nage)
1587 
1588  # Read Concentrations
1589  skip()
1590  cnt_i = getbin('i')
1591  dmp_i = getdump(cnt_i, 'i')
1592  skip()
1593  cnt_r = getbin('i')
1594  dmp_r = getdump(cnt_r)
1595  # print dmp_i, cnt_r, dmp_r, datagrid, ks-1, nage
1596  concgrid = dumpgrid(dmp_i, cnt_r, dmp_r, datagrid, ks - 1, nage)
1597 
1598  # G[H['ibtime']].append(concgrid)
1599  G[H['ibtime']] = concgrid
1600  f2.close()
1601  fail = 0
1602  else:
1603  print "\n\n INPUT ERROR: Could not find file: %s" % filename
1604  raise IOError('No file: %s' % filename)
1605 
1606  if useFortFlex == 1: G = zplot
1607  return G
flextest.flexread.pflexible._readgrid_noFF
def _readgrid_noFF
Grid Reading Routines ###############.
Definition: pflexible.py:1336
flextest.flexread.pflexible.BinaryFile
HELPER FUNCTIONS ##########.
Definition: pflexible.py:4484
def flextest.flexread.pflexible._readgridBF (   H,
  filename 
)
private 
Read grid using BinaryFile class

Definition at line 1610 of file pflexible.py.

1610 
1611 def _readgridBF(H, filename):
1612  """ Read grid using BinaryFile class"""
1613  # Utility functions
1614  skip = lambda n = 8 : f2.seek(n, 1)
1615  getbin = lambda dtype, n = 1 : f2.read(dtype, (n,))
1616 
1617 
1618  def getdump(n, fmt='f'):
1619  """ function to get the dump values for the sparse format """
1620  skip()
1621  # Dfmt=[fmt]*n
1622 # a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt]
1623  a = f2.read(fmt, n)
1624  ### DJM HACK - 2013-09-02 - without the following, 'a' ends up being
1625  ### a float rather than a numpy array
1626  if n==1:
1627  a = np.array([a])
1628 # dumplist=[a[j][0] for j in range(len(a))]
1629  # dumplist=[a[j] for j in range(len(a))]
1630  return a # dumplist
1631 
1632  def key2var(D, key):
1633  cmd = "global %s; %s = D['%s'];" % (key, key, key)
1634  exec(cmd)
1635 
1636  def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny):
1637  """ function to dump sparse elements into grid, fall back method if
1638  pflexcy.so module (cython) not available -- it is SLOW """
1639  conc = False
1640  if len(grd.shape) == 5:
1641  conc = True
1642  ii = 0
1643  fact = 1
1644  pos = 0
1645  '''
1646  print 'CNT_R: ' + str(cnt_r)
1647  print 'DMP_R: ' + str(dmp_r)
1648  print 'TYPE DMP_R: ' + str(type(dmp_r))
1649  '''
1650  for ir in range(cnt_r):
1651 
1652  if conc:
1653  if dmp_r[ir] * fact > 0:
1654  n = dmp_i[ii]
1655  ii = ii + 1
1656  fact = fact * -1.
1657  else:
1658  n = n + 1
1659 
1660  kz = n / (H.numxgrid * H.numygrid)
1661  jy = (n - kz * H.numxgrid * H.numygrid) / H.numxgrid
1662  ix = n - H.numxgrid * H.numygrid * kz - H.numxgrid * jy
1663  grd[ix, jy, kz - 1, k, nage] = abs(dmp_r[ir])
1664 
1665 #
1666 # print "n ==> ix,jy,kz,k,nage"
1667 # print "%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage)
1668 # print grd.shape
1669 # print grd[0,0,0,0,0]
1670 
1671  else:
1672  if dmp_r[ir] * fact > 0:
1673  n = dmp_i[ii]
1674  ii = ii + 1
1675  fact = fact * -1.
1676  else:
1677  n = n + 1
1678  jy = n / H.numxgrid
1679  ix = n - H.numxgrid * jy
1680  grd[ix, jy, k, nage] = abs(dmp_r[ir])
1681 
1682  return grd # flipud(grd.transpose())
1683 
1684  # # Import pflexcy.so (cython compiled version of dumpgrid)
1685  try:
1686 
1687  from sdspflexcy import dumpdatagrid, dumpdepogrid
1688  print 'using pflexcy'
1689  except:
1690  #print """WARNING: Using PURE Python to readgrid, execution will be slow.
1691  # Try compiling the FortFlex module or the pflexcy module
1692  # for your machine. For more information see the
1693  # pflexible/f2py_build directory or use cython with pflexcy.pyx
1694  #"""
1695  dumpdatagrid = _dumpgrid
1696  dumpdepogrid = _dumpgrid
1697 
1698  dat_cnt = 0
1699  nage = 1
1700  # #datagrid=np.zeros((numxgrid,numygrid,numzgrid[nspec-1],nspec,nageclass),np.float)
1701  wetgrid = np.zeros((H.numxgrid, H.numygrid, H.numpointspec, 1), np.float)
1702  drygrid = np.zeros((H.numxgrid, H.numygrid, H.numpointspec, 1), np.float)
1703  datagrid = np.zeros((H.numxgrid, H.numygrid, H.numzgrid, H.numpointspec, nage), np.float)
1704  # f = file(filename,'rb')
1705  # print filename
1706  f2 = BinaryFile(filename, order='fortran')
1707  # read data:
1708  skip(4)
1709  itime = getbin('i')
1710 
1711  for na in range(nage):
1712 
1713  for ks in range(H.numpointspec):
1714  # Read Wet Depostion
1715  skip()
1716  cnt_i = getbin('i')
1717  dmp_i = getdump(cnt_i, 'i')
1718  skip()
1719  cnt_r = getbin('i')
1720  dmp_r = getdump(cnt_r)
1721  if dmp_r.any():
1722  # print dmp_r, dmp_i
1723  wetgrid = dumpdepogrid(dmp_i, cnt_r, dmp_r, wetgrid, ks, na, H.numxgrid, H.numygrid)
1724 
1725  # Read Dry Deposition
1726  skip()
1727  cnt_i = getbin('i')
1728  dmp_i = getdump(cnt_i, 'i')
1729  skip()
1730  cnt_r = getbin('i')
1731  dmp_r = getdump(cnt_r)
1732  if dmp_r.any():
1733  # print dmp_r, dmp_i
1734  drygrid = dumpdepogrid(dmp_i, cnt_r, dmp_r, drygrid, ks, na, H.numxgrid, H.numygrid)
1735 
1736  # Read Concentrations
1737  skip()
1738  cnt_i = getbin('i')
1739  dmp_i = getdump(cnt_i, 'i')
1740  skip()
1741  cnt_r = getbin('i')
1742  dmp_r = getdump(cnt_r)
1743  # print len(dmp_r),len(dmp_i)
1744  # print cnt_r,cnt_i
1745  # print dmp_i
1746  # print type(dmp_i),type(cnt_r),type(dmp_r),type(datagrid),type(ks),type(na)
1747  # print type(H.numxgrid),type(H.numygrid)
1748  datagrid = dumpdatagrid(dmp_i, cnt_r, dmp_r, datagrid, ks, na, H.numxgrid, H.numygrid)
1749 
1750  # G[H['ibtime']].append(concgrid)
1751  f2.close()
1752 
1753  return datagrid, wetgrid, drygrid, itime
1754 
1755 
flextest.flexread.pflexible.BinaryFile
HELPER FUNCTIONS ##########.
Definition: pflexible.py:4484

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible._readV6 (   bf,
  h 
)
private 
Version 6, 7X read routines...

Definition at line 1275 of file pflexible.py.

1275 
1276 def _readV6(bf, h):
1277  """
1278  Version 6, 7X read routines...
1279 
1280  """
1281  # Utility functions
1282  skip = lambda n = 8 : bf.seek(n, 1)
1283  getbin = lambda dtype, n = 1 : bf.read(dtype, (n,))
1284 
1285  if bf:
1286  # bf.read('i')
1287  print h['numpoint']
1288  for i in range(h['numpoint']):
1289  # r2=getbin('i')
1290  i1 = getbin('i')
1291  i2 = getbin('i')
1292  # h['ireleasestart'].append( ss_start + datetime.timedelta(seconds=float(i1)) )
1293  # h['ireleaseend'].append( ss_start + datetime.timedelta(seconds=float(i2)) )
1294  h['ireleasestart'].append(i1)
1295  h['ireleaseend'].append(i2)
1296  h['kindz'][i] = getbin('i') # This is an int16, might need to to change something
1297  skip() # get xp, yp,...
1298 
1299  h['xp1'][i] = getbin('f')
1300  h['yp1'][i] = getbin('f')
1301  h['xp2'][i] = getbin('f')
1302  h['yp2'][i] = getbin('f')
1303  h['zpoint1'][i] = getbin('f')
1304  h['zpoint2'][i] = getbin('f')
1305 
1306  skip() # get n/mpart
1307  h['npart'][i] = getbin('i')
1308  h['mpart'][i] = getbin('i')
1309 
1310  getbin('i')
1311  # initialise release fields
1312 
1313  l = getbin('i') # get compoint length?
1314  gt = bf.tell() + l # create 'goto' point
1315  sp = ''
1316  sp = sp.join(getbin('c', l))
1317  bf.seek(gt) # skip ahead to gt point
1318 
1319  h['compoint'].append(sp) # species names in dictionary for each nspec
1320  # h['compoint'].append(''.join([getbin('c') for i in range(45)]))
1321  skip()
1322  # r1=getbin('i')
1323 
1324  # now loop for nspec to get xmass
1325  for v in range(h['nspec']):
1326  Dfmt = ['i', 'f', '2i', 'f', '2i', 'f', 'i']
1327  a = [bf.read(fmt) for fmt in Dfmt]
1328  h['xmass'][i, v] = a[1]
1329 
1330 
1331  return
1332 

Here is the caller graph for this function:

def flextest.flexread.pflexible._shout (   string,
  verbose = True 
)
private 
write a string to stdout if 'verbose' flag given

Definition at line 5059 of file pflexible.py.

5059 
5060 def _shout(string, verbose=True):
5061  """ write a string to stdout if 'verbose' flag given """
5062  w = sys.stdout.write
5063  if string[-1] != '\n':
5064  string = string + '\n'
5065  if verbose:
5066  w(string)
5067 

Here is the caller graph for this function:

def flextest.flexread.pflexible.add_nilu_logo (   fig,
  xo = 10,
  yo = 520,
  origin = 'upper' 
) 
Adds NILU logo to background of plot

Definition at line 4915 of file pflexible.py.

4915 
4916 def add_nilu_logo(fig, xo=10, yo=520, origin='upper'):
4917  """ Adds NILU logo to background of plot """
4918  im = image.imread('/xnilu_wrk/flex_wrk/WEB_APS/imgs/logo_nilu.png')
4919  im[:, :, -1] = 0.5
4920  fig.figimage(im, xo, yo, origin=origin)
4921  return fig
def flextest.flexread.pflexible.closest (   num,
  numlist 
) 
returns the index of the *closest* value in a list

Definition at line 4648 of file pflexible.py.

4648 
4649 def closest(num, numlist):
4650  """ returns the index of the *closest* value in a list """
4651  # check if we're using datetimes
4652  dates = False
4653  if isinstance(num, datetime.datetime):
4654  dates = True
4655  if dates:
4656  num = date2num(num)
4657  assert isinstance(numlist[0], datetime.datetime), \
4658  "num is date, numlist must be a list of dates"
4659  numlist = date2num(numlist)
4660 
4661  return (np.abs(numlist - num)).argmin()
4662 
4663 

Here is the caller graph for this function:

def flextest.flexread.pflexible.curtain_agltoasl (   H,
  curtain_agl,
  coords,
  below_gl = 0.0 
) 
converts the agl curtain to asl

    adds .asl_axis attribute to H

Definition at line 603 of file pflexible.py.

603 
604 def curtain_agltoasl(H, curtain_agl, coords, below_gl=0.0):
605  """ converts the agl curtain to asl
606 
607  adds .asl_axis attribute to H
608  """
609 
610  gl = groundlevel_for_line(H, H.longitude, H.latitude, coords)
611 
612  H.asl_axis = np.linspace(0, H.outheight[-1])
613 
614  xp = H.outheight - H.outheight[0]
615 
616  casl = np.zeros((len(H.asl_axis), len(coords)))
617 
618  for i in xrange(len(coords)):
619  casl[:, i] = np.interp(H.asl_axis, xp + gl[i], \
620  curtain_agl[:, i], left=below_gl)
621 
622  return casl
623 

Here is the call graph for this function:

def flextest.flexread.pflexible.curtain_for_flightrack (   H,
  flighttrack,
  nspec = 0,
  npspec_int = 0,
  index = 0,
  get_track = False 
) 
extracts curtain data along a given flight track

input:  H (a pf.Header with a FD data object from H.read_grid)
        if the read_grid method has not been called, we'll call
        it below.

        flighttrack
        (a numpy array with datetime, lon, lat, elv values)
        
        nspec: is optional for reference to the FD dict.
        
        index: is optional in case numpointspec > 1
        
        get_track: if True extracts the points along the flight
        track rather than the curtain.

output: curtain (a 2-d array with shape (len(flighttrack), len(H.outheight)
        
TODO::
    add interpolation, extract a track, not curtain (e.g. z points)

Definition at line 457 of file pflexible.py.

457 
458 def curtain_for_flightrack(H, flighttrack, nspec=0, npspec_int=0, index=0, get_track=False):
459  """
460  extracts curtain data along a given flight track
461 
462  input: H (a pf.Header with a FD data object from H.read_grid)
463  if the read_grid method has not been called, we'll call
464  it below.
465 
466  flighttrack
467  (a numpy array with datetime, lon, lat, elv values)
468 
469  nspec: is optional for reference to the FD dict.
470 
471  index: is optional in case numpointspec > 1
472 
473  get_track: if True extracts the points along the flight
474  track rather than the curtain.
475 
476  output: curtain (a 2-d array with shape (len(flighttrack), len(H.outheight)
477 
478  TODO::
479  add interpolation, extract a track, not curtain (e.g. z points)
480 
481  """
482 
483  # assert(isinstance(H, Header), 'ERROR: H Wrong Data Type')
484  if H.direction == 'forward':
485  if 'FD' in H.keys():
486  pass
487  else:
488  # need to call H.read_grid to fill read the grids
489  flighttimes = flighttrack[:, 0]
490  time_ret = H.closest_dates(flighttimes, take_set=True)
491  H.read_grid(nspec_ret=nspec, time_ret=time_ret,
492  npspec_int=npspec_int)
493 
494  elif H.direction == 'backward':
495  if 'C' in H.keys():
496  pass
497  else:
498  H.fill_backward()
499 
500  if get_track:
501  curtain = np.zeros(len(flighttrack))
502  else:
503  curtain = np.zeros((len(flighttrack), len(H.outheight)))
504 
505  for i, (t, lon, lat, elv) in enumerate(flighttrack):
506 
507 
508  # loop over FD tuples using t to extract the right time.
509  # assuming t is a datetime
510  idx = H.closest_date(t)
511  datestring = H.available_dates[idx]
512  if H.direction == 'forward':
513  grid = H.FD[(nspec, datestring)]['grid']
514  grid = grid[:, :, :, index]
515 
516  elif H.direction == 'backward':
517  grid = H.C[(0, idx)]['grid']
518 
519  I = closest(lon, H.longitude)
520  J = closest(lat, H.latitude)
521  if get_track:
522  K = closest(elv, H.outheight)
523  # print lon, I, H.longitude[I], '===========' , lat, J, H.latitude[J]
524  curtain[i] = grid[I, J, K]
525  else:
526  curtain[i] = grid[I, J, :]
527 
528  return curtain.T

Here is the call graph for this function:

def flextest.flexread.pflexible.curtain_for_line (   grid,
  X,
  Y,
  coords,
  index = 0 
) 
extracts curtain data from a grid given pairs of lon, lat

input:  H (a pf.Header with a FD data object from H.read_grid)
        if the read_grid method has not been called, we'll call
        it below.

        flighttrack
        (a numpy array with datetime, lon, lat, elv values)
        
        nspec: is optional for reference to the FD dict.
        
        index: is optional in case numpointspec > 1
        
        get_track: if True extracts the points along the flight
        track rather than the curtain.

output: curtain (a 2-d array with shape (len(flighttrack), len(H.outheight)
        
TODO::
    add interpolation, extract a track, not curtain (e.g. z points)

Definition at line 529 of file pflexible.py.

529 
530 def curtain_for_line(grid, X, Y, coords, index=0):
531  """
532  extracts curtain data from a grid given pairs of lon, lat
533 
534  input: H (a pf.Header with a FD data object from H.read_grid)
535  if the read_grid method has not been called, we'll call
536  it below.
537 
538  flighttrack
539  (a numpy array with datetime, lon, lat, elv values)
540 
541  nspec: is optional for reference to the FD dict.
542 
543  index: is optional in case numpointspec > 1
544 
545  get_track: if True extracts the points along the flight
546  track rather than the curtain.
547 
548  output: curtain (a 2-d array with shape (len(flighttrack), len(H.outheight)
549 
550  TODO::
551  add interpolation, extract a track, not curtain (e.g. z points)
552 
553  """
554  if len(grid.shape) == 4:
555  grid = grid[:, :, :, index]
556  else:
557  grid = grid
558 
559  curtain = np.zeros((len(coords), grid.shape[2]))
560  for i, (x, y) in enumerate(coords):
561 
562  I = closest(x, X)
563  J = closest(y, Y)
564  curtain[i] = grid[I, J, :]
565  curtain = np.nan_to_num(curtain)
566  # curtain = np.ma.fix_invalid(np.array(curtain,dtype=float))
567  return curtain.T
568 

Here is the call graph for this function:

def flextest.flexread.pflexible.data_range (   data,
  min = 'median' 
) 
return a data range for flexpart data

optional keyword min = ['median', 'mean', 'min']

Definition at line 4851 of file pflexible.py.

4851 
4852 def data_range(data, min='median'):
4853  """
4854  return a data range for flexpart data
4855 
4856  optional keyword min = ['median', 'mean', 'min']
4857  """
4858  dmax = np.nanmax(data)
4859  if np.isnan(dmax):
4860  dmax = 1e5
4861 
4862  if min == 'mean':
4863  dmin = np.mean(data[data.nonzero()])
4864  elif min == 'median':
4865  dmin = np.median(data[data.nonzero()])
4866  else:
4867  dmin = np.nanmin(data[data.nonzero()])
4868 
4869  if np.isnan(dmin):
4870  dmin = 1e-5
4871 
return [dmin, dmax]
def flextest.flexread.pflexible.fill_grids (   H,
  nspec = 0,
  FD = None,
  add_attributes = False 
) 
for backward runs, calculates the 20-day sensitivity at each release point.

Usage::

    > C = fill_backward(H,nspec=(0))


This will cycle through all available_dates and create the filled backward array
for each k in H.numpointspec.

Returns

    A dictionary keyed by a (species,k) tuple.
    OR
    C & FD attributes on H

Each element in the dictionary is a 3D array (x,y,z) for each species,k

.. note::
    USE with *caution*, this is **MEMORY** intensive!

Arguments


==============        ========================================
keyword               Description [default]
==============        ========================================
nspec                 the specied ID or a tuple of species IDs
FD                    FD can be passed if it is already read
add_attributes        will add C and FD as attributes to H,
                      rather than returning just C
==============        ========================================

.. todo::
    There's a lot of redundancy in the storage of attributes, maybe there is a
    better way to handle this.

Definition at line 2415 of file pflexible.py.

2415 
2416 def fill_grids(H, nspec=0, FD=None, add_attributes=False):
2417  """ for backward runs, calculates the 20-day sensitivity at each release point.
2418 
2419  Usage::
2420 
2421  > C = fill_backward(H,nspec=(0))
2422 
2423 
2424  This will cycle through all available_dates and create the filled backward array
2425  for each k in H.numpointspec.
2426 
2427  Returns
2428 
2429  A dictionary keyed by a (species,k) tuple.
2430  OR
2431  C & FD attributes on H
2432 
2433  Each element in the dictionary is a 3D array (x,y,z) for each species,k
2434 
2435  .. note::
2436  USE with *caution*, this is **MEMORY** intensive!
2437 
2438  Arguments
2439 
2440 
2441  ============== ========================================
2442  keyword Description [default]
2443  ============== ========================================
2444  nspec the specied ID or a tuple of species IDs
2445  FD FD can be passed if it is already read
2446  add_attributes will add C and FD as attributes to H,
2447  rather than returning just C
2448  ============== ========================================
2449 
2450  .. todo::
2451  There's a lot of redundancy in the storage of attributes, maybe there is a
2452  better way to handle this.
2453 
2454  """
2455 
2456 # assert H.direction == 'backward', "fill_backward is only valid for backward runs"
2457  # # make sure npsec is iterable
2458  if isinstance(nspec, int):
2459  species = [nspec]
2460  else:
2461  species = nspec
2462  assert iter(species), 'nspec must be iterable, or you can pass an int'
2463  # # initialize variables
2464  if FD is None:
2465  # # then we need to read the grids
2466  FD = read_grid(H, time_ret= -1, nspec_ret=species)
2467 
2468  C = Structure()
2469 
2470  if H.direction == 'backward':
2471 
2472  for s, k in itertools.product(species, range(H.numpointspec)):
2473  C[(s, k)] = Structure()
2474  C[(s, k)].grid = np.zeros((H.numxgrid, H.numygrid, H.numzgrid))
2475  C[(s, k)]['itime'] = None
2476  C[(s, k)]['timestamp'] = H.releasetimes[k]
2477  C[(s, k)]['species'] = H['species'][s]
2478  C[(s, k)]['gridfile'] = 'multiple'
2479  C[(s, k)]['rel_i'] = k
2480  C[(s, k)]['spec_i'] = s
2481 
2482  # read data grids and attribute/sum sensitivity
2483  print species
2484  for s in species:
2485 
2486  for d in FD.grid_dates:
2487  # cycle through all the date grids (20days back)
2488  for k in range(H.numpointspec):
2489  # cycle through each release point
2490  contribution = FD[(s, d)].grid[:, :, :, k]
2491  C[(s, k)].grid = C[(s, k)].grid + contribution
2492 
2493  # added this, not sure if it makes sense to have this for 'forward'
2494  # however, adding 'C' does not add to memory, as it points to same
2495  # memory location as FD
2496  if H.direction == 'forward':
2497  for s in species:
2498  for k in range(len(FD.grid_dates)):
2499  d = FD.grid_dates[k]
2500  C[(s, k)] = FD[(s, d)]
2501 
2502  for s, k in C:
2503  # # add total column
2504  C[(s, k)].slabs = get_slabs(H, C[(s, k)].grid)
2505  # # shape, min, max based on total column
2506  if H.direction == 'backward':
2507  C[(s, k)]['shape'] = C[(s, k)].grid[0].shape
2508  C[(s, k)]['max'] = C[(s, k)].grid[0].max()
2509  C[(s, k)]['min'] = C[(s, k)].grid[0].min()
2510 
2511 
2512  if add_attributes:
2513  H.C = C
2514  H.FD = FD
2515  else:
2516  return C
2517 
2518 

Here is the call graph for this function:

def flextest.flexread.pflexible.get_slabs (   H,
  G,
  index = None,
  normAreaHeight = True,
  scale = 1.0 
) 
Preps grid from readgridV8 for plotting.

Accepts an 3D or 4D GRID from readgrid with optional index.
shape ~= (360, 180, 3, 88, 1)

Usage::

    plotgrid = get_slabs(H,G)

Inputs

    H = A Header instance from a FLEXPART run
    G = A grid from the FLEXPARTDATA or FD dictionary returned from read_grid

Returns

    slabs : a dictionary of rank-2 arrays corresponding to vertical levels.

    **Note**: transposes the grid for plotting.
    Total Column as D[0]
    Level 1 as D[1]
    Level 2 as D[2]
    ...

Arguments

  .. tabularcolumns::  |l|L|

  ===============         ========================================
  keyword                 Description [default]
  ===============         ========================================
  index                   a release point index (k)
  normAreaHeight          [True], normalizes by the area/height
  ===============         ========================================

.. todo::
    Need to check the normalization and indexing.

Definition at line 2540 of file pflexible.py.

2540 
2541 def get_slabs(H, G, index=None, normAreaHeight=True, scale=1.0):
2542  """ Preps grid from readgridV8 for plotting.
2543 
2544  Accepts an 3D or 4D GRID from readgrid with optional index.
2545  shape ~= (360, 180, 3, 88, 1)
2546 
2547  Usage::
2548 
2549  plotgrid = get_slabs(H,G)
2550 
2551  Inputs
2552 
2553  H = A Header instance from a FLEXPART run
2554  G = A grid from the FLEXPARTDATA or FD dictionary returned from read_grid
2555 
2556  Returns
2557 
2558  slabs : a dictionary of rank-2 arrays corresponding to vertical levels.
2559 
2560  **Note**: transposes the grid for plotting.
2561  Total Column as D[0]
2562  Level 1 as D[1]
2563  Level 2 as D[2]
2564  ...
2565 
2566  Arguments
2567 
2568  .. tabularcolumns:: |l|L|
2569 
2570  =============== ========================================
2571  keyword Description [default]
2572  =============== ========================================
2573  index a release point index (k)
2574  normAreaHeight [True], normalizes by the area/height
2575  =============== ========================================
2576 
2577  .. todo::
2578  Need to check the normalization and indexing.
2579 
2580  """
2581  Heightnn = H['Heightnn']
2582  area = H['area']
2583  Slabs = Structure()
2584  grid_shape = G.shape
2585  if len(grid_shape) is 4:
2586  if index is None:
2587  if H.direction is 'forward':
2588  index = 0
2589  else:
2590  index = 0
2591  g = G[:, :, :, index]
2592  else:
2593  try:
2594  g = G[:, :, :, index]
2595  except:
2596  raise IOError( '######### ERROR: Which Release Point to get? ########## ')
2597 
2598  if len(grid_shape) is 3:
2599  if index is None:
2600  g = G
2601  else:
2602  g = G
2603 
2604 
2605  for i in range(g.shape[2]):
2606  # first time sum to create Total Column
2607  if i == 0:
2608  TC = np.sum(g, axis=2).T
2609  TC = TC * scale
2610  if normAreaHeight:
2611  data = g[:, :, i] / Heightnn[:, :, i]
2612  else:
2613  data = g[:, :, i]
2614 
2615  Slabs[i + 1] = data.T * scale
2616 
2617  Slabs[0] = TC
2618  return Slabs
2619 

Here is the caller graph for this function:

def flextest.flexread.pflexible.gridarea (   H) 
returns an array of area corresponding to each nx,ny,nz

Usage::

    > area = gridarea(H)


Returns
    OUT = array area corresponding to nx,ny,nz

Arguments
    H  = :class:`Header` object from readheader function.

Definition at line 746 of file pflexible.py.

746 
747 def gridarea(H):
748  """returns an array of area corresponding to each nx,ny,nz
749 
750  Usage::
751 
752  > area = gridarea(H)
753 
754 
755  Returns
756  OUT = array area corresponding to nx,ny,nz
757 
758  Arguments
759  H = :class:`Header` object from readheader function.
760 
761  """
762 
763 
764  pih = pi / 180.
765  r_earth = 6.371e6
766  cosfunc = lambda y : cos(y * pih) * r_earth
767 
768  nx = H['numxgrid']
769  ny = H['numygrid']
770  outlat0 = H['outlat0']
771  dyout = H['dyout']
772  dxout = H['dxout']
773  area = np.zeros((nx, ny))
774 
775  for iy in range(ny):
776  ylata = outlat0 + (float(iy) + 0.5) * dyout # NEED TO Check this, iy since arrays are 0-index
777  ylatp = ylata + 0.5 * dyout
778  ylatm = ylata - 0.5 * dyout
779  if (ylatm < 0 and ylatp > 0): hzone = dyout * r_earth * pih
780  else:
781  # cosfact = cosfunc(ylata)
782  cosfactp = cosfunc(ylatp)
783  cosfactm = cosfunc(ylatm)
784  if cosfactp < cosfactm:
785  hzone = sqrt(r_earth ** 2 - cosfactp ** 2) - sqrt(r_earth ** 2 - cosfactm ** 2)
786  else:
787  hzone = sqrt(r_earth ** 2 - cosfactm ** 2) - sqrt(r_earth ** 2 - cosfactp ** 2)
788 
789  gridarea = 2.*pi * r_earth * hzone * dxout / 360.
790  for ix in range(nx):
791  area[ix, iy] = gridarea
792 
793  return area

Here is the caller graph for this function:

def flextest.flexread.pflexible.groundlevel_for_line (   H,
  X,
  Y,
  coords,
  index = 0 
) 
extracts ground level from H.heightnn along a track of lon, lat

input:  H or H.heightnn (takes the lowest level)

        X, Y ,= H.longitude, H.latitude

        coords = zip(x, y) 
        
output: groundlevel (a 1-d array with shape (len(flighttrack)

Definition at line 569 of file pflexible.py.

569 
570 def groundlevel_for_line(H, X, Y, coords, index=0):
571  """
572  extracts ground level from H.heightnn along a track of lon, lat
573 
574  input: H or H.heightnn (takes the lowest level)
575 
576  X, Y ,= H.longitude, H.latitude
577 
578  coords = zip(x, y)
579 
580  output: groundlevel (a 1-d array with shape (len(flighttrack)
581 
582 
583  """
584  try:
585  hgt = H.Heightnn[:, :, 0]
586  except:
587  hgt = H
588 
589  # fix for hgt offset
590  hgt = hgt - hgt.min()
591 
592  grndlvl = np.zeros(len(coords))
593 
594  for i, (x, y) in enumerate(coords):
595 
596  I = closest(x, X)
597  J = closest(y, Y)
598  grndlvl[i] = hgt[I, J]
599  grndlvl = np.nan_to_num(grndlvl)
600 
601  return grndlvl - grndlvl.min()
602 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.main ( )

Below for use from command line ##########################################.

Definition at line 5141 of file pflexible.py.

5141 
5142 def main ():
5143 
5144  global options, args
5145  here = os.path.curdir
5146  dataDir = options.pathname
5147  H = read_header(dataDir, **{'readp':1}) # readheader and release pointspiu
5148  G = readgridV8(H, **{'unit':'time'})
5149  D = get_slabs(H, G)
5150 
5151  plt = mp.plot_map(D[0], H)
5152  plot_name = 'test'
5153  plt.title(plot_name, fontsize=10)
5154  plt.savefig(os.path.join(here, plot_name + '.png'))
5155  print plot_name
5156 
flextest.flexread.pflexible.main
def main
Below for use from command line ##########################################.
Definition: pflexible.py:5141

Here is the call graph for this function:

def flextest.flexread.pflexible.monthly_footprints (   H)

Definition at line 2401 of file pflexible.py.

2401 
2402 def monthly_footprints(H):
2403 
2404 
2405  footprints = np.zeros((H.ny, H.nx, len(H.C.keys())))
2406  for i, key in enumerate(H.C):
2407  footprints[:,:,i] = H.C[key].slabs[0]
2408 
2409  #footprints = np.average(footprints, axis=2)
2410 
2411  return footprints
2412 
2413 
2414 
def flextest.flexread.pflexible.plot_agespectra (   H,
  agespectra,
  plt_title = '',
  y_label = None,
  cbar_labels = None,
  FIGURE = None,
  y_datarange = None,
  web_version = False,
  continental = False,
  cum = False,
  bars = False,
  debug = False 
) 
plot an agespectra

Usage::

    > FIG = plot_agespectra(H,agespectra,*kwargs)

This creates an filled in between agespectra plot using information from
the header "H".

Returns
  A "mapping.py" ``FIGURE`` object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         ========================================
  keyword               Description [default]
  =============         ========================================
  runid                 an ID for the title
  yunits                units for the y-axis [None]
  FIGURE                A "FIGURE" object[None] (see mapping.py)
  data_range            y-axis data range
  web_version           Set to True if using the make_webpages
                        function. see: :func:`read_agespectrum`
  continental           Set to True if continental_spectrum
  cum                   Set to true if data should be cumulative
  bars                  Plot using stacked bars rather than a
                        filled line plot
  =============         ========================================

.. todo::
    There's a lot of redundancy in the storage of attributes, maybe there is a
    better way to handle this.

.. note::
    Required attributes of 'H' if you want to make a dummy H. Or just set
    H to "None" and it will be taken from the agespectra input.
    H.numageclasses
    H.releasetimes

Definition at line 2961 of file pflexible.py.

2961 
2962  bars=False, debug=False):
2963  """ plot an agespectra
2964 
2965  Usage::
2966 
2967  > FIG = plot_agespectra(H,agespectra,*kwargs)
2968 
2969  This creates an filled in between agespectra plot using information from
2970  the header "H".
2971 
2972  Returns
2973  A "mapping.py" ``FIGURE`` object.
2974 
2975  Arguments
2976 
2977  .. tabularcolumns:: |l|L|
2978 
2979  ============= ========================================
2980  keyword Description [default]
2981  ============= ========================================
2982  runid an ID for the title
2983  yunits units for the y-axis [None]
2984  FIGURE A "FIGURE" object[None] (see mapping.py)
2985  data_range y-axis data range
2986  web_version Set to True if using the make_webpages
2987  function. see: :func:`read_agespectrum`
2988  continental Set to True if continental_spectrum
2989  cum Set to true if data should be cumulative
2990  bars Plot using stacked bars rather than a
2991  filled line plot
2992  ============= ========================================
2993 
2994  .. todo::
2995  There's a lot of redundancy in the storage of attributes, maybe there is a
2996  better way to handle this.
2997 
2998  .. note::
2999  Required attributes of 'H' if you want to make a dummy H. Or just set
3000  H to "None" and it will be taken from the agespectra input.
3001  H.numageclasses
3002  H.releasetimes
3003 
3004 
3005  """
3006  # # make tick lables smaller
3007  mpl.rcParams['xtick.labelsize'] = 6
3008  mpl.rcParams['ytick.labelsize'] = 6
3009 
3010  if FIGURE == None:
3011  FIGURE = Structure()
3012  fig = plt.figure()
3013  FIGURE.fig = fig
3014  ax = fig.add_subplot(111)
3015  FIGURE.ax = ax
3016  else:
3017  fig = FIGURE.fig
3018  ax = FIGURE.ax
3019 
3020  if web_version:
3021  inspectra = agespectra.agespectrum[:, 4:]
3022  else:
3023  inspectra = agespectra[:]
3024 
3025  if continental:
3026  from nilu.pflexpart import emissions as em
3027  spectra_label = "Continents"
3028  spectra_type = "Continent"
3029  conts = em.Continents()
3030  cbar_labels = [i[1] for i in conts.flexpart_continents]
3031  else:
3032  spectra_label = "Ageclasses (Days)"
3033  spectra_type = "Ageclass"
3034 
3035  if not H:
3036  H = Structure()
3037  try:
3038  numageclasses = inspectra.numageclass
3039  releasetimes = inspectra.agespectrum[:, 0]
3040  except:
3041  numageclasses = inspectra.shape[1] - 1
3042  releasetimes = inspectra[:, 0]
3043  inspectra = inspectra[:, 1:]
3044  else:
3045  numageclasses = H.numageclasses
3046  releasetimes = H.releasetimes
3047 
3048  # if cum == 'norm':
3049  # ## Normalize the data so it fills to 100%
3050  # #spectra = np.zeros(inspectra.shape)
3051  # spectra = (inspectra.transpose()/np.sum(inspectra, axis=1)).transpose()
3052  # spectra = np.cumsum(spectra[:,:],axis=1)
3053  # #sums = np.sum(inspectra,axis=1)
3054  # #for i,elem in enumerate(inspectra):
3055  # # spectra[i,:] = elem/sums[i]
3056  # elif cum is True and cum != 'norm':
3057  # spectra = np.cumsum(inspectra[:,:],axis=1)
3058  # else:
3059  # spectra = inspectra
3060  if cum is not None:
3061  spectra = _cum_spec(inspectra, cum=cum)
3062 
3063  # Set up plotting environment colors
3064  Nc = np.array([float(i) / numageclasses for i in range(numageclasses)])
3065  norm = mpl.colors.normalize(Nc.min(), Nc.max())
3066  # jet = plt.cm.get_cmap('jet')
3067  jet = _gen_flexpart_colormap()
3068  plt.hold('on')
3069  facecolors = []
3070  # for i in range(0,H.numageclasses-1):
3071  for i in range(numageclasses):
3072  facecolors.append(jet(norm(Nc[i])))
3073  # ax.plot(ageclass[:,i])
3074  if i == 0:
3075  # create the baseline
3076  if bars:
3077  ax.bar(releasetimes, spectra[:, i], 0.03, color=facecolors[-1])
3078  else:
3079  ax.fill_between(releasetimes, np.zeros(len(spectra[:, i])), spectra[:, i],
3080  color=facecolors[-1], label='%s' % i)
3081  else:
3082  if bars:
3083  ax.bar(releasetimes, spectra[:, i], 0.03, color=facecolors[-1], bottom=spectra[:, i - 1])
3084  else:
3085  ax.fill_between(releasetimes, spectra[:, i - 1], spectra[:, i],
3086  color=facecolors[-1], label='%s' % (i))
3087  # facecolors.append(jet(norm(Nc[i+1])))
3088 
3089 
3090  # ax.set_yscale('log')
3091  if y_datarange:
3092  # print 'setting data range'
3093  ax.set_ylim(y_range)
3094  else:
3095  ax.set_ylim((0, spectra.max()))
3096  ax.grid(True)
3097  if y_label:
3098  ax.set_ylabel('%s' % y_label, rotation=90, ha='center', size='small')
3099 
3100  plt.title('%s' % (plt_title), size='small')
3101  fig.suptitle('Emissions by %s' % (spectra_type), size='small')
3102  fig.autofmt_xdate()
3103  # for xl in ax.get_xticklabels():
3104  # plt.setp(xl,size='x-small')
3105  # # ListedColormap
3106  pos = ax.get_position()
3107  l, b, w, h = getattr(pos, 'bounds', pos)
3108  # polygons = ax.collections
3109  # for i,p in enumerate(polygons): p.set_color(facecolors[i])
3110  # BoundaryNorm, and extended ends to show the "over" and "under"
3111  # value co
3112  ax2 = fig.add_axes([l, .08, w, 0.03])
3113  cmap = mpl.colors.ListedColormap(facecolors)
3114  # cmap.set_over('0.25')
3115  # cmap.set_under('0.75')
3116 
3117  # If a ListedColormap is used, the length of the bounds array must be
3118  # one greater than the length of the color list. The bounds must be
3119  # monotonically increasing.
3120  bounds = range(numageclasses + 1)
3121  # norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
3122  cb2 = mpl.colorbar.ColorbarBase(ax2,
3123  cmap=cmap,
3124  # norm=norm,
3125  # # to use 'extend', you must
3126  # # specify two extra boundaries:
3127  boundaries=bounds,
3128  # extend='both',
3129  # values=range(H.numageclasses+1),
3130  ticks=bounds, # optional
3131  spacing='proportional',
3132  orientation='horizontal')
3133  cb2.set_label(spectra_label, size='x-small')
3134  if cbar_labels:
3135  cb2.set_ticklabels(cbar_labels)
3136 
3137  if debug:
3138  plt.show()
3139 
3140  return FIGURE
3141 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.plot_at_level (   H,
  D,
  level = 1,
  ID = ' ',
  map_region = 5,
  projection = 'lcc',
  overlay = False,
  datainfo_str = None,
  log = True,
  data_range = None,
  coords = None,
  FIGURE = None,
  plot_title = None,
  units = None,
  kwargs 
) 
TODO:
-make units a function of H['species']

Definition at line 3625 of file pflexible.py.

3625 
3626  ):
3627  """
3628  TODO:
3629  -make units a function of H['species']
3630  """
3631  if units is None:
3632  units = H.output_unit
3633  if D is None:
3634  D = H.D
3635  rel_i = D.rel_i
3636  species = D.species
3637  timestamp = D.timestamp
3638  data = D.slabs[level]
3639 
3640  if level == 1:
3641  level_desc = 'Footprint'
3642  else:
3643  level_desc = H.outheight[level - 1]
3644 
3645  dmax = data.max()
3646  dmin = data.min()
3647  # print dmax,dmin
3648 
3649  if data_range == None:
3650  data_range = [dmin, dmax]
3651 
3652  if H.direction == 'backward' and H.options['readp']:
3653  zp1 = H['zpoint1'][rel_i]
3654  zp2 = H['zpoint2'][rel_i]
3655  if datainfo_str is None:
3656  # need to find a way to set m.a.s.l. or hPa here
3657  datainfo_str = """ Max Value: %.2g %s\n Release Z1: %.2f, Z2: %.2f (%s)\n""" \
3658  % (dmax, units, zp1, zp2, H.alt_unit)
3659  if plot_title is None:
3660  plot_title = """
3661  %s Sensitivity at %s %s: %s\n
3662  Release Start: %s, Release End: %s""" % \
3663  (ID, level_desc, H['alt_unit'], species, H['releasestart'][rel_i], H['releaseend'][rel_i])
3664  else:
3665  if datainfo_str is None:
3666  datainfo_str = """ Max Value: %.2g %s """ % (dmax, units)
3667  if plot_title is None:
3668  plot_title = """ %s Sensitivity at %s %s: %s \n %s """ % (ID, level_desc, H['alt_unit'], species, timestamp)
3669 
3670  FIGURE = plot_sensitivity(H, data, \
3671  data_range=data_range, \
3672  rel_i=rel_i, log=log,
3673  map_region=map_region, projection=projection,
3674  units=units, datainfo_str=datainfo_str,
3675  overlay=overlay,
3676  coords=coords, FIGURE=FIGURE, **kwargs)
3677 
3678 
3679  FIGURE.ax.set_title(plot_title, fontsize=10)
3680  return FIGURE

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_clusters (   H,
  T,
  rel_i = 0,
  ncluster = 0,
  sizescale = 10000,
  map_region = None,
  projection = None,
  coords = None,
  FIGURE = None,
  overlay = False,
  opacity = 0.7,
  draw_circles = True,
  draw_labels = True,
  MapPar = None 
) 
plots the clusters

Definition at line 3147 of file pflexible.py.

3147 
3148  MapPar=None):
3149  """ plots the clusters """
3150  trjs = T['Trajectories']
3151  labels = T['labels']
3152 
3153  # Set legend properties:
3154  p_legend = mpl.font_manager.FontProperties(size='8')
3155 
3156 
3157 
3158  # extract only releases of interest
3159  rel += 1 # account for zero indexing
3160  t = trjs[np.where(trjs[:, 0] == rel), :][0]
3161  if FIGURE == None:
3162  FIGURE = mp.get_FIGURE(map_region=map_region, projection=projection,
3163  coords=coords, MapPar=MapPar)
3164 
3165 
3166  try:
3167  ax = FIGURE.ax
3168  fig = FIGURE.fig
3169  m = FIGURE.m
3170  except:
3171  print 'problem getting ax,m, or fig.'
3172 
3173  # # Remove prior retroplume
3174  if not overlay:
3175  if len(ax.artists) > 1:
3176  try:
3177  art_ind = len(ax.artists)
3178  del ax.artists[:]
3179  # del previous texts, but not lat/lon labels
3180  del ax.texts[-art_ind:]
3181  except:
3182  'could not delete artists'
3183  pass
3184 
3185  # get cluster from data
3186  # set up cluster index, start column of T array
3187  cindx = [16, 20, 24, 28, 32, 36]
3188  if ncluster == 'all':
3189  # plot all clusters
3190  ncluster = range(5)
3191  elif isinstance(ncluster, int):
3192  ncluster = [ncluster]
3193 
3194  for nc in ncluster:
3195  # clstrindx = 16 + (5*(nc))
3196  clstrindx = cindx[nc]
3197  indx = [1] + range(clstrindx, clstrindx + 4)
3198  data = t[:, indx]
3199  # use ellipses
3200  ells = _genEllipse(data, m, sizescale=sizescale)
3201  texts = []
3202  for item in ells:
3203  e = item[0]
3204  xy = item[1]
3205  texts.append([xy, e.get_label()])
3206  e.set_clip_box(ax.bbox)
3207  e.set_alpha(opacity)
3208  if draw_circles:
3209  ax.add_artist(e)
3210  for txt in texts:
3211  x = txt[0][0]
3212  y = txt[0][1]
3213  # print x,y, t[1]
3214  if draw_labels:
3215  ax.text(x, y, txt[1])
3216  # plt.colorbar(ax)
3217  # ax.legend(ax.artists,(o.get_label() for o in ax.artists), prop=p_legend )
3218  FIGURE.ax = ax
3219  FIGURE.fig = fig
3220  FIGURE.m = m
3221 
3222  return FIGURE
3223 

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_curtain (   H,
  data,
  nx = None,
  ny = None,
  data_range = None,
  units = 'ppbv',
  datainfo_str = None,
  asl = True,
  plottitle = None,
  log = True,
  FIGURE = None,
  cax_title = None,
  method = 'contourf',
  figPar = None 
) 
plot_sensitivity: core function for plotting FLEXPART output.

Usage::

    > FIG = plot_sensitivity(H,data,*kwargs)

This returns the FIGURE object, and plots the sensitivity from the data contained in the "D"
array.

Inputs
  H = a :class:`Header` instance for a FLEXPART run.
  data = a 2d data array containing the sensitivity values to plot, this can be extracted from a
  grid instance (see :func:`readgridV8` and :func:`get_slabs`)

Returns
  A "mapping.py" ``FIGURE`` object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         ================================================
  keyword               Description [default]
  =============         ================================================
  data_range            range of data for scale bars, if None will
                        be taken from min/max of data
  asl                   [True] plot the units in m.a.s.l
  cax_title             string to be passed as colorbar title (units will
                        be passed to the format argument)
  units                 units for the scale bar
  datainfo_str          A string for labeling the scale bar.
  plottitle             Title for the plot.
  rel_i                 Release index to plot from the data array
  map_region                A map_region specified in mapping.py
  projection            [deprecated] use pre-defined map_regions.
  dropm                 Force creation of a new basemap instance
  coords                Used with autofit option. An array of lat,lon
                        values for the mapping module to autofit a
                        basemap instance to.
  autofit               Try to generate map_region automagically (flakey)
  overlay               Force removal of previous figure elements.
  transform             For use with imshow method, if your data is not
                        in same coordinates as projection, try to transform
                        the data to the basemap projection.
  log                   Create a logarithmic color scale.
  FIGURE                A FIGURE instance from mapping module get_FIGURE
  MapPar                A Structure of paramters to be passed to the
                        basemap class when creating an instance.
  method                The method to use for plotting array data. May be
                        one of: [pcolormesh], imshow, or contourf
  lsmask                set to True to draw a grey landseamask [False]
  =============         ================================================

.. todo::
    A lot!! There are some problems here and it is sensitive to options.
    lsmask = True seems to only work with certain projections (POLARCAT)

.. note::
    This is the primary main function for creating plots of flexpart output.
    Most the other routines are simply wrappers to this function, passing
    arguments in with some predefined settings. For more information on the
    mechanics of this function, see the mapping.py module and the matplotlib
    basemap toolkit.

Definition at line 4215 of file pflexible.py.

4215 
4216  figPar=None):
4217  """ plot_sensitivity: core function for plotting FLEXPART output.
4218 
4219  Usage::
4220 
4221  > FIG = plot_sensitivity(H,data,*kwargs)
4222 
4223  This returns the FIGURE object, and plots the sensitivity from the data contained in the "D"
4224  array.
4225 
4226  Inputs
4227  H = a :class:`Header` instance for a FLEXPART run.
4228  data = a 2d data array containing the sensitivity values to plot, this can be extracted from a
4229  grid instance (see :func:`readgridV8` and :func:`get_slabs`)
4230 
4231  Returns
4232  A "mapping.py" ``FIGURE`` object.
4233 
4234  Arguments
4235 
4236  .. tabularcolumns:: |l|L|
4237 
4238  ============= ================================================
4239  keyword Description [default]
4240  ============= ================================================
4241  data_range range of data for scale bars, if None will
4242  be taken from min/max of data
4243  asl [True] plot the units in m.a.s.l
4244  cax_title string to be passed as colorbar title (units will
4245  be passed to the format argument)
4246  units units for the scale bar
4247  datainfo_str A string for labeling the scale bar.
4248  plottitle Title for the plot.
4249  rel_i Release index to plot from the data array
4250  map_region A map_region specified in mapping.py
4251  projection [deprecated] use pre-defined map_regions.
4252  dropm Force creation of a new basemap instance
4253  coords Used with autofit option. An array of lat,lon
4254  values for the mapping module to autofit a
4255  basemap instance to.
4256  autofit Try to generate map_region automagically (flakey)
4257  overlay Force removal of previous figure elements.
4258  transform For use with imshow method, if your data is not
4259  in same coordinates as projection, try to transform
4260  the data to the basemap projection.
4261  log Create a logarithmic color scale.
4262  FIGURE A FIGURE instance from mapping module get_FIGURE
4263  MapPar A Structure of paramters to be passed to the
4264  basemap class when creating an instance.
4265  method The method to use for plotting array data. May be
4266  one of: [pcolormesh], imshow, or contourf
4267  lsmask set to True to draw a grey landseamask [False]
4268  ============= ================================================
4269 
4270  .. todo::
4271  A lot!! There are some problems here and it is sensitive to options.
4272  lsmask = True seems to only work with certain projections (POLARCAT)
4273 
4274  .. note::
4275  This is the primary main function for creating plots of flexpart output.
4276  Most the other routines are simply wrappers to this function, passing
4277  arguments in with some predefined settings. For more information on the
4278  mechanics of this function, see the mapping.py module and the matplotlib
4279  basemap toolkit.
4280 
4281 
4282  """
4283 
4284 
4285 
4286  methods = ['imshow', 'pcolormesh', 'contourf', 'contour', 'None']
4287  assert method in methods, "method keyword must be one of: %s" % methods
4288 
4289  if FIGURE is None:
4290  FIGURE = Structure()
4291 
4292 
4293  fig = plt.figure(**figPar)
4294  ax = fig.add_subplot(111)
4295 
4296  FIGURE['fig'] = fig
4297  FIGURE['ax'] = ax
4298 
4299 
4300  fig = FIGURE.fig
4301  ax = FIGURE.ax
4302 
4303  # # make the figure current
4304  plt.figure(fig.number)
4305  plt.axes(ax)
4306 
4307  # # get min/max range
4308  if data_range != None:
4309  dat_min = data_range[0]
4310  dat_max = data_range[1]
4311  else:
4312  dat_min, dat_max = data_range(data)
4313 
4314  if log:
4315  clevs = _log_clevs(dat_min, dat_max)
4316  else:
4317  clevs = [i for i in np.arange(dat_min, dat_max, (dat_max - dat_min) / 100)]
4318 
4319  # # Set up the IMAGE
4320  # # cmapnames = ['jet', 'hsv', 'gist_ncar', 'gist_rainbow', 'cool', 'spectral']
4321  colmap = _gen_flexpart_colormap()
4322  colmap.set_over(color='k', alpha=0.8)
4323  topodat = data
4324 
4325 
4326  if method == 'imshow':
4327  im = plt.imshow(np.flipud(topodat), cmap=colmap, zorder= -1,
4328  norm=mpl.colors.LogNorm(vmin=clevs[0],
4329  vmax=clevs[-1]))
4330 
4331  if method == 'pcolormesh':
4332  im = plt.pcolormesh(nx, ny, topodat, cmap=colmap,
4333  norm=mpl.colors.LogNorm(vmin=clevs[0],
4334  vmax=clevs[-1]))
4335  if method == 'contourf':
4336  im = plt.contourf(nx, ny, topodat, cmap=colmap, levels=clevs,
4337  norm=mpl.colors.LogNorm(vmin=clevs[0],
4338  vmax=clevs[-1]))
4339 
4340  if method == 'contour':
4341  im = plt.contour(nx, ny, topodat, cmap=colmap,
4342  norm=mpl.colors.LogNorm(vmin=clevs[0],
4343  vmax=clevs[-1]))
4344 
4345  # # Get the current axes, and properties for use later
4346  pos = ax.get_position()
4347  l, b, w, h = pos.bounds
4348 
4349  # # CREATE COLORBAR
4350  # # Note, with upgrades to matplotlib and basemap had to make some
4351  # # changes here... no more 'ghost' axes
4352  # # does a colorbar already exist?
4353  try:
4354  cb = FIGURE.cb
4355  cax = FIGURE.cax
4356  cb.update_normal(im)
4357  except:
4358  # # make a copy of the image object, change
4359  # # colormap to linear version of the precip colormap.
4360  # # create new axis for colorbar.
4361  h = 0.8 * h
4362  l = l + w + .02
4363  b = 0.5 - (h / 2)
4364  w = 0.025
4365  cax = plt.axes([l, b, w, h])
4366  # # using im2, not im (hack to prevent colors from being
4367  # # too compressed at the low end on the colorbar - results
4368  # # from highly nonuniform colormap)
4369  cb = fig.colorbar(im, cax=cax) # , format='%3.2g') # draw colorbar
4370  FIGURE.cax = cax
4371  FIGURE.cb = cb
4372 
4373 
4374  # # set colorbar label and ticks
4375  p_cax = mpl.font_manager.FontProperties(size='6')
4376  clabels = list(clevs[::10]) # #clevs, by 10 steps
4377  clabels.append(clevs[-1]) # # add the last label
4378  # cax.set_yticks(np.linspace(clabels[0],clabels[-1],len(clabels)))
4379  cax.set_yticks(np.linspace(0, 1, len(clabels)))
4380  cax.set_yticklabels(['%3.2g' % cl for cl in clabels])
4381  # fontproperties=p_cax)
4382 
4383  if cax_title:
4384  cax.set_title(cax_title.format(units), fontproperties=p_cax)
4385  else:
4386  cax.set_title('sensitivity\n({0})'.format(units),
4387  fontproperties=p_cax)
4388 
4389  # # make the original axes current again
4390  plt.axes(ax)
4391  plt.grid(True)
4392  FIGURE.ax = ax
4393  FIGURE.fig = fig
4394 
4395  if plottitle != None:
4396  FIGURE.ax.set_title(plottitle, fontsize=10)
4397 
4398  return FIGURE
4399 

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_markers (   H,
  lon,
  lat,
  zsize = None,
  zlevel = None,
  FIGURE = None,
  map_region = None,
  projection = None,
  coords = None,
  overlay = True,
  draw_circles = True,
  draw_labels = False,
  cbar2 = True,
  cbar2_title = None,
  MapPar = None,
  color = 'blue',
  edgecolor = 'none',
  zorder = 10,
  alpha = 0.85 
) 
Plot a group of x,y pairs, with optional size and color information.

Usage::

    > FIG = plot_trajectory(H,RT,rel_id,*kwargs)

.. note::
    You can substitude "None" for the :class:`Header` instance if you haven't
    created one


Returns
  A :mod:`mapping` "FIGURE" object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         =============================================
  keyword               Description [default]
  =============         =============================================
  rel_i                 **required** release index
  FIGURE                A "FIGURE" object[None] (see mapping.py)
  projection            A projection pre-defined in :mod:`mapping`
  coords                A set of lon,lat coords for autosetting
                        map map_region (not really working).
  overlay               [True] will overlay the trajectory 
                        on top of another map instance.
  draw_circles          [True] will mark the trajectory with
                        circles. If [False] nothing is drawn.
  draw_labels           [False] unimplemented
  cbar2                 [True] draws the scale bar as a second axis.
  cbar2_title            Optional argument to overide the cbar title.
  MapPar                A Structure of mapping parameters to pass
                        to the basemap instance if desired.
  =============         =============================================

.. todo::
    Who knows?! Could probably handle the overlaying better.

.. note::
    Just set H to "None" if you haven't already created a "Header" instance.

Definition at line 3305 of file pflexible.py.

3305 
3306  zorder=10, alpha=0.85):
3307 
3308  """Plot a group of x,y pairs, with optional size and color information.
3309 
3310  Usage::
3311 
3312  > FIG = plot_trajectory(H,RT,rel_id,*kwargs)
3313 
3314  .. note::
3315  You can substitude "None" for the :class:`Header` instance if you haven't
3316  created one
3317 
3318 
3319  Returns
3320  A :mod:`mapping` "FIGURE" object.
3321 
3322  Arguments
3323 
3324  .. tabularcolumns:: |l|L|
3325 
3326  ============= =============================================
3327  keyword Description [default]
3328  ============= =============================================
3329  rel_i **required** release index
3330  FIGURE A "FIGURE" object[None] (see mapping.py)
3331  projection A projection pre-defined in :mod:`mapping`
3332  coords A set of lon,lat coords for autosetting
3333  map map_region (not really working).
3334  overlay [True] will overlay the trajectory
3335  on top of another map instance.
3336  draw_circles [True] will mark the trajectory with
3337  circles. If [False] nothing is drawn.
3338  draw_labels [False] unimplemented
3339  cbar2 [True] draws the scale bar as a second axis.
3340  cbar2_title Optional argument to overide the cbar title.
3341  MapPar A Structure of mapping parameters to pass
3342  to the basemap instance if desired.
3343  ============= =============================================
3344 
3345  .. todo::
3346  Who knows?! Could probably handle the overlaying better.
3347 
3348  .. note::
3349  Just set H to "None" if you haven't already created a "Header" instance.
3350 
3351 
3352  """
3353 
3354  if H:
3355  pass
3356 
3357  # # Set up the FIGURE
3358  if FIGURE == None:
3359  FIGURE = mp.get_FIGURE(map_region=map_region, projection=projection,
3360  coords=coords, MapPar=MapPar)
3361  # #Get fig info and make active
3362  fig = FIGURE.fig
3363  m = FIGURE.m
3364  ax = FIGURE.fig.axes[0]
3365  plt.figure(fig.number)
3366  plt.axes(ax)
3367 
3368  # # prepare the data
3369  cx, cy = m(lon, lat)
3370  if zlevel is None:
3371  zlevel = np.ones(len(lon)) * 10.
3372  if zsize is None:
3373  zsize = np.ones(len(lon)) * 1
3374  marker = 'o'
3375 
3376 
3377  # # clear the previous track
3378  if 'circles' in FIGURE.keys():
3379  del FIGURE['circles']
3380  if overlay is False:
3381  del ax.collections[FIGURE.indices.collections:]
3382  del ax.texts[FIGURE.indices.texts:]
3383 
3384  # # plot the track
3385  if draw_circles:
3386  cmap = plt.get_cmap('jet')
3387  circles = m.scatter(cx, cy, zsize, zlevel, cmap=cmap,
3388  marker=marker, edgecolor=edgecolor,
3389  zorder=zorder, alpha=alpha)
3390 
3391  try:
3392  FIGURE.circles = circles
3393  except:
3394  pass
3395  # # make the figure active again
3396  plt.figure(fig.number)
3397  # # draw the legend and title
3398  # # CREATE COLORBAR
3399  # # does a colorbar already exist?
3400  # # Get the current axes, and properties for use later
3401  ax0 = fig.axes[0]
3402  pos = ax0.get_position()
3403  l, b, w, h = pos.bounds
3404  if cbar2:
3405  try:
3406  cb2 = FIGURE.cb2
3407  cax2 = FIGURE.cax2
3408  cb2.update_normal(circles)
3409  except:
3410  cax2 = plt.axes([l + w + 0.08, b, 0.02, h])
3411  # # using im2, not im (hack to prevent colors from being
3412  # # too compressed at the low end on the colorbar - results
3413  # # from highly nonuniform colormap)
3414  cb2 = fig.colorbar(circles, cax=cax2) # , format='%3.2g') # draw colorbar
3415  FIGURE.cax2 = cax2
3416  FIGURE.cb2 = cb2
3417  # # check if a colorbar legend exists (this will cause
3418  # # problems for subplot routines!)
3419  p_leg = mpl.font_manager.FontProperties(size='6')
3420  if cbar2_title:
3421  cax2.set_title(cbar2_title, fontproperties=p_leg)
3422  else:
3423  cax2.set_title('altitude\n(m)', fontproperties=p_leg)
3424  else:
3425  pass
3426  # cax2 = plt.axes([l+w+0.03, b, 0.025, h-0.2])
3427  # cb2 = fig.colorbar(jnkmap,cax=cax2) # draw colorbar
3428 
3429  # # delete the ghost instance
3430  # plt.close(jnkfig.number)
3431  # del jnkax, jnkfig, jnkmap
3432 
3433  plt.axes(ax)
3434  plt.setp(ax, xticks=[], yticks=[])
3435 
3436 
3437  FIGURE.fig = fig
3438  FIGURE.m = m
3439  FIGURE.ax = ax
3440  try:
3441  FIGURE.cax2 = cax2
3442  FIGURE.cb2 = cb2
3443  except:
3444  pass
3445  return FIGURE

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_METDATA (   METDATA,
  FIGURE,
  date = None,
  level = None 
) 
plots met data returned from :module:`pflexible.mapping.get_OpenDAP`

Definition at line 4400 of file pflexible.py.

4400 
4401 def plot_METDATA(METDATA, FIGURE, date=None, level=None):
4402  """ plots met data returned from :module:`pflexible.mapping.get_OpenDAP`
4403 
4404  """
4405 
4406 
4407  fig = FIGURE.fig
4408  m = FIGURE.m
4409  ax = FIGURE.ax
4410 
4411  # # make the figure current
4412  plt.figure(fig.number)
4413  plt.axes(ax)
4414 
4415  datelabels = METDATA['extracted_dates']
4416  slpin = METDATA['prmslmsl']['data'] * .01
4417  uin = METDATA['ugrdprs']['data']
4418  vin = METDATA['vgrdprs']['data']
4419  longitudes = METDATA['longitudes']
4420  latitudes = METDATA['latitudes']
4421 
4422  if date:
4423  cnt = datelabels.index(date)
4424  else:
4425  cnt = 0
4426  # # FOR OPENDAP OVERLAY
4427  slp_ = slpin[cnt, :, :]
4428  u_ = uin[cnt, :, :]
4429  v_ = vin[cnt, :, :]
4430 
4431  # plot wind vectors on projection grid (looks better).
4432  # first, shift grid so it goes from -180 to 180 (instead of 0 to 360
4433  # in longitude). Otherwise, interpolation is messed uplt.
4434 
4435  if longitudes[-1] - longitudes[0] < 360.:
4436  longarray = np.array(longitudes)
4437  slp, cyclon = addcyclic(slp_, longarray)
4438  u, cyclon = addcyclic(u_, longarray)
4439  v, cyclon = addcyclic(v_, longarray)
4440  ugrid, newlons = shiftgrid(180., u, cyclon, start=False)
4441  vgrid, newlons = shiftgrid(180., v, cyclon, start=False)
4442  slpgrid, newlons = shiftgrid(180., slp, cyclon, start=False)
4443  # transform vectors to projection grid.
4444 
4445  lons, lats = np.meshgrid(newlons, latitudes)
4446  # set desired contour levels.
4447  clevs = np.arange(960, 1061, 5)
4448  # compute native x,y coordinates of grid.
4449  x, y = m(lons, lats)
4450  # pos = FIG.ax.get_position()
4451  # l, b, w, h = pos.bounds
4452 
4453  CS = m.contour(x, y, slpgrid, clevs, linewidths=1, colors='k', animated=True)
4454  # CS = FIG.m.contourf(x,y,slpgrid,clevs,cmap=plt.cm.RdBu_r,animated=True,alpha=0.7)
4455  # CS = FIG.m.contour(x,y,slp[0,:,:],clevs,linewidths=0.5,colors='k',animated=True)
4456 # CS = FIG.m.contourf(x,y,slp[0,:,:],clevs,cmap=plt.cm.RdBu_r,animated=True,alpha=0.7)
4457  # plt.clabel(CS,inline=1,fontsize=10)
4458  # plot wind vectors over maplt.
4459  urot, vrot, xx, yy = m.transform_vector(ugrid, vgrid, newlons, latitudes, 51, 51,
4460  returnxy=True, masked=True)
4461  Q = m.quiver(xx, yy, urot, vrot, scale=500)
4462 # # make quiver key.
4463  qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W')
4464 
4465 
4466  FIGURE.ax = ax
4467  FIGURE.m = m
4468  FIGURE.fig = fig
4469 
4470  return FIGURE
4471 
4472 
4473 
4474 
4475 
4476 
4477 
4478 
4479 
4480 
4481 
def flextest.flexread.pflexible.plot_releases (   R,
  FIGURE = None,
  threedim = False,
  map_region = None,
  projection = None,
  coords = None,
  overlay = True,
  draw_circles = True,
  draw_labels = True,
  cbar2 = True,
  MapPar = None 
) 
plot the llon,llat,elv1 of the releases.

Usage::
    >F = plot_releases(R)

See the :func:`read_releases` function for information on how to generate "R"

Definition at line 2681 of file pflexible.py.

2681 
2682  MapPar=None):
2683  """ plot the llon,llat,elv1 of the releases.
2684 
2685  Usage::
2686  >F = plot_releases(R)
2687 
2688  See the :func:`read_releases` function for information on how to generate "R"
2689  """
2690 
2691  if threedim:
2692  try:
2693  from mpl_toolkits.mplot3d import Axes3D
2694  except:
2695  raise ImportError('mplot3d not available from mpl_toolkits!')
2696 
2697  fig = plt.figure()
2698  ax = Axes3D(fig)
2699 
2700  ax.plot(R.lllon, R.lllat, R.elv1)
2701 
2702  return fig
2703  else:
2704  # # Set up the FIGURE
2705  if FIGURE == None:
2706  FIGURE = mp.get_FIGURE(map_region=map_region, projection=projection,
2707  coords=coords, MapPar=MapPar)
2708  # #Get fig info and make active
2709  fig = FIGURE.fig
2710  m = FIGURE.m
2711  ax = FIGURE.ax
2712  plt.figure(fig.number)
2713  plt.axes(ax)
2714 
2715 
2716  # # prepare the data
2717  lon = R.lllon
2718  lat = R.lllat
2719  zlevel = R.elv1
2720  zsize = np.ones(len(lon)) * 30
2721  marker = 'o'
2722 
2723 
2724  # # clear the previous track
2725  if 'circles' in FIGURE.keys():
2726  del FIGURE['circles']
2727  if overlay is False:
2728  del ax.collections[FIGURE.indices.collections:]
2729  del ax.texts[FIGURE.indices.texts:]
2730 
2731  # # plot the track
2732  cx, cy = m(lon, lat)
2733  if draw_circles:
2734  cmap = plt.get_cmap('gist_gray')
2735  circles = m.scatter(cx, cy, zsize, zlevel, cmap=cmap,
2736  marker=marker, edgecolor=None,
2737  zorder=10, alpha=0.85)
2738  # # m.scatter has no color bar,
2739  # # so create a ghost 'scatter' instance:
2740  pos = ax.get_position()
2741  l, b, w, h = getattr(pos, 'bounds', pos)
2742  jnkfig = plt.figure()
2743  jnkax = jnkfig.add_axes([l, b, w, h], frameon=False)
2744  jnkmap = jnkax.scatter(cx, cy, zsize, zlevel,
2745  cmap=cmap, marker=marker,
2746  edgecolor=None)
2747 
2748  try:
2749  FIGURE.circles = circles
2750  except:
2751  pass
2752  # # make the figure active again
2753  plt.figure(fig.number);
2754  # # draw the legend and title
2755  # # check if a colorbar legend exists (this will cause
2756  # # problems for subplot routines!)
2757  if cbar2 is True:
2758  cax = plt.axes([l + w + 0.12, b, 0.02, h - 0.035])
2759  else:
2760  cax = plt.axes([l + w + 0.03, b, 0.025, h - 0.035])
2761  cb = fig.colorbar(jnkmap, cax=cax) # draw colorbar
2762  p_leg = mpl.font_manager.FontProperties(size='6')
2763  cax.set_title('altitude\n(m)', fontproperties=p_leg)
2764 
2765  # # delete the ghost instance
2766  plt.close(jnkfig.number)
2767  del jnkax, jnkfig, jnkmap
2768 
2769  plt.axes(ax)
2770  plt.setp(ax, xticks=[], yticks=[])
2771 
2772  if draw_labels:
2773  p_cax = mpl.font_manager.FontProperties(size='8',
2774  style='italic',
2775  weight='bold',
2776  )
2777 
2778  FIGURE.fig = fig
2779  FIGURE.m = m
2780  FIGURE.ax = ax
2781  return FIGURE
2782 

Here is the caller graph for this function:

def flextest.flexread.pflexible.plot_sensitivity (   H,
  data,
  data_range = None,
  units = 'ns m^2 / kg',
  datainfo_str = None,
  plottitle = None,
  rel_i = None,
  map_region = None,
  projection = None,
  dropm = None,
  coords = None,
  overlay = False,
  transform = True,
  log = True,
  FIGURE = None,
  MapPar = None,
  FigPar = None,
  cax_title = None,
  autofit = False,
  method = 'contourf',
  lsmask = False 
) 
plot_sensitivity: core function for plotting FLEXPART output.

Usage::

    > FIG = plot_sensitivity(H,data,*kwargs)

This returns the FIGURE object, and plots the sensitivity from the data contained in the "D"
array.

Inputs
  H = a :class:`Header` instance for a FLEXPART run.
  data = a 2d data array containing the sensitivity values to plot, this can be extracted from a
  grid instance (see :func:`readgridV8` and :func:`get_slabs`)

Returns
  A "mapping.py" ``FIGURE`` object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         ================================================
  keyword               Description [default]
  =============         ================================================
  data_range            range of data for scale bars, if None will
                        be taken from min/max of data
  cax_title             string to be passed as colorbar title (units will
                        be passed to the format argument)
  units                 units for the scale bar
  datainfo_str          A string for labeling the scale bar.
  plottitle             Title for the plot.
  rel_i                 Release index to plot from the data array
  map_region                A map_region specified in mapping.py
  projection            [deprecated] use pre-defined map_regions.
  dropm                 Force creation of a new basemap instance
  coords                Used with autofit option. An array of lat,lon
                        values for the mapping module to autofit a
                        basemap instance to.
  autofit               Try to generate map_region automagically (flakey)
  overlay               Force removal of previous figure elements.
  transform             For use with imshow method, if your data is not
                        in same coordinates as projection, try to transform
                        the data to the basemap projection.
  log                   Create a logarithmic color scale.
  FIGURE                A FIGURE instance from mapping module get_FIGURE
  MapPar                A Structure of paramters to be passed to the
                        basemap class when creating an instance.
  method                The method to use for plotting array data. May be
                        one of: [pcolormesh], imshow, or contourf
  lsmask                set to True to draw a grey landseamask [False]
  =============         ================================================

.. todo::
    A lot!! There are some problems here and it is sensitive to options.
    lsmask = True seems to only work with certain projections (POLARCAT)

.. note::
    This is the primary main function for creating plots of flexpart output.
    Most the other routines are simply wrappers to this function, passing
    arguments in with some predefined settings. For more information on the
    mechanics of this function, see the mapping.py module and the matplotlib
    basemap toolkit.

Definition at line 3894 of file pflexible.py.

3894 
3895  autofit=False, method='contourf', lsmask=False):
3896  # autofit=False,method='imshow'):
3897  """ plot_sensitivity: core function for plotting FLEXPART output.
3898 
3899  Usage::
3900 
3901  > FIG = plot_sensitivity(H,data,*kwargs)
3902 
3903  This returns the FIGURE object, and plots the sensitivity from the data contained in the "D"
3904  array.
3905 
3906  Inputs
3907  H = a :class:`Header` instance for a FLEXPART run.
3908  data = a 2d data array containing the sensitivity values to plot, this can be extracted from a
3909  grid instance (see :func:`readgridV8` and :func:`get_slabs`)
3910 
3911  Returns
3912  A "mapping.py" ``FIGURE`` object.
3913 
3914  Arguments
3915 
3916  .. tabularcolumns:: |l|L|
3917 
3918  ============= ================================================
3919  keyword Description [default]
3920  ============= ================================================
3921  data_range range of data for scale bars, if None will
3922  be taken from min/max of data
3923  cax_title string to be passed as colorbar title (units will
3924  be passed to the format argument)
3925  units units for the scale bar
3926  datainfo_str A string for labeling the scale bar.
3927  plottitle Title for the plot.
3928  rel_i Release index to plot from the data array
3929  map_region A map_region specified in mapping.py
3930  projection [deprecated] use pre-defined map_regions.
3931  dropm Force creation of a new basemap instance
3932  coords Used with autofit option. An array of lat,lon
3933  values for the mapping module to autofit a
3934  basemap instance to.
3935  autofit Try to generate map_region automagically (flakey)
3936  overlay Force removal of previous figure elements.
3937  transform For use with imshow method, if your data is not
3938  in same coordinates as projection, try to transform
3939  the data to the basemap projection.
3940  log Create a logarithmic color scale.
3941  FIGURE A FIGURE instance from mapping module get_FIGURE
3942  MapPar A Structure of paramters to be passed to the
3943  basemap class when creating an instance.
3944  method The method to use for plotting array data. May be
3945  one of: [pcolormesh], imshow, or contourf
3946  lsmask set to True to draw a grey landseamask [False]
3947  ============= ================================================
3948 
3949  .. todo::
3950  A lot!! There are some problems here and it is sensitive to options.
3951  lsmask = True seems to only work with certain projections (POLARCAT)
3952 
3953  .. note::
3954  This is the primary main function for creating plots of flexpart output.
3955  Most the other routines are simply wrappers to this function, passing
3956  arguments in with some predefined settings. For more information on the
3957  mechanics of this function, see the mapping.py module and the matplotlib
3958  basemap toolkit.
3959 
3960 
3961  """
3962 
3963 
3964 
3965  methods = ['imshow', 'pcolormesh', 'contourf', 'contour', 'None']
3966  assert method in methods, "method keyword must be one of: %s" % methods
3967 
3968  if FIGURE is None:
3969  if map_region is None:
3970  if MapPar is None:
3971  if autofit:
3972  MapPar = _gen_MapPar_fromHeader(H)
3973 
3974  FIGURE = mp.get_FIGURE(map_region=map_region,
3975  projection=projection, coords=coords,
3976  MapPar=MapPar, FigPar=FigPar)
3977  else:
3978  if FIGURE.m is None:
3979  FIGURE = mp.get_FIGURE(fig=FIGURE.fig, ax=FIGURE.ax, map_region=map_region,
3980  projection=projection, coords=coords,
3981  MapPar=MapPar, FigPar=FigPar)
3982 
3983  if overlay is False:
3984  del FIGURE.ax.images[FIGURE.indices.images:]
3985  del FIGURE.ax.collections[FIGURE.indices.collections:]
3986  del FIGURE.ax.lines[FIGURE.indices.lines:]
3987 
3988  if dropm != None:
3989  try:
3990  del m
3991  plt.close('all')
3992  except:
3993  print 'could not drop m'
3994 
3995  # # make tick lables smaller
3996  mpl.rcParams['xtick.labelsize'] = 6
3997  mpl.rcParams['ytick.labelsize'] = 6
3998 
3999  fig = FIGURE.fig
4000  m = FIGURE.m
4001  ax = FIGURE.ax
4002 
4003  # # make the figure current
4004  plt.figure(fig.number)
4005  plt.axes(ax)
4006 
4007  # # set up transformations for the data array
4008  if method == 'imshow':
4009  if m.projection not in ['cyl', 'merc', 'mill']:
4010  lats = np.arange(H.outlat0, (H.outlat0 + (H.numygrid * H.dyout)), H.dyout)[:-1]
4011  lons = np.arange(H.outlon0, (H.outlon0 + (H.numxgrid * H.dxout)), H.dxout)[:-1]
4012  data = data[:-1, :-1]
4013  else:
4014  lats = np.arange(H.outlat0, (H.outlat0 + (H.numygrid * H.dyout)), H.dyout)
4015  lons = np.arange(H.outlon0, (H.outlon0 + (H.numxgrid * H.dxout)), H.dxout)
4016 
4017  # # transform to nx x ny regularly spaced native projection grid
4018  if transform:
4019  dx = 2.*np.pi * m.rmajor / len(lons)
4020  nx = int((m.xmax - m.xmin) / dx) + 1; ny = int((m.ymax - m.ymin) / dx) + 1
4021  if nx is 1:
4022  topodat = data
4023  else:
4024  topodat = m.transform_scalar(data, lons, lats, nx, ny)
4025  else:
4026  topodat = data
4027 
4028  if method != 'imshow':
4029  # # Check to see if a cyclic wraparound is required
4030 
4031  lons = np.arange(H.outlon0, H.outlon0 + (H.dxout * H.numxgrid), H.dxout)
4032  lats = np.arange(H.outlat0, H.outlat0 + (H.dyout * H.numygrid), H.dyout)
4033  # # if required add polar coordinates
4034  if 'npstere' in m.projection:
4035  if lats[-1] != 90.:
4036  npole = np.ones(len(lons)).T * data[0, :]
4037  npole = np.reshape(npole, (1, -1))
4038  data = np.vstack((npole, data))
4039  lats = np.hstack((lats, [lats[-1] + H.dyout]))
4040 
4041  if m.projection != 'merc':
4042  if lons[-1] - lons[0] < 360.:
4043  topodat, lons = addcyclic(data, lons)
4044 
4045  if m.projection == 'merc':
4046  topodat = data
4047 
4048  # # get min/max range
4049  if data_range != None:
4050  dat_min = data_range[0]
4051  dat_max = data_range[1]
4052  else:
4053  dat_min, dat_max = data_range(data)
4054 
4055 
4056  if log:
4057  clevs = _log_clevs(dat_min, dat_max)
4058 
4059  else:
4060  clevs = [i for i in np.arange(dat_min, dat_max, (dat_max - dat_min) / 100)]
4061 
4062  # # draw land sea mask
4063  # m.fillcontinents(zorder=0)
4064  if lsmask:
4065  m.drawlsmask(ocean_color='grey', zorder= -10)
4066 
4067  # # Plot Release Location if points were read
4068  if H.options['readp']:
4069  if rel_i:
4070  releaselocation = (H.xpoint[rel_i], H.ypoint[rel_i])
4071  xpt, ypt = m(releaselocation[0], releaselocation[1])
4072  # # Remove prior location point
4073  try:
4074  del ax.lines[-1]
4075  except:
4076  pass
4077  location, = m.plot([xpt], [ypt], 'bx', linewidth=6, markersize=20, zorder=1000)
4078 
4079  # # Plot the footprint
4080 
4081  # # Set up the IMAGE
4082  # # cmapnames = ['jet', 'hsv', 'gist_ncar', 'gist_rainbow', 'cool', 'spectral']
4083  # colmap = plt.get_cmap('jet')
4084  colmap = _gen_flexpart_colormap()
4085  colmap.set_over(color='k', alpha=0.8)
4086  # # Plotting METHODS (pcolormesh now default, imshow is smoother)
4087  # print topodat.max(), topodat.min(), topodat.shape
4088  if method == 'imshow':
4089  im = m.imshow(topodat, cmap=colmap, zorder= -1,
4090  norm=mpl.colors.LogNorm(vmin=clevs[0],
4091  vmax=clevs[-1]))
4092 
4093  if method == 'pcolormesh':
4094  nx, ny = m(*np.meshgrid(lons, lats))
4095  im = m.pcolormesh(nx, ny, topodat, cmap=colmap,
4096  norm=mpl.colors.LogNorm(vmin=clevs[0],
4097  vmax=clevs[-1]))
4098  if method == 'contourf':
4099  # # Trying some fancier scaling
4100  # cnts,bins = np.histogram(topodat,bins=100)
4101  # topodat = np.ma.masked_where(topodat< .05* np.average((0,bins[1])),topodat)
4102  nx, ny = m(*np.meshgrid(lons, lats))
4103  im = m.contourf(nx, ny, topodat, cmap=colmap, levels=clevs,
4104  norm=mpl.colors.LogNorm(vmin=clevs[0],
4105  vmax=clevs[-1]))
4106 
4107  if method == 'contour':
4108  nx, ny = m(*np.meshgrid(lons, lats))
4109  im = m.contour(nx, ny, topodat, cmap=colmap,
4110  norm=mpl.colors.LogNorm(vmin=clevs[0],
4111  vmax=clevs[-1]))
4112 
4113  # # Get the current axes, and properties for use later
4114  pos = ax.get_position()
4115  l, b, w, h = pos.bounds
4116 
4117  # # CREATE COLORBAR
4118  # # Note, with upgrades to matplotlib and basemap had to make some
4119  # # changes here... no more 'ghost' axes
4120  # # does a colorbar already exist?
4121  try:
4122  cb = FIGURE.cb
4123  cax = FIGURE.cax
4124  cb.update_normal(im)
4125  except:
4126  # # make a copy of the image object, change
4127  # # colormap to linear version of the precip colormap.
4128  # pdb.set_trace()
4129  # im2 = copy.copy(im)
4130  # im2.set_cmap(colmap)
4131  # # create new axis for colorbar.
4132  h = 0.5 * h
4133  l = l + w + .03
4134  b = 0.5 - (h / 2)
4135  w = 0.025
4136  cax = plt.axes([l, b, w, h])
4137  # # using im2, not im (hack to prevent colors from being
4138  # # too compressed at the low end on the colorbar - results
4139  # # from highly nonuniform colormap)
4140  cb = fig.colorbar(im, cax=cax) # , format='%3.2g') # draw colorbar
4141  FIGURE.cax = cax
4142  FIGURE.cb = cb
4143  # cb.update_normal(im2)
4144 
4145 
4146 
4147 
4148  # # set colorbar label and ticks
4149  # pdb.set_trace()
4150  p_cax = mpl.font_manager.FontProperties(size='6')
4151  clabels = list(clevs[::10]) # #clevs, by 10 steps
4152  clabels.append(clevs[-1]) # # add the last label
4153  # cax.set_yticks(np.linspace(clabels[0],clabels[-1],len(clabels)))
4154  cax.set_yticks(np.linspace(0, 1, len(clabels)))
4155  cax.set_yticklabels(['%3.2g' % cl for cl in clabels])
4156  # fontproperties=p_cax)
4157  if H.direction == 'forward':
4158  cax.set_title('%s' % units,
4159  fontproperties=p_cax)
4160  else:
4161  if cax_title:
4162  cax.set_title(cax_title.format(units), fontproperties=p_cax)
4163  else:
4164  cax.set_title('sensitivity\n({0})'.format(units),
4165  fontproperties=p_cax)
4166 
4167  # # make the original axes current again
4168  plt.axes(ax)
4169 
4170  # # write text information block on plot
4171  # # first try to remove prior text by accessing
4172  # # the last text element added to the axes from
4173  # # the prior iteration.
4174  # # This is tricky when using together with plot_clusters...
4175  # # need to figure out how to resolve the indexing
4176  # # of what texts, collections, etc to delete, when iterating.
4177  try:
4178  del ax.texts[FIGURE.indices.texts:]
4179  del ax.artists[FIGURE.indices.artists:]
4180  except:
4181  pass
4182  if datainfo_str:
4183  plt.text(l, b + 1000,
4184  datainfo_str,
4185  fontsize=10,
4186  bbox=dict(boxstyle="round",
4187  ec=(1., 0.5, 0.5),
4188  fc=(1., 0.8, 0.8),
4189  alpha=0.8
4190  )
4191  )
4192 
4193  FIGURE.ax = ax
4194  FIGURE.m = m
4195  FIGURE.fig = fig
4196 
4197  if plottitle != None:
4198  # plt.title(plottitle,fontproperties=p_cax)
4199  # plt = plt
4200  FIGURE.ax.set_title(plottitle, fontsize=10)
4201  return FIGURE
4202 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.plot_sourcecontribution (   H,
  D,
  rel_i = None,
  s = None,
  ID = ' ',
  map_region = 'POLARCAT',
  data_range = None,
  coords = None,
  datainfo_str = None,
  units = None,
  FIGURE = None,
  overlay = False,
  cax_title = None,
  kwargs 
) 
plot_sourcecontribution: a wrapper around :func:`plot_sensitivity` for
plotting a source contribution field.

.. warning:
    This function is still quite finicky, I'm working on resolving it.

Usage::

    > FIG = plot_sourcecontribution(H,D,*kwargs)

Inputs
  H = a :class:`Header` instance from a FLEXPART run
  D = a grid instance from a FLEXDATA dictionary. See :func:`readgridV8`

  .. note::
      I am trying to create this so D can be a number of different types
      from a 2-d array to the grid instance type.


Returns
  A "mapping.py" ``FIGURE`` object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         ================================================
  keyword               Description [default]
  =============         ================================================
  ID                    Can be passed to prefix the title.
  data_range            range of data for scale bars, if None will
                        be taken from min/max of data
  datainfo_str          A string for labeling the scale bar.
  cax_title             A string to pass to plot_sensitivity for colorbar
                        title (units will be passed as format argument)
  k                     Release index to plot from the data array (aka rel_i)
  s                     Species index
  map_region                A map_region specified in mapping.py
  projection            [deprecated] use pre-defined map_regions.
  overlay               Force removal of previous figure elements.
  FIGURE                A FIGURE instance from mapping module get_FIGURE
  =============         ================================================

.. todo::
    A lot!! There are some problems here and it is sensitive to options.

.. note::
    k and rel_i are used throughout pflexible, should be more consistent in future.

Definition at line 3755 of file pflexible.py.

3755 
3756  **kwargs):
3757  """ plot_sourcecontribution: a wrapper around :func:`plot_sensitivity` for
3758  plotting a source contribution field.
3759 
3760  .. warning:
3761  This function is still quite finicky, I'm working on resolving it.
3762 
3763  Usage::
3764 
3765  > FIG = plot_sourcecontribution(H,D,*kwargs)
3766 
3767  Inputs
3768  H = a :class:`Header` instance from a FLEXPART run
3769  D = a grid instance from a FLEXDATA dictionary. See :func:`readgridV8`
3770 
3771  .. note::
3772  I am trying to create this so D can be a number of different types
3773  from a 2-d array to the grid instance type.
3774 
3775 
3776  Returns
3777  A "mapping.py" ``FIGURE`` object.
3778 
3779  Arguments
3780 
3781  .. tabularcolumns:: |l|L|
3782 
3783  ============= ================================================
3784  keyword Description [default]
3785  ============= ================================================
3786  ID Can be passed to prefix the title.
3787  data_range range of data for scale bars, if None will
3788  be taken from min/max of data
3789  datainfo_str A string for labeling the scale bar.
3790  cax_title A string to pass to plot_sensitivity for colorbar
3791  title (units will be passed as format argument)
3792  k Release index to plot from the data array (aka rel_i)
3793  s Species index
3794  map_region A map_region specified in mapping.py
3795  projection [deprecated] use pre-defined map_regions.
3796  overlay Force removal of previous figure elements.
3797  FIGURE A FIGURE instance from mapping module get_FIGURE
3798  ============= ================================================
3799 
3800  .. todo::
3801  A lot!! There are some problems here and it is sensitive to options.
3802 
3803  .. note::
3804  k and rel_i are used throughout pflexible, should be more consistent in future.
3805 
3806 
3807  """
3808  if D is None:
3809  D = H.D[(0, 0)]
3810  data = D.slabs[1]
3811  elif isinstance(D, list):
3812  data = D[0]
3813  elif isinstance(D, np.ndarray):
3814  data = D
3815  else:
3816  try:
3817  D = D[(s, rel_i)]
3818  data = D.slabs[0] # take total column
3819  except:
3820  raise IOError("Unrecognized 'grid' type passed")
3821 
3822  if s is None:
3823  try:
3824  species = D.species
3825  except:
3826  species = 'unknown'
3827  else:
3828  species = H.species[s]
3829 
3830  try:
3831  timestamp = D.timestamp
3832  except:
3833  if rel_i is None:
3834  try:
3835  timestamp = D.timestamp
3836  rel_i = D.rel_i
3837  except:
3838  raise IOError("Either a release index (rel_i) must be passed, \
3839  or the grid must have a timestamp attribute")
3840  else:
3841  timestamp = H.releasetimes[rel_i]
3842 
3843  if units is None:
3844  units = 'emission units'
3845 
3846  dmax = data.max()
3847  dmin = data.min()
3848 
3849  if data_range == None:
3850  data_range = [dmin, dmax]
3851 
3852  if H.direction == 'backward':
3853  zp1 = H['zpoint1'][rel_i]
3854  zp2 = H['zpoint2'][rel_i]
3855  if datainfo_str is None:
3856  datainfo_str = """ Max Value: %.2g %s\n Release Z1: %.2f , Z2: %.2f (%s)\n""" % \
3857  (dmax, units, zp1, zp2, H.alt_unit)
3858  plot_title = """
3859  %s Total Column Sensitivity: %s\n
3860  Release Start: %s, Release End: %s""" % \
3861  (ID, species, H['releasestart'][rel_i], H['releaseend'][rel_i])
3862  else:
3863  if datainfo_str is None:
3864  datainfo_str = """ Max Value: %.2g %s\n Release Z1: %.2f , Z2: %.2f (%s)\n""" % \
3865  (dmax, units, zp1, zp2, H.alt_unit)
3866  plot_title = """
3867  %s Total Column Sensitivity: %s\n %s""" % (ID, species, timestamp)
3868 
3869  FIGURE = plot_sensitivity(H, data, \
3870  data_range=data_range, \
3871  rel_i=rel_i, map_region=map_region, \
3872  units=units, \
3873  datainfo_str=datainfo_str, coords=coords, \
3874  FIGURE=FIGURE,
3875  cax_title=cax_title,
3876  overlay=overlay, **kwargs)
3877  FIGURE.ax.set_title(plot_title, fontsize=10)
3878  return FIGURE

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_spectra (   inspectra,
  plt_title = '',
  fig_title = '',
  y_label = None,
  spectra_label = 'bins',
  FIGURE = None,
  y_datarange = None,
  cum = False,
  labels = [],
  bars = False,
  debug = False 
) 
plot a spectra

Usage::

    > FIG = plot_agespectra(spectra,*kwargs)


Returns
  A "mapping.py" ``FIGURE`` object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         ========================================
  keyword               Description [default]
  =============         ========================================
  runid                 an ID for the title
  yunits                units for the y-axis [None]
  FIGURE                A "FIGURE" object[None] (see mapping.py)
  data_range            y-axis data range
  web_version           Set to True if using the make_webpages
                        function. see: :func:`read_agespectrum`
  continental           Set to True if continental_spectrum
  cum                   Set to true if data should be cumulative
  bars                  Plot using stacked bars rather than a
                        filled line plot
  =============         ========================================

.. todo::
    There's a lot of redundancy in the storage of attributes, maybe there is a
    better way to handle this.

Definition at line 2789 of file pflexible.py.

2789 
2790  bars=False, debug=False):
2791  """ plot a spectra
2792 
2793  Usage::
2794 
2795  > FIG = plot_agespectra(spectra,*kwargs)
2796 
2797 
2798  Returns
2799  A "mapping.py" ``FIGURE`` object.
2800 
2801  Arguments
2802 
2803  .. tabularcolumns:: |l|L|
2804 
2805  ============= ========================================
2806  keyword Description [default]
2807  ============= ========================================
2808  runid an ID for the title
2809  yunits units for the y-axis [None]
2810  FIGURE A "FIGURE" object[None] (see mapping.py)
2811  data_range y-axis data range
2812  web_version Set to True if using the make_webpages
2813  function. see: :func:`read_agespectrum`
2814  continental Set to True if continental_spectrum
2815  cum Set to true if data should be cumulative
2816  bars Plot using stacked bars rather than a
2817  filled line plot
2818  ============= ========================================
2819 
2820  .. todo::
2821  There's a lot of redundancy in the storage of attributes, maybe there is a
2822  better way to handle this.
2823 
2824 
2825  """
2826  # # make tick lables smaller
2827  mpl.rcParams['xtick.labelsize'] = 6
2828  mpl.rcParams['ytick.labelsize'] = 6
2829 
2830  if FIGURE == None:
2831  FIGURE = Structure()
2832  fig = plt.figure(figsize=(8, 6))
2833  FIGURE.fig = fig
2834  ax = fig.add_subplot(111) # ,pos=[0.1,0.2,.8,.7])
2835  FIGURE.ax = ax
2836  else:
2837  fig = FIGURE.fig
2838  ax = FIGURE.ax
2839 
2840 
2841  try:
2842  numageclasses = inspectra.numageclass
2843  releasetimes = inspectra.releasetimes
2844  inspectra = inspectra[:]
2845  except:
2846  numageclasses = inspectra.shape[1] - 1
2847  releasetimes = inspectra[:, 0]
2848  inspectra = inspectra[:, 1:]
2849 
2850  if cum == 'norm':
2851  # # Normalize the data so it fills to 100%
2852  # spectra = np.zeros(inspectra.shape)
2853  spectra = (inspectra.transpose() / np.sum(inspectra, axis=1)).transpose()
2854  spectra = np.cumsum(spectra[:, :], axis=1)
2855  # sums = np.sum(inspectra,axis=1)
2856  # for i,elem in enumerate(inspectra):
2857  # spectra[i,:] = elem/sums[i]
2858  elif cum is True and cum != 'norm':
2859  spectra = np.cumsum(inspectra[:, :], axis=1)
2860  else:
2861  spectra = inspectra
2862 
2863  # Set up plotting environment colors
2864  Nc = np.array([float(i) / numageclasses for i in range(numageclasses)])
2865  norm = mpl.colors.normalize(Nc.min(), Nc.max())
2866  # jet = plt.cm.get_cmap('jet')
2867  jet = _gen_flexpart_colormap()
2868  plt.hold('on')
2869  facecolors = []
2870  # for i in range(0,H.numageclasses-1):
2871  for i in range(numageclasses):
2872  facecolors.append(jet(norm(Nc[i])))
2873  if labels:
2874  lbl = labels[i]
2875  else:
2876  lbl = i
2877  # ax.plot(ageclass[:,i])
2878  if i == 0:
2879  # create the baseline
2880  if bars:
2881  ax.bar(releasetimes, spectra[:, i], 0.03, color=facecolors[-1])
2882  else:
2883  ax.fill_between(releasetimes, np.zeros(len(spectra[:, i])), spectra[:, i],
2884  color=facecolors[-1], label='%s' % lbl)
2885  else:
2886  if bars:
2887  ax.bar(releasetimes, spectra[:, i], 0.03, color=facecolors[-1], bottom=spectra[:, i - 1])
2888  else:
2889  ax.fill_between(releasetimes, spectra[:, i - 1], spectra[:, i],
2890  color=facecolors[-1], label='%s' % (lbl))
2891  # facecolors.append(jet(norm(Nc[i+1])))
2892 
2893 
2894  # ax.set_yscale('log')
2895  if y_datarange:
2896  print 'setting data range'
2897  ax.set_ylim(y_datarange)
2898  else:
2899  ax.set_ylim((0, spectra.max()))
2900  ax.grid(True)
2901  if y_label:
2902  ax.set_ylabel('%s' % y_label, rotation=90, ha='center', size='small')
2903 
2904  if plt_title:
2905  plt.title('%s' % (plt_title), size='small')
2906  if fig_title:
2907  fig.suptitle('Emissions by %s' % (spectra_type), size='small')
2908 
2909  fig.autofmt_xdate()
2910  # for xl in ax.get_xticklabels():
2911  # plt.setp(xl,size='x-small')
2912  # # ListedColormap
2913  pos = ax.get_position()
2914  l, b, w, h = getattr(pos, 'bounds', pos)
2915  # polygons = ax.collections
2916  # for i,p in enumerate(polygons): p.set_color(facecolors[i])
2917  # BoundaryNorm, and extended ends to show the "over" and "under"
2918  # value co
2919  ax2 = fig.add_axes([l, .08, w, 0.03])
2920  cmap = mpl.colors.ListedColormap(facecolors)
2921  # cmap.set_over('0.25')
2922  # cmap.set_under('0.75')
2923 
2924  # If a ListedColormap is used, the length of the bounds array must be
2925  # one greater than the length of the color list. The bounds must be
2926  # monotonically increasing.
2927  bounds = range(numageclasses + 1)
2928  # norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2929  cb2 = mpl.colorbar.ColorbarBase(ax2,
2930  cmap=cmap,
2931  # norm=norm,
2932  # # to use 'extend', you must
2933  # # specify two extra boundaries:
2934  boundaries=bounds,
2935  # extend='both',
2936  # values=range(H.numageclasses+1),
2937  ticks=bounds, # optional
2938  spacing='proportional',
2939  orientation='horizontal')
2940  if spectra_label:
2941  cb2.set_label(spectra_label, size='x-small')
2942  if labels:
2943  ax2.set_xticklabels(labels, va='center', ha='left', rotation=0, size='xx-small')
2944 
2945  if debug:
2946  plt.show()
2947 
2948  # make main axes current
2949  fig.add_subplot(ax)
2950 
2951  FIGURE.ax = ax
2952  FIGURE.fig = fig
2953 
2954  return FIGURE
2955 

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_totalcolumn (   H,
  D = None,
  ID = ' ',
  map_region = 5,
  projection = 'lcc',
  data_range = None,
  coords = None,
  FIGURE = None,
  overlay = False,
  datainfo_str = None,
  kwargs 
)

Definition at line 3688 of file pflexible.py.

3688 
3689  datainfo_str=None, **kwargs):
3690 
3691 
3692 
3693  if D is None:
3694  D = H.D[(0, 0)]
3695 
3696  if 'units' in kwargs:
3697  units = kwargs.pop('units')
3698  else:
3699  units = 'ns m kg-1'
3700 
3701  rel_i = D.rel_i
3702  species = D.species
3703  timestamp = D.timestamp
3704  data = D.slabs[0]
3705 
3706 
3707  if data_range == None:
3708  dmax = data.max()
3709  dmin = data.min()
3710  data_range = [dmin, dmax]
3711  else:
3712  dmin, dmax , = data_range
3713  # print dmin,dmax
3714 
3715 
3716  if H.direction == 'backward':
3717  rel_i = D.rel_i
3718  zp1 = H['zpoint1'][rel_i]
3719  zp2 = H['zpoint2'][rel_i]
3720  if datainfo_str is None:
3721  datainfo_str = """ Max Value: %.2g %s\n Release Z1: %.2f, Z2: %.2f (%s)\n""" % \
3722  (dmax, units, zp1, zp2, H.alt_unit)
3723  plot_title = """
3724  %s Total Column Sensitivity: %s\n
3725  Release Start: %s, Release End: %s""" % \
3726  (ID, species, H['releasestart'][rel_i], H['releaseend'][rel_i])
3727  else:
3728  if datainfo_str is None:
3729  datainfo_str = """ Max Value: %.2g %s""" % (dmax, units)
3730  plot_title = """
3731  %s Total Column Sensitivity: %s\n %s """ % (ID, species, timestamp)
3732 
3733 
3734  FIGURE = plot_sensitivity(H, data, \
3735  data_range=data_range, \
3736  rel_i=rel_i, map_region=map_region, \
3737  projection=projection, units=units, \
3738  datainfo_str=datainfo_str, coords=coords, \
3739  FIGURE=FIGURE, overlay=overlay, **kwargs)
3740 
3741  FIGURE.ax.set_title(plot_title, fontsize=10)
3742 
3743  return FIGURE
3744 

Here is the call graph for this function:

def flextest.flexread.pflexible.plot_trajchar (   H,
  R,
  numrelease = 1,
  varindex = [4,
  FIGURE = None,
  map_region = None,
  projection = None,
  coords = None 
) 
plot FLEXPART trajectory (or set of trajectories) characteristics
    R is a dictionary returned from pflexible.read_trajectories

    varindex is a list of the variables to plot from the trajectories

Definition at line 2637 of file pflexible.py.

2637 
2638  FIGURE=None, map_region=None, projection=None, coords=None):
2639  """ plot FLEXPART trajectory (or set of trajectories) characteristics
2640  R is a dictionary returned from pflexible.read_trajectories
2641 
2642  varindex is a list of the variables to plot from the trajectories
2643  """
2644  if FIGURE == None:
2645  FIGURE = mp.get_FIGURE(getm=False)
2646 
2647  try:
2648  ax = FIGURE.ax
2649  fig = FIGURE.fig
2650  m = FIGURE.m
2651  except:
2652  print 'problem getting ax,m, or fig.'
2653 
2654  try:
2655  del fig.axes[0]
2656  except:
2657  print 'ax not removed'
2658 
2659 
2660  T = R['Trajectories']
2661  labels = R['labels']
2662  t = T[np.where(T[:, 0] == numrelease), :][0]
2663  numplot = len(varindex)
2664 
2665  for i in range(numplot):
2666  sp = fig.add_subplot(numplot, 1, i + 1)
2667  sp.plot(t[:, 1], t[:, varindex[i]])
2668  ax = plt.gca()
2669  if i + 1 != numplot: plt.setp(ax, xticklabels=[])
2670  plt.ylabel(labels[varindex[i]])
2671  plt.grid('on')
2672 
2673  return FIGURE
2674 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.plot_trajectory (   H,
  T,
  rel_i,
  FIGURE = None,
  map_region = None,
  projection = None,
  coords = None,
  overlay = True,
  draw_circles = True,
  draw_labels = True,
  days_back = 20,
  cbar2 = True,
  cbar2_title = None,
  MapPar = None 
) 
Plot the center trajectory of the plume on the map

Usage::

    > FIG = plot_trajectory(H,RT,rel_id,*kwargs)

.. note::
    You can substitude "None" for the :class:`Header` instance if you haven't
    created one


Returns
  A :mod:`mapping` "FIGURE" object.

Arguments

  .. tabularcolumns::  |l|L|

  =============         =============================================
  keyword               Description [default]
  =============         =============================================
  rel_i                 **required** release index
  FIGURE                A "FIGURE" object[None] (see mapping.py)
  projection            A projection pre-defined in :mod:`mapping`
  coords                A set of lon,lat coords for autosetting
                        map map_region (not really working).
  overlay               [True] will overlay the trajectory 
                        on top of another map instance.
  draw_circles          [True] will mark the trajectory with
                        circles. If [False] nothing is drawn.
  draw_labels           [True] will draw the 'day' numbers on the
                        trajectories.
  days_back             For how many days back should the labels be
                        shown? [20]
  cbar2                 [True] draws the scale bar as a second axis.
  cbar2_title            Optional argument to overide the cbar title.
  MapPar                A Structure of mapping parameters to pass
                        to the basemap instance if desired.
  =============         =============================================

.. todo::
    Who knows?! Could probably handle the overlaying better.

.. note::
    Just set H to "None" if you haven't already created a "Header" instance.

Definition at line 3453 of file pflexible.py.

3453 
3454  MapPar=None):
3455  """Plot the center trajectory of the plume on the map
3456 
3457  Usage::
3458 
3459  > FIG = plot_trajectory(H,RT,rel_id,*kwargs)
3460 
3461  .. note::
3462  You can substitude "None" for the :class:`Header` instance if you haven't
3463  created one
3464 
3465 
3466  Returns
3467  A :mod:`mapping` "FIGURE" object.
3468 
3469  Arguments
3470 
3471  .. tabularcolumns:: |l|L|
3472 
3473  ============= =============================================
3474  keyword Description [default]
3475  ============= =============================================
3476  rel_i **required** release index
3477  FIGURE A "FIGURE" object[None] (see mapping.py)
3478  projection A projection pre-defined in :mod:`mapping`
3479  coords A set of lon,lat coords for autosetting
3480  map map_region (not really working).
3481  overlay [True] will overlay the trajectory
3482  on top of another map instance.
3483  draw_circles [True] will mark the trajectory with
3484  circles. If [False] nothing is drawn.
3485  draw_labels [True] will draw the 'day' numbers on the
3486  trajectories.
3487  days_back For how many days back should the labels be
3488  shown? [20]
3489  cbar2 [True] draws the scale bar as a second axis.
3490  cbar2_title Optional argument to overide the cbar title.
3491  MapPar A Structure of mapping parameters to pass
3492  to the basemap instance if desired.
3493  ============= =============================================
3494 
3495  .. todo::
3496  Who knows?! Could probably handle the overlaying better.
3497 
3498  .. note::
3499  Just set H to "None" if you haven't already created a "Header" instance.
3500 
3501 
3502  """
3503 
3504  if H:
3505  pass
3506 
3507  # # Set up the FIGURE
3508  if FIGURE == None:
3509  FIGURE = mp.get_FIGURE(map_region=map_region, projection=projection,
3510  coords=coords, MapPar=MapPar)
3511  # #Get fig info and make active
3512  fig = FIGURE.fig
3513  m = FIGURE.m
3514  ax = FIGURE.fig.axes[0]
3515  plt.figure(fig.number)
3516  plt.axes(ax)
3517 
3518  # # prepare the data
3519  trjs = T['Trajectories']
3520  rel = rel_i + 1 # # account for zero indexing
3521 
3522  # #extract only releases of interest
3523  t = trjs[np.where(trjs[:, 0] == rel), :][0]
3524 
3525  # # Get the data for the days_back we're interested in
3526  day_labels = _gen_daylabels(t[:days_back, 1])
3527  lon = t[:days_back, 2]
3528  lat = t[:days_back, 3]
3529  zlevel = t[:days_back, 4]
3530  zsize = np.ones(len(lon)) * 50
3531  marker = 'o'
3532 
3533 
3534  # # clear the previous track
3535  if 'circles' in FIGURE.keys():
3536  del FIGURE['circles']
3537  if overlay is False:
3538  del ax.collections[FIGURE.indices.collections:]
3539  del ax.texts[FIGURE.indices.texts:]
3540 
3541  # # plot the track
3542  cx, cy = m(lon, lat)
3543  if draw_circles:
3544  cmap = plt.get_cmap('gist_gray')
3545  circles = m.scatter(cx, cy, zsize, zlevel, cmap=cmap,
3546  marker=marker, edgecolor=None,
3547  zorder=10, alpha=0.85)
3548 
3549  try:
3550  FIGURE.circles = circles
3551  except:
3552  pass
3553  # # make the figure active again
3554  plt.figure(fig.number)
3555  # # draw the legend and title
3556  # # CREATE COLORBAR
3557  # # does a colorbar already exist?
3558  # # Get the current axes, and properties for use later
3559  ax0 = fig.axes[1]
3560  pos = ax0.get_position()
3561  l, b, w, h = pos.bounds
3562  if cbar2:
3563  try:
3564  cb2 = FIGURE.cb2
3565  cax2 = FIGURE.cax2
3566  cb2.update_normal(circles)
3567  except:
3568  cax2 = plt.axes([l + w + 0.08, b, 0.02, h])
3569  # # using im2, not im (hack to prevent colors from being
3570  # # too compressed at the low end on the colorbar - results
3571  # # from highly nonuniform colormap)
3572  cb2 = fig.colorbar(circles, cax=cax2) # , format='%3.2g') # draw colorbar
3573  FIGURE.cax2 = cax2
3574  FIGURE.cb2 = cb2
3575  # # check if a colorbar legend exists (this will cause
3576  # # problems for subplot routines!)
3577  else:
3578  cax2 = plt.axes([l + w + 0.03, b, 0.025, h - 0.2])
3579  cb2 = fig.colorbar(jnkmap, cax=cax2) # draw colorbar
3580 
3581  p_leg = mpl.font_manager.FontProperties(size='6')
3582  if cbar2_title:
3583  cax.set_title(cbar2_title, fontproperties=p_leg)
3584  else:
3585  cax2.set_title('altitude\n(m)', fontproperties=p_leg)
3586  # # delete the ghost instance
3587  # plt.close(jnkfig.number)
3588  # del jnkax, jnkfig, jnkmap
3589 
3590  plt.axes(ax)
3591  plt.setp(ax, xticks=[], yticks=[])
3592 
3593  if draw_labels:
3594  p_cax = mpl.font_manager.FontProperties(size='10',
3595  style='italic',
3596  weight='bold',
3597  )
3598 
3599 
3600 
3601  for i, (p, x, y) in enumerate(zip(day_labels, cx, cy)):
3602  if x > m.llcrnrx and x < m.urcrnrx:
3603  if y > m.llcrnry and y < m.urcrnry:
3604  ax.text(x, y + y * .02, '{0}'.format(p), va='bottom', ha='left',
3605  fontproperties=p_cax, zorder=11,
3606  color='white', bbox=dict(facecolor='green', alpha=0.5)
3607  )
3608 
3609  FIGURE.fig = fig
3610  FIGURE.m = m
3611  FIGURE.ax = ax
3612  FIGURE.cax2 = cax2
3613  FIGURE.cb2 = cb2
3614  return FIGURE
3615 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.plot_trajectory_ellipses (   H,
  T,
  rel_i = 0,
  ncluster = 0,
  sizescale = 10000,
  map_region = None,
  projection = None,
  coords = None,
  FIGURE = None,
  overlay = False,
  opacity = 0.7,
  draw_circles = True,
  draw_labels = True,
  MapPar = None 
) 
Plots trajectories from FLEXPART output.

Definition at line 3229 of file pflexible.py.

3229 
3230  MapPar=None):
3231  """
3232  Plots trajectories from FLEXPART output.
3233 
3234  """
3235 
3236  trjs = T['Trajectories']
3237  labels = T['labels']
3238 
3239  # Set legend properties:
3240  p_legend = mpl.font_manager.FontProperties(size='8')
3241 
3242  # extract only releases of interest according to rel_i
3243  rel = rel_i + 1 # account for zero indexing
3244  t = trjs[np.where(trjs[:, 0] == rel), :][0]
3245  if FIGURE == None:
3246  FIGURE = mp.get_FIGURE(map_region=map_region, projection=projection,
3247  coords=coords, MapPar=MapPar)
3248 
3249  try:
3250  ax = FIGURE.ax
3251  fig = FIGURE.fig
3252  m = FIGURE.m
3253  except:
3254  print 'problem getting ax,m, or fig.'
3255 
3256  # # Remove prior retroplume
3257  if overlay:
3258  if len(ax.artists) > 1:
3259  try:
3260  art_ind = len(ax.artists)
3261  del ax.artists[:]
3262  # del previous texts, but not lat/lon labels
3263  del ax.texts[-art_ind:]
3264  except:
3265  'could not delete artists'
3266  pass
3267 
3268  # get trajectory from data
3269  indx = [1, 2, 3, 4] # return, j,x,y,z
3270  data = t[:, indx]
3271  # use ellipses
3272  ells = _genEllipse(data, m, sizescale=sizescale)
3273  texts = []
3274  for item in ells:
3275  e = item[0]
3276  xy = item[1]
3277  texts.append([xy, e.get_label()])
3278  e.set_clip_box(ax.bbox)
3279  e.set_alpha(opacity)
3280  if draw_circles:
3281  ax.add_artist(e)
3282  for txt in texts:
3283  x = txt[0][0]
3284  y = txt[0][1]
3285  # print x,y, t[1]
3286  if draw_labels:
3287  ax.text(x, y, txt[1])
3288  plt.draw()
3289  FIGURE.ax = ax
3290  FIGURE.fig = fig
3291  FIGURE.m = m
3292  FIGURE.circles = ells
3293 
3294  return FIGURE

Here is the call graph for this function:

def flextest.flexread.pflexible.read_agespectrum (   filename,
  part = False,
  ndays = 20 
) 
Reads the spectrum.txt files generated from the "make_webpages" scripts.

.. note::
    This functionality depends on having run the make_webpages scripts
    in-house. It is not supported in the public API.

USAGE::

    > A = read_agespectrum(filename)

Returns a dictionary containing ageclass information

  .. tabularcolumns::  |l|L|

  =============       ========================================
  Keys                Description
  =============       ========================================
  agespectrum         2d-array
                      (dt, lat, lon, alt, tracer[:numageclass]
                      for each release)
  ageclasses          1d-array(ageclass ages)
  numpspec            number of species
  numageclass         number of ageclassesA
  part                if it is from make webpages without the
                      header information.
  ndays               The number of ageclasses
  =============       ========================================

Definition at line 624 of file pflexible.py.

624 
625 def read_agespectrum(filename, part=False, ndays=20):
626  """
627  Reads the spectrum.txt files generated from the "make_webpages" scripts.
628 
629  .. note::
630  This functionality depends on having run the make_webpages scripts
631  in-house. It is not supported in the public API.
632 
633  USAGE::
634 
635  > A = read_agespectrum(filename)
636 
637  Returns a dictionary containing ageclass information
638 
639  .. tabularcolumns:: |l|L|
640 
641  ============= ========================================
642  Keys Description
643  ============= ========================================
644  agespectrum 2d-array
645  (dt, lat, lon, alt, tracer[:numageclass]
646  for each release)
647  ageclasses 1d-array(ageclass ages)
648  numpspec number of species
649  numageclass number of ageclassesA
650  part if it is from make webpages without the
651  header information.
652  ndays The number of ageclasses
653  ============= ========================================
654 
655  """
656 
657  f = file(filename, 'r').readlines()
658 
659  line1 = f[0].strip().split()
660  if part:
661  numageclass = ndays
662  ageclasses = np.arange(1, ndays * 1) * 60 * 60 * 24
663  numspec = 1
664  else:
665  numageclass = int(line1[0])
666  ageclasses = np.array([int(i) for i in line1[1:]])
667  numspec = int(f[1].strip().split()[0])
668  A = Structure()
669  D = []
670  for line in f[2:]:
671  line = line.strip().split()
672  # convert ibdate, ibtime to datetime
673  dt = datetime.datetime.strptime(line[0] + line[1].zfill(6), '%Y%m%d%H%M%S')
674  data = [float(i) for i in line[2:]]
675  # reordering the array, dt object will be '0'
676  data = [data[1]] + [data[0]] + data[2:]
677  D.append([dt] + data)
678  D = np.array(D)
679  A['agespectrum'] = D
680  A['numageclass'] = numageclass
681  A['ageclasses'] = ageclasses
682  A['numspec'] = numspec
683  A['filename'] = filename
684  A['info'] = \
685  """
686  A dictionary containing ageclass information:
687  keys:
688  agespectrum = 2d-array(dt, lon, lat, alt, tracer[:numageclass] for each release)
689  ageclasses = 1d-array(ageclass ages)
690  numpspec = number of species
691  numageclass = number of ageclasses
692  info = this message
693  """
694  return A
695 

Here is the call graph for this function:

def flextest.flexread.pflexible.read_command (   path,
  headerrows = 7 
)

Functions for reading FLEXPART output #####. 

Reads a FLEXPART COMMAND file.

.. note::
    Assumes releases file has a header of 7 lines. Use
    option "headerrows" to override this.

USAGE::

    > A = read_command(filepath)

where filepath is either a file object or a path.


Returns
a dictionary with the following keys

================    =================================================================
fields              description
================    =================================================================
SIM_DIR             Simulation direction
SIM_START           Text str of YYYYMMDD HHMMSS
SIM_END             Text str of YYYYMMDD HHMMSS
AVG_CNC_INT         Average concentrations are calculated every SSSSS seconds
AVG_CNC_TAVG        The average concentrations are time averages of SSSSS sec
CNC_SAMP_TIME       The concentrations are sampled every SSSSS seconds to
                    calculate the time average concentration.
T_PARTSPLIT         Time constant for particle splitting.
SYNC                All processes are synchronized with this time interval
CTL                 --   
IFINE               IFINE=Reduction factor for time step used for vertical wind
IOUT                IOUT determines how the output shall be made: concentration
                    (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume
                    trajectory mode, or concentration + plume trajectory mode.
IPOUT               IPOUT determines whether particle positions are outputted
                    (in addition to the gridded concentrations or mixing ratios)
                    or not. 0=no output, 1 output every output interval, 2 only
                    at end of the simulation
LSUBGRID            Switch on/off subgridscale terrain parameterization
                    (increase of mixing heights due to subgridscale orog. var
LCONVECTION         Switch on/off the convection parameterization
LAGESPECTRA         Switch on/off the calculation of age spectra: if yes, the
                    file AGECLASSES must be available
IPIN                If IPIN=1, a file "partposit_end" from a previous run must
                    be available in the output directory. Particle positions
                    are read in and previous simulation is continued. If
                    IPIN=0, no particles from a previous run are used
IOFR                 Switch on/off writing out each release.
IFLUX               If IFLUX is set to 1, fluxes of each species through each
                    of the output boxes are calculated. Six fluxes,
                    corresponding to northward, southward, eastward, westward,
                    upward and downward are calculated for each grid cell of
                    the output grid. The control surfaces are placed in the
                    middle of each output grid cell. If IFLUX is set to 0,
                    no fluxes are determined.
MDOMAINFILL         If MDOMAINFILL is set to 1, the first box specified in file
                    RELEASES is used as the domain where domain-filling
                    trajectory calculations are to be done. Particles are
                    initialized uniformly distributed (according to the air mass
                    distribution) in that domain at the beginning of the
                    simulation, and are created at the boundaries throughout
                    the simulation period
IND_SOURCE          IND_SOURCE switches between different units for
                    concentrations at the source. NOTE that in backward
                    simulations the release of computational particles takes
                    place at the "receptor" and the sampling of particles at
                    the "source".  1=mass units (for bwd-runs = concentration)
                    2=mass mixing ratio units'''],

IND_RECEPTOR        IND_RECEPTOR switches between different units for
                    concentrations at the receptor 1=mass units (concentrations)
                    2=mass mixing ratio units

MQUASILAG           MQUASILAG indicates whether particles shall be numbered
                    consecutively (1) or with their release location number (0).
                    The first option allows tracking of individual particles
                    using the partposit output files

NESTED_OUTPUT       NESTED_OUTPUT decides whether model output shall be made
                    also for a nested output field (normally with higher resolution)
LINIT_COND          For Backward Runs, sets initial conditions:
                    [0]=No, 1=Mass Unit, 2=Mass Mixing

================    =================================================================

Definition at line 90 of file pflexible.py.

90 
91 def read_command(path, headerrows=7):
92  """
93  Reads a FLEXPART COMMAND file.
94 
95  .. note::
96  Assumes releases file has a header of 7 lines. Use
97  option "headerrows" to override this.
98 
99  USAGE::
100 
101  > A = read_command(filepath)
102 
103  where filepath is either a file object or a path.
104 
105 
106  Returns
107  a dictionary with the following keys
108 
109  ================ =================================================================
110  fields description
111  ================ =================================================================
112  SIM_DIR Simulation direction
113  SIM_START Text str of YYYYMMDD HHMMSS
114  SIM_END Text str of YYYYMMDD HHMMSS
115  AVG_CNC_INT Average concentrations are calculated every SSSSS seconds
116  AVG_CNC_TAVG The average concentrations are time averages of SSSSS sec
117  CNC_SAMP_TIME The concentrations are sampled every SSSSS seconds to
118  calculate the time average concentration.
119  T_PARTSPLIT Time constant for particle splitting.
120  SYNC All processes are synchronized with this time interval
121  CTL --
122  IFINE IFINE=Reduction factor for time step used for vertical wind
123  IOUT IOUT determines how the output shall be made: concentration
124  (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume
125  trajectory mode, or concentration + plume trajectory mode.
126  IPOUT IPOUT determines whether particle positions are outputted
127  (in addition to the gridded concentrations or mixing ratios)
128  or not. 0=no output, 1 output every output interval, 2 only
129  at end of the simulation
130  LSUBGRID Switch on/off subgridscale terrain parameterization
131  (increase of mixing heights due to subgridscale orog. var
132  LCONVECTION Switch on/off the convection parameterization
133  LAGESPECTRA Switch on/off the calculation of age spectra: if yes, the
134  file AGECLASSES must be available
135  IPIN If IPIN=1, a file "partposit_end" from a previous run must
136  be available in the output directory. Particle positions
137  are read in and previous simulation is continued. If
138  IPIN=0, no particles from a previous run are used
139  IOFR Switch on/off writing out each release.
140  IFLUX If IFLUX is set to 1, fluxes of each species through each
141  of the output boxes are calculated. Six fluxes,
142  corresponding to northward, southward, eastward, westward,
143  upward and downward are calculated for each grid cell of
144  the output grid. The control surfaces are placed in the
145  middle of each output grid cell. If IFLUX is set to 0,
146  no fluxes are determined.
147  MDOMAINFILL If MDOMAINFILL is set to 1, the first box specified in file
148  RELEASES is used as the domain where domain-filling
149  trajectory calculations are to be done. Particles are
150  initialized uniformly distributed (according to the air mass
151  distribution) in that domain at the beginning of the
152  simulation, and are created at the boundaries throughout
153  the simulation period
154  IND_SOURCE IND_SOURCE switches between different units for
155  concentrations at the source. NOTE that in backward
156  simulations the release of computational particles takes
157  place at the "receptor" and the sampling of particles at
158  the "source". 1=mass units (for bwd-runs = concentration)
159  2=mass mixing ratio units'''],
160 
161  IND_RECEPTOR IND_RECEPTOR switches between different units for
162  concentrations at the receptor 1=mass units (concentrations)
163  2=mass mixing ratio units
164 
165  MQUASILAG MQUASILAG indicates whether particles shall be numbered
166  consecutively (1) or with their release location number (0).
167  The first option allows tracking of individual particles
168  using the partposit output files
169 
170  NESTED_OUTPUT NESTED_OUTPUT decides whether model output shall be made
171  also for a nested output field (normally with higher resolution)
172  LINIT_COND For Backward Runs, sets initial conditions:
173  [0]=No, 1=Mass Unit, 2=Mass Mixing
174 
175  ================ =================================================================
176 
177  """
178  lines = file(path, 'r').readlines()
179  command_vals = [i.strip() for i in lines[headerrows:]] # clean line ends
180  COMMAND_KEYS = (
181  'SIM_DIR',
182  'SIM_START',
183  'SIM_END',
184  'AVG_CNC_INT',
185  'AVG_CNC_TAVG',
186  'CNC_SAMP_TIME',
187  'T_PARTSPLIT',
188  'SYNC',
189  'CTL',
190  'IFINE',
191  'IOUT',
192  'IPOUT',
193  'LSUBGRID',
194  'LCONVECTION',
195  'LAGESPECTRA',
196  'IPIN',
197  'OUTPUTFOREACHRELEASE',
198  'IFLUX',
199  'MDOMAINFILL',
200  'IND_SOURCE',
201  'IND_RECEPTOR',
202  'MQUASILAG',
203  'NESTED_OUTPUT',
204  'LINIT_COND')
205  float_keys = ['CTL']
206  date_keys = ['SIM_START', 'SIM_END']
207 
208  COMMAND = {}
209 
210  for i, key in enumerate(COMMAND_KEYS):
211  val = command_vals[i].split()
212  if key in date_keys:
213  val = val[:2]
214  elif key in float_keys:
215  val = float(val[0])
216  else:
217  val = int(val[0])
218 
219  COMMAND[key] = val
220  return COMMAND
221 
222 
223 
flextest.flexread.pflexible.read_command
def read_command
Functions for reading FLEXPART output #####.
Definition: pflexible.py:90
def flextest.flexread.pflexible.read_emissions (   emissionsfile,
  E = None,
  maxemissions = 1 
) 
Use reademissions.so module to read emissions file

Definition at line 2519 of file pflexible.py.

2519 
2520 def read_emissions(emissionsfile, E=None, maxemissions=1):
2521  """ Use reademissions.so module to read emissions file """
2522  try:
2523  from FortFlex import reademissions
2524  except:
2525  raise ImportError("Cannot find FortFlex module or missing reademissions")
2526  if not E:
2527  E = Structure()
2528  # set defaults for global 0.5 degree emissions
2529  defaults = {'nxmax':720, 'nymax':360, 'outlon0':-180, 'outlat0':-90, \
2530  'numxgrid':720, 'numygrid':360, 'dxout':0.5, 'dyout':0.5}
2531  for k, v in defaults.iteritems():
2532  exec("E.%s = %s" % (k, v))
2533 
2534  emissions = reademissions(emissionsfile, maxemissions, E.nxmax, E.nymax, \
2535  E.outlon0, E.outlat0, E.numxgrid, E.numygrid, E.dxout, E.dyout)
2536  E.grid = emissions
2537  E.filename = emissionsfile
2538  return E
2539 
def flextest.flexread.pflexible.read_grid (   H,
  kwargs 
) 
Accepts a header object as input, returns dictionary of Grid values
keyed by species and datestring from H['available_dates']. A grid instance from this
dictionary may be passed to the get_slabs function to return a dictionary of slabs 
suitable for plotting with the plotting routines. See below for more information.

**DEPENDENCY**
    Requires FortFlex.so module compiled using f2py. See FortFlex.f for more details.

Usage::

    > FLEXDATA = read_grid(H,**kwargs)

Returns:

    A grid dictionary key by tuples with (species,available_dates), where species is
    and integer. See grid.keys()

    FLEXDATA[(s,datestring)]['grid']
    FLEXDATA[(s,datestring)]['itime']
    FLEXDATA[(s,datestring)]['shape']
    FLEXDATA[(s,datestring)]['max']
    FLEXDATA[(s,datestring)]['min']
    FLEXDATA[(s,datestring)]['timestamp']
    FLEXDATA[(s,datestring)]['species']
    FLEXDATA[(s,datestring)]['gridfile']
    FLEXDATA[(s,datestring)]['rel_i']
    FLEXDATA[(s,datestring)]['spec_i']

Arguments

  .. tabularcolumns::  |l|L|

  =============         ========================================
  keyword               Description [default]
  =============         ========================================
  date                  which yyyymmddhhmmss from available_dates
                        or use (time_ret)
  time_ret              index to time
  unit                  'conc', 'pptv', ['time'], 'footprint'
  nspec_ret             numspecies
  pspec_ret             index to ???
  age_ret               index to ageclass
  nested                obtained from H['nested']
  BinaryFile            Use BinaryFile vs. FortFlex [False]
  getwet                True, [False]
  getdry                True, [False]
  scaledepo             A float value to scale deposition [1.0]
  scaleconc             A float value to scale conc [1.0]
  decayconstant         A float for decay const. [9e6]
  calcfoot              Will cause footprint to be calculated
                        [False]
  verbose               more output
  version               A string 'V8' or 'V6', ['V8']
  =============         ========================================


.. note::
    most arguments are able to be extracted fro the header "H"

Definition at line 1804 of file pflexible.py.

1804 
1805 def read_grid(H, **kwargs):
1806  """
1807  Accepts a header object as input, returns dictionary of Grid values
1808  keyed by species and datestring from H['available_dates']. A grid instance from this
1809  dictionary may be passed to the get_slabs function to return a dictionary of slabs
1810  suitable for plotting with the plotting routines. See below for more information.
1811 
1812  **DEPENDENCY**
1813  Requires FortFlex.so module compiled using f2py. See FortFlex.f for more details.
1814 
1815  Usage::
1816 
1817  > FLEXDATA = read_grid(H,**kwargs)
1818 
1819  Returns:
1820 
1821  A grid dictionary key by tuples with (species,available_dates), where species is
1822  and integer. See grid.keys()
1823 
1824  FLEXDATA[(s,datestring)]['grid']
1825  FLEXDATA[(s,datestring)]['itime']
1826  FLEXDATA[(s,datestring)]['shape']
1827  FLEXDATA[(s,datestring)]['max']
1828  FLEXDATA[(s,datestring)]['min']
1829  FLEXDATA[(s,datestring)]['timestamp']
1830  FLEXDATA[(s,datestring)]['species']
1831  FLEXDATA[(s,datestring)]['gridfile']
1832  FLEXDATA[(s,datestring)]['rel_i']
1833  FLEXDATA[(s,datestring)]['spec_i']
1834 
1835  Arguments
1836 
1837  .. tabularcolumns:: |l|L|
1838 
1839  ============= ========================================
1840  keyword Description [default]
1841  ============= ========================================
1842  date which yyyymmddhhmmss from available_dates
1843  or use (time_ret)
1844  time_ret index to time
1845  unit 'conc', 'pptv', ['time'], 'footprint'
1846  nspec_ret numspecies
1847  pspec_ret index to ???
1848  age_ret index to ageclass
1849  nested obtained from H['nested']
1850  BinaryFile Use BinaryFile vs. FortFlex [False]
1851  getwet True, [False]
1852  getdry True, [False]
1853  scaledepo A float value to scale deposition [1.0]
1854  scaleconc A float value to scale conc [1.0]
1855  decayconstant A float for decay const. [9e6]
1856  calcfoot Will cause footprint to be calculated
1857  [False]
1858  verbose more output
1859  version A string 'V8' or 'V6', ['V8']
1860  ============= ========================================
1861 
1862 
1863  .. note::
1864  most arguments are able to be extracted fro the header "H"
1865 
1866  """
1867 
1868  if H.version == 'V9':
1869  return readgridV9(H, **kwargs)
1870  if H.version == 'V8':
1871  return readgridV8(H, **kwargs)
1872  if H.version == 'V6':
1873  return readgridV6(H, **kwargs)
1874  else:
1875  raise IOError("No version attribute defined for Header.")
1876 
1877 

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.read_header (   pathname,
  kwargs 
) 
The readheader function returns a special class (Structure) which behaves
like a dictionary. It contains all the metadata from the simulation which
is contained in the "header" or "header_nest" binary files from the model 
output.

.. warning::
    It is recommended to use the :class:`Header` class: H = pf.Header(path)

This version is using the BinaryFile class rather than FortFlex. 

Usage::

    > H = read_header(pathname) #Don't include header filename

Returns a dictionary

    H = dictionary like object with all the run metadata.
Arguments

  .. tabularcolumns::  |l|L|

  =============       ========================================
  keyword             Description [default]
  =============       ========================================
  pathname            FLEXPART run output directory
  readp               read release points 0=no, [1]=y
  nested              nested output [False] or True
  headerfile          provide custom name for header file
  datefile            provide a custom name for the date file
  verbose             print information while loading header
  =============       ========================================

.. note::
    **This function is in development**

    Please report any bugs found.

.. TODO::

    probably a lot of things, among which ... 
    
   [] choose skip/getbin or direct seek/read 
   [] define output keys in docstring
   
   
.. note::
    The user is no longer required to indicate which version of FLEXPART
    the header is from. Checks are made, and the header is read accordingly.
    Please report any problems...

Definition at line 813 of file pflexible.py.

813 
814 def read_header(pathname, **kwargs):
815  """
816  The readheader function returns a special class (Structure) which behaves
817  like a dictionary. It contains all the metadata from the simulation which
818  is contained in the "header" or "header_nest" binary files from the model
819  output.
820 
821  .. warning::
822  It is recommended to use the :class:`Header` class: H = pf.Header(path)
823 
824  This version is using the BinaryFile class rather than FortFlex.
825 
826  Usage::
827 
828  > H = read_header(pathname) #Don't include header filename
829 
830  Returns a dictionary
831 
832  H = dictionary like object with all the run metadata.
833  Arguments
834 
835  .. tabularcolumns:: |l|L|
836 
837  ============= ========================================
838  keyword Description [default]
839  ============= ========================================
840  pathname FLEXPART run output directory
841  readp read release points 0=no, [1]=y
842  nested nested output [False] or True
843  headerfile provide custom name for header file
844  datefile provide a custom name for the date file
845  verbose print information while loading header
846  ============= ========================================
847 
848  .. note::
849  **This function is in development**
850 
851  Please report any bugs found.
852 
853  .. TODO::
854 
855  probably a lot of things, among which ...
856 
857  [] choose skip/getbin or direct seek/read
858  [] define output keys in docstring
859 
860 
861  .. note::
862  The user is no longer required to indicate which version of FLEXPART
863  the header is from. Checks are made, and the header is read accordingly.
864  Please report any problems...
865 
866  """
867 
868  OPS = Structure()
869  OPS.readp = True
870  OPS.nested = False
871  OPS.ltopo = 1 # 1 for AGL, 0 for ASL
872  OPS.verbose = False
873  OPS.headerfile = None
874  OPS.datefile = None
875 
876  # # add keyword overides and options to header
877  for k in kwargs.keys():
878  if k not in OPS.keys():
879  print("WARNING: {0} not a valid input option.".format(k))
880 
881  # BW compat fixes
882  if 'nest' in kwargs.keys():
883  raise IOError("nest is no longer a valid keyword, see docs. \n Now use nested=True or nested=False")
884 
885 
886  if 'nested' in kwargs.keys():
887  if kwargs['nested'] is 1:
888  print("WARNING, use of nested=1, deprecated converting to nested=True")
889  kwargs['nested'] = True
890 
891  if 'nested' in kwargs.keys():
892  if kwargs['nested'] is 0:
893  print("WARNING, use of nested=0, deprecated converting to nested=False")
894  kwargs['nested'] = False
895 
896  OPS.update(kwargs)
897 
898  if OPS.verbose:
899  print "Reading Header with:\n"
900 
901  for o in OPS:
902  print "%s ==> %s" % (o, OPS[o])
903 
904  # Define utility functions for reading binary file
905  skip = lambda n = 8 : bf.seek(n, 1)
906  getbin = lambda dtype, n = 1 : bf.read(dtype, (n,))
907 
908 
909  h = Structure()
910 
911 
912  if OPS.headerfile:
913  filename = os.path.join(pathname, OPS.headerfile);
914 
915  elif OPS.nested is True:
916  filename = os.path.join(pathname, 'header_nest')
917  h['nested'] = True;
918  else:
919  filename = os.path.join(pathname, 'header')
920  h['nested'] = False;
921 
922  # Open header file in binary format
923  if not os.path.exists(filename):
924  raise IOError("No such file: {0}".format(filename))
925  else:
926  try:
927  bf = BinaryFile(filename, order="fortran")
928  except:
929  raise IOError("Error opening: {0} with BinaryFile class".format(filename))
930 
931  # Get available_dates from dates file in same directory as header
932  if OPS.datefile:
933  datefile = os.path.join(pathname, OPS.datefile)
934  else:
935  datefile = os.path.join(pathname, 'dates')
936 
937  if not os.path.exists(datefile):
938  raise IOError("No DATEFILE: {0}".format(datefile))
939  else:
940  try:
941  fd = file(datefile, 'r').readlines()
942  except:
943  raise IOError("Could not read datefile: {0}".format(datefile))
944 
945  # get rid of any duplicate dates (a fix for the forecast system)
946  fd = sorted(list(set(fd)))
947  h['available_dates'] = [d.strip('\n') for d in fd]
948 
949 
950  # which version format is header file:
951  version = _get_header_version(bf)
952 
953 
954  # required containers
955  junk = [] # for catching unused output
956  h['nz_list'] = []
957  h['species'] = []
958  h['wetdep'] = []
959  h['drydep'] = []
960  h['ireleasestart'] = []
961  h['ireleaseend'] = []
962  h['compoint'] = []
963 
964  # Define Header format and create Dictionary Keys
965  I = {0:'_0', 1:'ibdate', 2:'ibtime', 3:'flexpart', \
966  4:'_1', 5:'loutstep', 6:'loutaver', 7:'loutsample', \
967  8:'_2', 9:'outlon0', 10:'outlat0', 11:'numxgrid', \
968  12:'numygrid', 13:'dxout', 14:'dyout', 15:'_3', 16:'numzgrid', \
969  }
970  # format for binary reading first part of the header file
971  Dfmt = ['i', 'i', 'i', '13S', '2i', 'i', 'i', 'i', '2i', 'f', 'f', 'i', 'i', 'f', 'f', '2i', 'i']
972  if bf:
973  a = [bf.read(fmt) for fmt in Dfmt]
974  for j in range(len(a)):
975  h[I[j]] = a[j]
976 
977  h['outheight'] = np.array([bf.read('f') for i in range(h['numzgrid'])])
978  junk.append(bf.read('2i'))
979  h['jjjjmmdd'] = bf.read('i')
980  h['hhmmss'] = bf.read('i')
981  junk.append(bf.read('2i'))
982  h['nspec'] = bf.read('i') / 3 # why!?
983  h['numpointspec'] = bf.read('i')
984  junk.append(bf.read('2i'))
985 
986  # Read in the species names and levels for each nspec
987  # temp dictionaries
988  for i in range(h['nspec']):
989  if 'V6' in version:
990 
991  h['wetdep'].append(''.join([bf.read('c') for i in range(10)]).strip())
992  junk.append(bf.read('2i'))
993  junk.append(bf.read('i'))
994  h['drydep'].append(''.join([bf.read('c') for i in range(10)]).strip())
995  junk.append(bf.read('2i'))
996  h['nz_list'].append(getbin('i'))
997  h['species'].append(''.join([getbin('c') for i in range(10)]).strip())
998 
999  else:
1000 
1001  junk.append(bf.read('i'))
1002  h['wetdep'].append(''.join([bf.read('c') for i in range(10)]).strip())
1003  junk.append(bf.read('2i'))
1004  junk.append(bf.read('i'))
1005  h['drydep'].append(''.join([bf.read('c') for i in range(10)]).strip())
1006  junk.append(bf.read('2i'))
1007  h['nz_list'].append(bf.read('i'))
1008  h['species'].append(''.join([bf.read('c') for i in range(10)]).strip())
1009  junk.append(bf.read('2i'))
1010 
1011  if 'V6' in version:
1012  bf.seek(8, 1)
1013  # pdb.set_trace()
1014  h['numpoint'] = bf.read('i')
1015 
1016  # read release info if requested
1017  # if OPS.readp pass has changed, we cycle through once,
1018  # then break the loop if OPS.readp is false. This is done
1019  # in order to get some of the information into the header.
1020  before_readp = bf.tell()
1021 
1022  # initialise release fields
1023  I = {2:'kindz', 3:'xp1', 4:'yp1', 5:'xp2', \
1024  6:'yp2', 7:'zpoint1', 8:'zpoint2', 9:'npart', 10:'mpart'}
1025 
1026  for k, v in I.iteritems():
1027  h[v] = np.zeros(h['numpoint']) # create zero-filled lists in H dict
1028 
1029  h['xmass'] = np.zeros((h['numpoint'], h['nspec']))
1030 
1031 
1032  if 'V6' in version:
1033  skip()
1034  _readV6(bf, h)
1035  else:
1036  junk.append(bf.read('i'))
1037  for i in range(h['numpoint']):
1038  junk.append(bf.read('i'))
1039 
1040  h['ireleasestart'].append(bf.read('i'))
1041  h['ireleaseend'].append(bf.read('i'))
1042  h['kindz'][i] = bf.read('h') # This is an int16, might need to to change something
1043  junk.append(bf.read('2i'))
1044 
1045  h['xp1'][i] = bf.read('f')
1046  h['yp1'][i] = bf.read('f')
1047  h['xp2'][i] = bf.read('f')
1048  h['yp2'][i] = bf.read('f')
1049  h['zpoint1'][i] = bf.read('f')
1050  h['zpoint2'][i] = bf.read('f')
1051 
1052  junk.append(bf.read('2i'))
1053  h['npart'][i] = bf.read('i')
1054  h['mpart'][i] = bf.read('i')
1055 
1056  junk.append(bf.read('i'))
1057  # initialise release fields
1058 
1059  l = bf.read('i') # get compoint length?
1060  gt = bf.tell() + l # create 'goto' point
1061  sp = ''
1062  while re.search("\w", bf.read('c')): # collect the characters for the compoint
1063  bf.seek(-1, 1)
1064  sp = sp + bf.read('c')
1065 
1066  bf.seek(gt) # skip ahead to gt point
1067 
1068  h['compoint'].append(sp) # species names in dictionary for each nspec
1069  # h['compoint'].append(''.join([bf.read('c') for i in range(45)]))
1070 
1071  junk.append(bf.read('i'))
1072  # now loop for nspec to get xmass
1073  for v in range(h['nspec']):
1074  Dfmt = ['i', 'f', '2i', 'f', '2i', 'f', 'i']
1075  a = [bf.read(fmt) for fmt in Dfmt]
1076  h['xmass'][i, v] = a[1]
1077 
1078  if OPS.readp is False:
1079  """
1080  We get the first set of release points here in order
1081  to get some information, but skip reading the rest
1082  """
1083  bf.seek(before_readp)
1084  if 'V6' in version:
1085  bf.seek(157 * h['numpoint'], 1);
1086  else:
1087  bf.seek(119 * h['numpoint'] + (h['nspec'] * 36) * h['numpoint'] + 4, 1)
1088  break
1089 
1090 
1091 
1092  if 'V6' in version:
1093  h['method'] = bf.read('i')
1094 
1095  else:
1096  junk.append(bf.read('i'))
1097  junk.append(bf.read('i'))
1098 
1099  h['lsubgrid'] = bf.read('i')
1100  h['lconvection'] = bf.read('i')
1101  h['ind_source'] = bf.read('i')
1102  h['ind_receptor'] = bf.read('i')
1103 
1104  if 'V6' in version:
1105  pass
1106  else:
1107  junk.append(bf.read('2i'))
1108 
1109  h['nageclass'] = bf.read('i')
1110 
1111  Lage_fmt = ['i'] * h.nageclass
1112  h['lage'] = [bf.read(fmt) for fmt in Lage_fmt]
1113 
1114  # Orography
1115  nx = h['numxgrid']
1116  ny = h['numygrid']
1117  Dfmt = ['f'] * nx
1118  h['oro'] = np.zeros((nx, ny), np.float)
1119  junk.append(bf.read('2i'))
1120 
1121  for ix in range(nx):
1122  h['oro'][ix] = bf.read('f', ny)
1123  bf.seek(8, 1)
1124 
1125  # Why was this? / deprecated.
1126  # if h['loutstep'] < 0:
1127  # h['nspec'] = h['numpoint']
1128 
1129  bf.close()
1130 
1131 
1132 
1133 
1134  ############# ADDITIONS ###########
1135  # attributes to header that can be added after reading below here
1136  h['pathname'] = pathname
1137  h['decayconstant'] = 0
1138  h.path = pathname
1139 
1140 
1141  # Calculate Height (outheight + topography)
1142  # There is an offset issue here related to the 0-indexing. Be careful.
1143  oro = h['oro'] # z is a numpy array
1144  nx = h['numxgrid']
1145  ny = h['numygrid']
1146  nz = h['numzgrid']
1147 
1148  Heightnn = np.zeros((nx, ny, nz), np.float)
1149  for ix in range(nx):
1150  if OPS.ltopo == 1:
1151  Heightnn[ix, :, 0] = oro[ix, :]
1152  else:
1153  Heightnn[ix, :, 0] = np.zeros(ny)
1154 
1155  for iz in range(nz):
1156  if OPS.ltopo == 1:
1157  Heightnn[ix, :, iz] = [h['outheight'][iz] + oro[ix, y] for y in range(ny)]
1158  else:
1159  Heightnn[ix, :, iz] = h['outheight'][iz]
1160 
1161  h['Area'] = gridarea(h)
1162  h['Heightnn'] = Heightnn
1163  h['nx'] = nx
1164  h['ny'] = ny
1165  h['nz'] = nz
1166 
1167  # Convert ireleasestart and ireleaseend to datetimes, fix issue #10
1168  start_day = datetime.datetime.strptime(str(h['ibdate']), '%Y%m%d')
1169  H, M, S , = [int(str(h['ibtime']).zfill(6)[i:i + 2]) for i in range(0, 6, 2)]
1170  start_time = datetime.timedelta(hours=H, minutes=M, seconds=S)
1171  h['simulationstart'] = start_day + start_time
1172 
1173  if OPS.readp:
1174  releasestart, releaseend = [], []
1175  for i in range(h.numpointspec):
1176  releasestart.append(h.simulationstart + \
1177  datetime.timedelta(seconds=int(h.ireleasestart[i])))
1178  releaseend.append(h.simulationstart + \
1179  datetime.timedelta(seconds=int(h.ireleaseend[i])))
1180  h.releasestart = releasestart
1181  h.releaseend = releaseend[:h.numpointspec]
1182  h.releasetimes = [b - ((b - a) / 2) for a, b in zip(h.releasestart, h.releaseend)]
1183 
1184  # Add datetime objects for dates
1185  available_dates_dt = []
1186  for i in h.available_dates:
1187  available_dates_dt.append(datetime.datetime(
1188  int(i[:4]), int(i[4:6]), int(i[6:8]), int(i[8:10]), int(i[10:12]), int(i[12:])))
1189  h.available_dates_dt = available_dates_dt
1190  h.first_date = available_dates_dt[0]
1191  h.last_date = available_dates_dt[-1]
1192  h.ageclasses = np.array([act - h.simulationstart for act in h.available_dates_dt])
1193  h.numageclasses = len(h.ageclasses)
1194 
1195  # Add other helpful attributes
1196  h.nxmax = h.numxgrid
1197  h.nymax = h.numygrid
1198  h.nzmax = h.numzgrid
1199  h.maxspec = h.nspec
1200  h.maxpoint = h.numpoint
1201  h.maxageclass = h.numageclasses
1202 
1203  h.area = h.Area # fix an annoyance
1204 
1205  if OPS.readp:
1206  h.xpoint = h.xp1
1207  h.ypoint = h.yp1
1208 
1209  # Add release unit derived from kindz
1210  if 'kindz' not in h.keys():
1211  h.kindz = [0]
1212  h.alt_unit = 'unkn.'
1213  if 3 in h.kindz:
1214  h.alt_unit = 'hPa'
1215  elif 2 in h.kindz:
1216  h.alt_unit = 'm.a.s.l.'
1217  elif 1 in h.kindz:
1218  h.alt_unit = 'm.a.g.l.'
1219 
1220  #
1221  if h.loutstep > 0:
1222  h.direction = 'forward'
1223  h.unit = 'conc' # could be pptv
1224  h.plot_unit = 'ppb'
1225  else:
1226  h.direction = 'backward'
1227  h.unit = 'time'
1228  h.plot_unit = 'ns / kg' # Not sure about this
1229 
1230  # Units based on Table 1, ACP 2005
1231  if h.direction == 'forward':
1232  if h.ind_source == 1:
1233  if h.ind_receptor == 1:
1234  h.output_unit = 'ng m-3'
1235  if h.ind_receptor == 2:
1236  h.output_unit = 'pptm'
1237  if h.ind_source == 2:
1238  if h.ind_receptor == 1:
1239  h.output_unit = 'ng m-3'
1240  if h.ind_receptor == 2:
1241  h.output_unit = 'pptm'
1242  if h.direction == 'backward':
1243  if h.ind_source == 1:
1244  if h.ind_receptor == 1:
1245  h.output_unit = 's'
1246  if h.ind_receptor == 2:
1247  h.output_unit = 's m^3 kg-1'
1248  if h.ind_source == 2:
1249  if h.ind_receptor == 1:
1250  h.output_unit = 's kg m-3'
1251  if h.ind_receptor == 2:
1252  h.output_unit = 's'
1253 
1254 
1255 
1256  # Add layer thickness
1257  layerthickness = [h.outheight[0]]
1258  for i, lh in enumerate(h.outheight[1:]):
1259  layerthickness.append(lh - h.outheight[i])
1260  h.layerthickness = layerthickness
1261 
1262  h.options = OPS
1263  h.fp_version = h.flexpart
1264  h.junk = junk # the extra bits that were read... more for debugging
1265 
1266 
1267  print('Read {0} Header: {1}'.format(h.fp_version, filename))
1268 
1269  return h
1270 
# BW Compatability
flextest.flexread.pflexible.BinaryFile
HELPER FUNCTIONS ##########.
Definition: pflexible.py:4484

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.read_releases (   path,
  headerrows = 11 
) 
Reads a FLEXPART releases file.

.. note::
    Assumes releases file has a header of 11 lines. Use
    option "headerrows" to override this.

USAGE::

    > A = read_releases(filepath)

where filepath is either a file object or a path.


Returns
  a record array with fields:

  ============      ==========================
  fields            description
  ============      ==========================
  start_time        datetime start
  end_time          datetime end
  lllon             lower left longitude
  llat              lower left latitude
  urlon             upper right longitude
  urlat             upper right latitude
  altunit           1=magl, 2=masl, 3=hPa
  elv1              lower z level
  elv2              upper z level
  numpart           numparticles
  mass              mass for each spec, so the
                    array will actually have
                    fields: mass0, mass1,..
  id                ID for each release
  ============      ==========================

Definition at line 224 of file pflexible.py.

224 
225 def read_releases(path, headerrows=11):
226  """
227  Reads a FLEXPART releases file.
228 
229  .. note::
230  Assumes releases file has a header of 11 lines. Use
231  option "headerrows" to override this.
232 
233  USAGE::
234 
235  > A = read_releases(filepath)
236 
237  where filepath is either a file object or a path.
238 
239 
240  Returns
241  a record array with fields:
242 
243  ============ ==========================
244  fields description
245  ============ ==========================
246  start_time datetime start
247  end_time datetime end
248  lllon lower left longitude
249  llat lower left latitude
250  urlon upper right longitude
251  urlat upper right latitude
252  altunit 1=magl, 2=masl, 3=hPa
253  elv1 lower z level
254  elv2 upper z level
255  numpart numparticles
256  mass mass for each spec, so the
257  array will actually have
258  fields: mass0, mass1,..
259  id ID for each release
260  ============ ==========================
261 
262  """
263  lines = _getfile_lines(path)
264  lines = [i.strip() for i in lines] # clean line ends
265 
266  # we know nspec is at line 11
267  nspec = int(lines[headerrows])
268  blocklength = headerrows + nspec
269  spec = []
270  for i in range(nspec):
271  spec.append(int(lines[headerrows + 1 + i]))
272  indx = headerrows + 1 + (i + 2)
273  blocks = []
274  for i in range(indx, len(lines), blocklength + 1):
275  blocks.append(lines[i:i + blocklength])
276  for b in blocks:
277  b[0] = datetime.datetime.strptime(b[0], '%Y%m%d %H%M%S')
278  b[1] = datetime.datetime.strptime(b[1], '%Y%m%d %H%M%S')
279  b[2:6] = [np.float(i) for i in b[2:6]]
280  b[6] = int(b[6])
281  b[7] = float(b[7])
282  b[8] = float(b[8])
283  b[9] = int(b[9])
284  for i in range(nspec):
285  b[10 + i] = float(b[10 + i])
286  b = tuple(b)
287 
288 
289  names = ['start_time', 'end_time', 'lllon', 'lllat', 'urlon', 'urlat', \
290  'altunit', 'elv1', 'elv2', 'numpart']
291  # formats = [object, object, np.float, np.float, np.float, np.float,\
292  # int, np.float, np.float, int]
293  for i in range(nspec):
294  names.append('mass%s' % i)
295  # formats.append(np.float)
296  names.append('id')
297  # formats.append('S30')
298 
299  # dtype = {'names':names, 'formats':formats}
300  # RELEASES = np.rec.array(blocks,dtype=dtype)
301  return np.rec.fromrecords(blocks, names=names)
302 

Here is the call graph for this function:

def flextest.flexread.pflexible.read_trajectories (   H,
  trajfile = 'trajectories.txt',
  ncluster = 5,
  ageclasses = 20 
) 
Reads the trajectories.txt file in a FLEXPART run output directory.

Based on output from plumetraj.f::

    ********************************************************************************
    *                                                                              *
    * Determines a plume centroid trajectory for each release site, and manages    *
    * clustering of particle locations. Certain parameters (average PV,            *
    * tropopause height, etc., are provided along the plume trajectories.          *
    * At the end, output is written to file 'trajectories.txt'.                    *
    *                                                                              *
    *     Author: A. Stohl                                                         *
    *                                                                              *
    *     24 January 2002                                                          *
    *                                                                              *
    * Variables:                                                                   *
    * fclust          fraction of particles belonging to each cluster              *
    * hmixcenter      mean mixing height for all particles                         *
    * ncluster        number of clusters to be used                                *
    * pvcenter        mean PV for all particles                                    *
    * pvfract         fraction of particles with PV<2pvu                           *
    * rms             total horizontal rms distance after clustering               *
    * rmsdist         total horizontal rms distance before clustering              *
    * rmsclust        horizontal rms distance for each individual cluster          *
    * topocenter      mean topography underlying all particles                     *
    * tropocenter     mean tropopause height at the positions of particles         *
    * tropofract      fraction of particles within the troposphere                 *
    * zrms            total vertical rms distance after clustering                 *
    * zrmsdist        total vertical rms distance before clustering                *
    * xclust,yclust,  Cluster centroid positions                                   *
    * zclust                                                                       *
    *                                                                              *
    ********************************************************************************

USAGE::

    > T = read_trajectories(H_OR_path_to_directory, **kwargs)

.. note::
    The first argument is either a path to a trajectory file, or simply a :class:`Header`
    instance.

Returns a dictionary of the trajectories for each release.

  .. tabularcolumns::  |l|L|

  =============       ========================================
  Keys                Description
  =============       ========================================
  Trajectories        array_of_floats(j,it1,xi,yi,zi,topoi,hmixi,tropoi,pvi,
                      rmsdisti,rmsi,zrmsdisti,zrmsi,hfri,pvfri,trfri,
                      (xclusti(k),yclusti(k),zclusti(k),fclusti(k),rmsclusti(k),k=1,5))
  RELEASE_ID          1array(i1,i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart)
  numpspec            number of species
  numageclass         number of ageclasses
  =============       ========================================


Arguments

  .. tabularcolumns::  |l|L|

  =============       ========================================
  keyword             Description [default]
  =============       ========================================
  trajfile            over the name of the input 
                      file ['trajectories.txt']
  ncluster            number of clusters [5]
  ageclasses          number of ageclasses [20]
  =============       ========================================

Definition at line 305 of file pflexible.py.

305 
306  ageclasses=20):
307  """
308  Reads the trajectories.txt file in a FLEXPART run output directory.
309 
310  Based on output from plumetraj.f::
311 
312  ********************************************************************************
313  * *
314  * Determines a plume centroid trajectory for each release site, and manages *
315  * clustering of particle locations. Certain parameters (average PV, *
316  * tropopause height, etc., are provided along the plume trajectories. *
317  * At the end, output is written to file 'trajectories.txt'. *
318  * *
319  * Author: A. Stohl *
320  * *
321  * 24 January 2002 *
322  * *
323  * Variables: *
324  * fclust fraction of particles belonging to each cluster *
325  * hmixcenter mean mixing height for all particles *
326  * ncluster number of clusters to be used *
327  * pvcenter mean PV for all particles *
328  * pvfract fraction of particles with PV<2pvu *
329  * rms total horizontal rms distance after clustering *
330  * rmsdist total horizontal rms distance before clustering *
331  * rmsclust horizontal rms distance for each individual cluster *
332  * topocenter mean topography underlying all particles *
333  * tropocenter mean tropopause height at the positions of particles *
334  * tropofract fraction of particles within the troposphere *
335  * zrms total vertical rms distance after clustering *
336  * zrmsdist total vertical rms distance before clustering *
337  * xclust,yclust, Cluster centroid positions *
338  * zclust *
339  * *
340  ********************************************************************************
341 
342  USAGE::
343 
344  > T = read_trajectories(H_OR_path_to_directory, **kwargs)
345 
346  .. note::
347  The first argument is either a path to a trajectory file, or simply a :class:`Header`
348  instance.
349 
350  Returns a dictionary of the trajectories for each release.
351 
352  .. tabularcolumns:: |l|L|
353 
354  ============= ========================================
355  Keys Description
356  ============= ========================================
357  Trajectories array_of_floats(j,it1,xi,yi,zi,topoi,hmixi,tropoi,pvi,
358  rmsdisti,rmsi,zrmsdisti,zrmsi,hfri,pvfri,trfri,
359  (xclusti(k),yclusti(k),zclusti(k),fclusti(k),rmsclusti(k),k=1,5))
360  RELEASE_ID 1array(i1,i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart)
361  numpspec number of species
362  numageclass number of ageclasses
363  ============= ========================================
364 
365 
366  Arguments
367 
368  .. tabularcolumns:: |l|L|
369 
370  ============= ========================================
371  keyword Description [default]
372  ============= ========================================
373  trajfile over the name of the input
374  file ['trajectories.txt']
375  ncluster number of clusters [5]
376  ageclasses number of ageclasses [20]
377  ============= ========================================
378 
379  """
380 
381  if isinstance(H, str):
382  try:
383  alltraj = file(H, 'r').readlines()
384  except:
385  raise IOError('Could not open file: %s' % H)
386  else:
387  path = H.path
388  alltraj = file(os.path.join(path, trajfile), 'r').readlines()
389 
390 
391  try:
392  ibdate, ibtime, model, version = alltraj[0].strip().split()[:4]
393  except:
394  ibdate, ibtime = alltraj[0].strip().split()[:2]
395  model = 'Flexpart'
396  version = 'V.x'
397  dt = datetime.datetime.strptime(ibdate + ibtime.zfill(6), '%Y%m%d%H%M%S')
398  numpoint = int(alltraj[2].strip())
399  # Fill a dictionary with the Release points and information keyed by name
400  # RelTraj['RelTraj_ID'] = (i1,i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart)
401  RelTraj = Structure()
402  Trajectories = []
403 
404  for i in range(3, 3 + (numpoint * 2), 2):
405  i1, i2, xp1, yp1, xp2, yp2, zp1, zp2, k, npart , = \
406  tuple([float(j) for j in alltraj[i].strip().split()])
407  itimerel1 = dt + datetime.timedelta(seconds=i1)
408  itimerel2 = dt + datetime.timedelta(seconds=i2)
409  Xp = (xp1 + xp2) / 2
410  Yp = (yp1 + yp2) / 2
411  Zp = (zp1 + zp2) / 2
412  RelTraj[alltraj[i + 1].strip()] = np.array((itimerel1, itimerel2, Xp, Yp, Zp, k, npart))
413 
414  for i in range(3 + (numpoint * 2), len(alltraj)):
415  raw = alltraj[i]
416  FMT = [0, 5, 8, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 6, 6] + ncluster * [8, 8, 7, 6, 8]
417 
418  data = [raw[sum(FMT[:ii]):sum(FMT[:ii + 1])] for ii in range(1, len(FMT) - 1)] + \
419  [raw[sum(FMT[:-1]):]]
420  ### FIX ###
421  # ## To get rid of '******' that is now in trajectories.txt
422  data = [float(r.replace('********', 'NaN')) for r in data]
423 
424  Trajectories.append(data)
425 
426  data = np.array(Trajectories)
427  RelTraj['version'] = model + ' ' + version
428  RelTraj['date'] = dt
429  RelTraj['Trajectories'] = data
430  RelTraj['labels'] = \
431  ['release number', 'seconds prior to release', 'lon', 'lat', 'height', 'mean topography', \
432  'mean mixing height', 'mean tropopause height', 'mean PV index', \
433  'rms distance', 'rms', 'zrms distance', 'zrms', \
434  'fraction height??', 'fraction PV<2pvu', 'fraction in troposphere'] + \
435  ncluster * ['xcluster', 'ycluster', 'zcluster', 'fcluster', 'rmscluster']
436  RelTraj['info'] = \
437  """
438  Returns a dictionary:
439  R['Trajectories'] = array_of_floats(
440  releasenum,it1,xi,yi,zi,topoi,hmixi,tropoi,pvi,
441  rmsdisti,rmsi,zrmsdisti,zrmsi,hfri,pvfri,trfri,
442  (xclusti(k),yclusti(k),zclusti(k),fclusti(k),rmsclusti(k),k=1,5))
443 
444  R['RELEASE_ID'] = (dt_i1,dt_i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart)
445  R['info'] = this message
446 
447  To plot a trajectory timeseries:
448  RT = read_trajectories(H)
449  T = RT['Trajectories']
450  rel = 1
451  t = T[np.where(T[:,0]==rel),:][0]
452  plt.plot(t[:,1],t[:,14])
453  plt.savefig('trajectories.png')
454  """
455  return RelTraj
456 

Here is the caller graph for this function:

def flextest.flexread.pflexible.readgridV6 (   H,
  kwargs 
) 
Accepts a header object as input, and selects appropriate readgrid function
to use for reading in data from the flexpart binary Fortran files.

See the :func:`read_grid` for more information on keyword arguments

This is the 'V6' version of the function.

Definition at line 2149 of file pflexible.py.

2149 
2150 def readgridV6(H, **kwargs):
2151 
2152  """
2153  Accepts a header object as input, and selects appropriate readgrid function
2154  to use for reading in data from the flexpart binary Fortran files.
2155 
2156  See the :func:`read_grid` for more information on keyword arguments
2157 
2158  This is the 'V6' version of the function.
2159 
2160  """
2161  # # OPS is the options Structure, sets defaults, then update w/ kwargs
2162  OPS = Structure()
2163  OPS.unit = 'time'
2164  OPS.nspec_ret = 0
2165  OPS.pspec_ret = 0
2166  OPS.age_ret = 0
2167  OPS.time_ret = 0
2168  OPS.scaledepo = 1.0
2169  OPS.scaleconc = 1.0
2170  OPS.decayconstant = 9e6
2171  OPS.date = None
2172  OPS.calcfoot = False
2173  OPS.verbose = False
2174  OPS.BinaryFile = False
2175  # # add keyword overides and options to header
2176  OPS.update(kwargs)
2177  # H.update(OPS)
2178 
2179  # # set up the return dictionary (FLEXDATA updates fd, fd is returned)
2180  FLEXDATA = {}
2181  fd = Structure()
2182  fd.options = Structure()
2183  # # add the requests to the fd object to be returned
2184  fd.options.update(OPS)
2185 
2186  # # What direction is the run?
2187  unit = OPS.unit
2188  if H['loutstep'] > 0:
2189  forward = True
2190  if unit == 'time':
2191  # # default forward unit
2192  unit = 'conc'
2193  else:
2194  forward = False
2195 
2196  # # What species to return?
2197  nspec_ret = OPS.nspec_ret
2198  if isinstance(nspec_ret, int):
2199  nspec_ret = [nspec_ret]
2200  assert iter(nspec_ret), "nspec_ret must be iterable."
2201 
2202  # # get times to return
2203  get_dates = None
2204  if OPS.time_ret is not None:
2205  get_dates = []
2206  time_ret = OPS.time_ret
2207  if isinstance(time_ret, int) == True:
2208  time_ret = [time_ret]
2209 
2210  if time_ret[0] < 0:
2211  if forward == False:
2212  # # get all dates for calculating footprint.
2213  time_ret = np.arange(len(H.available_dates))
2214  else:
2215  raise ValueError("Must enter a positive time_ret for forward runs")
2216 
2217  for t in time_ret:
2218  get_dates.append(H.available_dates[t])
2219 
2220 
2221  # # define what dates to extract if user has explicitly defined a 'date'
2222  if OPS.date != None:
2223  date = OPS.date
2224  if time_ret is not None:
2225  Warning("overwriting time_ret variable, date was requested")
2226  get_dates = []
2227  if not isinstance(date, list):
2228  date = date.strip().split(',')
2229  for d in date:
2230  try:
2231  get_dates.append(H.available_dates[H.available_dates.index(d)])
2232  time_ret = None
2233  except:
2234  _shout("Cannot find date: %s in H['available_dates']\n" % d)
2235 
2236  if get_dates is None:
2237  raise ValueError("Must provide either time_ret or date value.")
2238  else:
2239  # # assign grid dates for indexing fd
2240  fd.grid_dates = get_dates[:]
2241 
2242  #print 'getting grid for: ', get_dates
2243  # Some predifinitions
2244  fail = 0
2245  # set filename prefix
2246  prefix = ['grid_conc_', 'grid_pptv_', \
2247  'grid_time_', 'footprint_', 'footprint_total', \
2248  'grid_conc_nest_', 'grid_pptv_nest_', \
2249  'grid_time_nest_', 'footprint_nest_', 'footprint_total_nest'
2250  ]
2251 
2252  units = ['conc', 'pptv', 'time', 'footprint', 'footprint_total']
2253  unit_i = units.index(unit)
2254  # Determine what module to read, try to use FortFlex, then dumpgrid, lastly pure Python
2255  # import the FortFlex / Fortran module
2256  try:
2257  from FortFlex import readgrid_v6 as readgrid
2258  from FortFlex import sumgrid
2259  useFortFlex = True
2260  OPS.useFortFlex = useFortFlex
2261  print 'using FortFlex VERSION 6'
2262  except:
2263  # get the original module (no memory allocation)
2264  useFortFlex = False
2265  raise Warning('Cannot import FortFlex trying to use BinaryFile')
2266  if not useFortFlex:
2267  readgrid = _readgridBF
2268  OPS.BinaryFile = True
2269 
2270  # cheat: use key2var function to get values from header dict, H
2271  # need to change this to using H.xxxxx
2272  # headervars = ['nested','nspec','numxgrid','numygrid','numzgrid','nageclass',\
2273  # 'dates','pathname','decayconstant','numpoint','numpointspec',
2274  # 'area','Heightnn','lage']
2275  # for k in headervars:
2276  # key2var(H,k)
2277 
2278 
2279 
2280  # reserve output fields
2281  #print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint
2282 
2283  # -------------------------------------------------
2284 
2285  # if forward:
2286  # numpointspec = 1
2287  # else:
2288  # numpointspec = numpoint
2289  # numpointspec = H['numpointspec']
2290  # if unit_i == 4:
2291  # grid=np.empty((numxgrid,numygrid,1,nspec_ret,pspec_ret,age_ret,len(get_dates)))
2292  # else:
2293  # grid=np.empty((numxgrid,numygrid,numzgrid,nspec_ret,pspec_ret,age_ret,len(get_dates)))
2294  # zplot=np.empty((numxgrid,numygrid,numzgrid,numpointspec))
2295 
2296  #--------------------------------------------------
2297  # Loop over all times, given in field H['available_dates']
2298  #--------------------------------------------------
2299 
2300  for date_i in range(len(get_dates)):
2301  datestring = get_dates[date_i]
2302  #print datestring
2303  FLEXDATA[datestring] = {}
2304  for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1):
2305  total_footprint = False
2306  FLEXDATA[(s, datestring)] = Structure()
2307  # spec_fid = '_'+str(s+1).zfill(3)
2308 
2309  if unit_i != 4:
2310  filename = os.path.join(H['pathname'], \
2311  prefix[(unit_i) + (H.nested * 5)] + datestring)
2312  H.zdims = H.numzgrid
2313 
2314  else:
2315  # grid total footprint
2316  print "Total footprint"
2317  total_footprint = True
2318  filename = os.path.join(H['pathname'], \
2319  prefix[(unit_i) + (H.nested * 5)])
2320  H.zdims = 1
2321 
2322  if os.path.exists(filename):
2323  H.filename = filename
2324  # print 'reading: ' + filename
2325  if OPS.verbose:
2326  print 'with values:'
2327  inputvars = ['filename', 'numxgrid', 'numygrid',
2328  'zdims', 'numpoint', 'nageclass', \
2329  'scaledepo', 'scaleconc',
2330  'decayconstant', 'numpointspec']
2331  for v in inputvars:
2332  print v, " ==> ", H[v]
2333 
2334 
2335  if OPS.BinaryFile:
2336  # print 'Using BinaryFile'
2337  gridT, wetgrid, drygrid, itime = _readgridBF(H, filename)
2338  else:
2339  gridT, wetgrid, drygrid, itime = readgrid(filename, \
2340  H.numxgrid, H.numygrid,
2341  H.zdims, H.numpointspec, H.nageclass, \
2342  OPS.scaledepo, OPS.scaleconc, H.decayconstant)
2343 
2344 
2345  if forward:
2346  zplot = gridT[:, :, :, :, 0]
2347  else:
2348  zplot = gridT[:, :, :, :, 0]
2349 
2350  # if total_footprint:
2351  # zplot = np.squeeze(gridT)
2352 
2353  if OPS.calcfoot:
2354  zplot = sumgrid(zplot, gridT, \
2355  H.area, H.Heightnn)
2356 
2357 
2358  # # get the total column and prep the grid
2359  if H.direction == 'forward':
2360  # not trying to do anything here... must be done
2361  # after retrieving the grid
2362  # D = get_slabs(H,np.squeeze(zplot))
2363  # rel_i = H.available_dates.index(datestring)
2364  D = zplot
2365  rel_i = 'k'
2366  else:
2367  D = zplot
2368  rel_i = 'k'
2369 
2370  # # NOTE:
2371  # # If you're changing things here, you might want to change
2372  # # them in fill_backward as well, yes I know... something is
2373  # # poorly designed ;(
2374  FLEXDATA[(s, datestring)]['grid'] = D # zplot
2375  FLEXDATA[(s, datestring)]['itime'] = itime
2376  FLEXDATA[(s, datestring)]['shape'] = zplot.shape
2377  FLEXDATA[(s, datestring)]['max'] = zplot.max()
2378  FLEXDATA[(s, datestring)]['min'] = zplot.min()
2379  FLEXDATA[(s, datestring)]['timestamp'] = datetime.datetime.strptime(datestring, '%Y%m%d%H%M%S')
2380  FLEXDATA[(s, datestring)]['species'] = H['species'][s]
2381  FLEXDATA[(s, datestring)]['gridfile'] = filename
2382  FLEXDATA[(s, datestring)]['rel_i'] = rel_i
2383  FLEXDATA[(s, datestring)]['spec_i'] = s
2384 
2385 
2386  else:
2387  _shout('***ERROR: file %s not found! \n' % filename)
2388  fail = 1
2389  fd.set_with_dict(FLEXDATA)
2390 
2391  try:
2392  # just for testing, set the first available grid as a shortcut
2393  # this will be removed.
2394  qind = (0, fd.grid_dates[0])
2395  fd.grid = fd[qind][fd[qind].keys()[0]].grid
2396  except:
2397  pass
2398 
2399 
2400  return fd

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.readgridV8 (   H,
  kwargs 
) 
Accepts a header object as input, and selects appropriate readgrid function
to use for reading in data from the flexpart binary Fortran files.

See the :func:`read_grid` for more information on keyword arguments

This is the 'V8' version of the function.

Definition at line 1878 of file pflexible.py.

1878 
1879 def readgridV8(H, **kwargs):
1880 
1881  """
1882  Accepts a header object as input, and selects appropriate readgrid function
1883  to use for reading in data from the flexpart binary Fortran files.
1884 
1885  See the :func:`read_grid` for more information on keyword arguments
1886 
1887  This is the 'V8' version of the function.
1888 
1889  """
1890  # # OPS is the options Structure, sets defaults, then update w/ kwargs
1891  OPS = Structure()
1892  OPS.unit = H.unit
1893  OPS.getwet = False
1894  OPS.getdry = False
1895  OPS.nspec_ret = 0
1896  OPS.npspec_int = False # allows to select an index of npsec when calling readgrid
1897  OPS.pspec_ret = 0
1898  OPS.age_ret = 0
1899  OPS.time_ret = 0
1900  OPS.scaledepo = 1.0
1901  OPS.scaleconc = 1.0
1902  OPS.decayconstant = 9e6
1903  OPS.date = None
1904  OPS.calcfoot = False
1905  OPS.verbose = False
1906  OPS.BinaryFile = False
1907  OPS.version = 'V8'
1908  # # add keyword overides and options to header
1909  OPS.update(kwargs)
1910  # H.update(OPS)
1911 
1912  # # set up the return dictionary (FLEXDATA updates fd, fd is returned)
1913  FLEXDATA = {}
1914  fd = Structure()
1915  fd.options = Structure()
1916 
1917  # # What direction is the run?
1918  unit = OPS.unit
1919 
1920  if H['loutstep'] > 0:
1921  forward = True
1922  if unit == 'time':
1923  # # default forward unit
1924  unit = 'conc'
1925  OPS.unit = unit
1926  else:
1927  forward = False
1928 
1929  # # What species to return?
1930  nspec_ret = OPS.nspec_ret
1931  if isinstance(nspec_ret, int):
1932  nspec_ret = [nspec_ret]
1933  assert iter(nspec_ret), "nspec_ret must be iterable."
1934 
1935  # # get times to return
1936  get_dates = None
1937  if OPS.time_ret is not None:
1938  get_dates = []
1939  time_ret = OPS.time_ret
1940  if isinstance(time_ret, int) == True:
1941  time_ret = [time_ret]
1942 
1943  if time_ret[0] < 0:
1944  #if forward == False:
1945  # # get all dates for calculating footprint.
1946  time_ret = np.arange(len(H.available_dates))
1947  #else:
1948  # raise ValueError("Must enter a positive time_ret for forward runs")
1949 
1950  for t in time_ret:
1951  get_dates.append(H.available_dates[t])
1952 
1953 
1954  # # define what dates to extract if user has explicitly defined a 'date'
1955  if OPS.date != None:
1956  date = OPS.date
1957  if time_ret is not None:
1958  Warning("overwriting time_ret variable, date was requested")
1959  get_dates = []
1960  if not isinstance(date, list):
1961  date = date.strip().split(',')
1962  for d in date:
1963  try:
1964  get_dates.append(H.available_dates[H.available_dates.index(d)])
1965  time_ret = None
1966  except:
1967  _shout("Cannot find date: %s in H['available_dates']\n" % d)
1968 
1969  if get_dates is None:
1970  raise ValueError("Must provide either time_ret or date value.")
1971  else:
1972  # # assign grid dates for indexing fd
1973  fd.grid_dates = get_dates[:]
1974 
1975  #print 'getting grid for: ', get_dates
1976  # Some predifinitions
1977  fail = 0
1978  # set filename prefix
1979  prefix = ['grid_conc_', 'grid_pptv_', \
1980  'grid_time_', 'footprint_', 'footprint_total', \
1981  'grid_conc_nest_', 'grid_pptv_nest_', \
1982  'grid_time_nest_', 'footprint_nest_', 'footprint_total_nest'
1983  ]
1984 
1985  units = ['conc', 'pptv', 'time', 'footprint', 'footprint_total']
1986  unit_i = units.index(unit)
1987 
1988  # Determine what module to read, try to use FortFlex, then dumpgrid, lastly pure Python
1989  # import the FortFlex / Fortran module
1990  try:
1991  if OPS.version == 'V6':
1992  from FortFlex import readgrid_v6 as readgrid
1993  from FortFlex import sumgrid
1994  useFortFlex = True
1995  print 'using FortFlex VERSION 6'
1996  else:
1997  #print('Assumed V8 Flexpart')
1998  from FortFlex import readgrid, sumgrid
1999  useFortFlex = True
2000  except:
2001  # get the original module (no memory allocation)
2002  try:
2003  from nilu.pflexpart.FortFlex import readgrid, sumgrid
2004  useFortFlex = True
2005  #print 'using nilu.pflexpart FortFlex'
2006  except:
2007  useFortFlex = False
2008  #print('Cannot load FortFlex, reverting to BinaryFile.')
2009  if not useFortFlex:
2010  readgrid = _readgridBF
2011  OPS.BinaryFile = True
2012 
2013 
2014  # reserve output fields
2015  #print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint
2016 
2017  # -------------------------------------------------
2018 
2019  # # add the requests to the fd object to be returned
2020  OPS.unit = unit
2021  fd.options.update(OPS)
2022 
2023  #--------------------------------------------------
2024  # Loop over all times, given in field H['available_dates']
2025  #--------------------------------------------------
2026 
2027  for date_i in range(len(get_dates)):
2028  datestring = get_dates[date_i]
2029  #print datestring
2030  for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1):A
2031 
2032  FLEXDATA[(s, datestring)] = Structure()
2033  spec_fid = '_' + str(s + 1).zfill(3)
2034 
2035  if unit_i != 4:
2036  filename = os.path.join(H['pathname'], \
2037  prefix[(unit_i) + (H.nested * 5)] + datestring + spec_fid)
2038  H.zdims = H.numzgrid
2039 
2040  else:
2041  # grid total footprint
2042  print "Total footprint"
2043  filename = os.path.join(H['pathname'], \
2044  prefix[(unit_i) + (H.nested * 5)] + spec_fid)
2045  H.zdims = 1
2046 
2047  if os.path.exists(filename):
2048  H.filename = filename
2049  # print 'reading: ' + filename
2050  if OPS.verbose:
2051  print 'with values:'
2052  inputvars = ['filename', 'numxgrid', 'numygrid',
2053  'zdims', 'numpoint', 'nageclass', \
2054  'scaledepo', 'scaleconc',
2055  'decayconstant', 'numpointspec']
2056  for v in inputvars:
2057  print v, " ==> ", H[v]
2058 
2059 
2060  if OPS.BinaryFile:
2061  print("Reading {0} with BinaryFile".format(filename))
2062  gridT, wetgrid, drygrid, itime = _readgridBF(H, filename)
2063  else:
2064  # # Quick fix for Sabine's Ship releases, added nspec_int so that only one
2065  # # field of the nspec dimension is actually read
2066  if OPS.npspec_int is not False:
2067  npspec_int = OPS.npspec_int
2068  numpointspec = 1
2069  else:
2070  npspec_int = 0
2071  numpointspec = H.numpointspec
2072 
2073  gridT, wetgrid, drygrid, itime = readgrid(filename, \
2074  H.numxgrid, H.numygrid,
2075  H.zdims, numpointspec, H.nageclass, \
2076  OPS.scaledepo, OPS.scaleconc, H.decayconstant, npspec_int)
2077 
2078  if OPS.getwet:
2079  wet = wetgrid.squeeze()
2080  if OPS.getdry:
2081  dry = drygrid.squeeze()
2082  if forward:
2083  zplot = gridT[:, :, :, :, 0]
2084  else:
2085  zplot = gridT[:, :, :, :, 0]
2086 
2087  if OPS.calcfoot:
2088 
2089  zplot = sumgrid(zplot, gridT, \
2090  H.area, H.Heightnn)
2091 
2092 
2093  # # get the total column and prep the grid
2094  if H.direction == 'forward':
2095  # not trying to do anything here... must be done
2096  # after retrieving the grid
2097  # D = get_slabs(H,np.squeeze(zplot))
2098  rel_i = 0 #H.available_dates.index(datestring)
2099  D = zplot
2100 
2101  else:
2102  D = zplot
2103  rel_i = 'k'
2104 
2105  # # NOTE:
2106  # # If you're changing things here, you might want to change
2107  # # them in fill_backward as well, yes I know... something is
2108  # # poorly designed ;(
2109 
2110  FLEXDATA[(s, datestring)]['grid'] = D # zplot
2111 
2112  FLEXDATA[(s, datestring)]['itime'] = itime
2113 
2114  FLEXDATA[(s, datestring)]['shape'] = zplot.shape
2115 
2116  FLEXDATA[(s, datestring)]['max'] = zplot.max()
2117 
2118  FLEXDATA[(s, datestring)]['min'] = zplot.min()
2119  FLEXDATA[(s, datestring)]['timestamp'] = datetime.datetime.strptime(datestring, '%Y%m%d%H%M%S')
2120  FLEXDATA[(s, datestring)]['species'] = H['species'][s]
2121  FLEXDATA[(s, datestring)]['gridfile'] = filename
2122  FLEXDATA[(s, datestring)]['rel_i'] = rel_i
2123  FLEXDATA[(s, datestring)]['spec_i'] = s
2124  if OPS.getwet:
2125  FLEXDATA[(s, datestring)]['wet'] = wet
2126  else:
2127  FLEXDATA[(s, datestring)]['wet'] = None
2128 
2129  if OPS.getdry:
2130  FLEXDATA[(s, datestring)]['dry'] = dry
2131  else:
2132  FLEXDATA[(s, datestring)]['dry'] = None
2133 
2134  else:
2135  _shout('***ERROR: file %s not found! \n' % filename)
2136  fail = 1
2137 
2138  fd.set_with_dict(FLEXDATA)
2139  try:
2140  # just for testing, set the first available grid as a shortcut
2141  # this will be removed.
2142  qind = (nspec_ret[0], fd.grid_dates[0])
2143  fd.grid = fd[qind][fd[qind].keys()[0]].grid
2144  except:
2145  pass
2146 
2147 
2148  return fd

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.pflexible.save_spectrum (   outf,
  H,
  agespectra,
  spectype = 'agespec',
  header = '## AGECLASS File' 
) 
Save an ageclass or continents spectrum to an outfile.

Definition at line 697 of file pflexible.py.

697 
698  header='## AGECLASS File'):
699  """ Save an ageclass or continents spectrum to an outfile. """
700 
701  if spectype == 'agespec':
702  # ftype = 'AGECLASS'
703  try:
704  T = H.releasetimes
705  spectrum = agespectra
706  nClasses = spectrum.shape[1]
707  except:
708  # assume H is a list or None
709  if H:
710  nClasses = H[0]
711  T = H[1]
712  spectrum = agespectra
713  else:
714  T = agespectra[:, 0]
715  nClasses = agespectra.shape[1] - 1
716  spectrum = agespectra[:, 1:]
717  elif spectype == 'contspec':
718  # assume H is a list or None
719  # ftype = 'SPECTRUM'
720  try:
721  T = H.releasetimes
722  spectrum = agespectra
723  nClasses = spectrum.shape[1]
724  header = '## Continental Spectrum File'
725  except:
726  # assume H is a list or None
727  if H:
728  nClasses = H[0]
729  T = H[1]
730  spectrum = agespectra
731  else:
732  T = agespectra[:, 0]
733  nClasses = agespectra.shape[1] - 1
734  spectrum = agespectra[:, 1:]
735 
736  T = np.array(T)
737  T = np.reshape(T, (len(T), 1))
738  xout = np.hstack((T, spectrum))
739  fmt = '%s ' + nClasses * '%10.5f '
740  outf.write("%s \n" % (header))
741  np.savetxt(outf, xout, fmt=fmt)
742  outf.close()
743 
744 
745 

Here is the caller graph for this function:

Variable Documentation

tuple flextest.flexread.pflexible.__path__ = os.path.abspath(os.curdir)

Definition at line 85 of file pflexible.py.

string flextest.flexread.pflexible.__version__ = '0.9.5'

Definition at line 84 of file pflexible.py.

string flextest.flexread.pflexible.default = 'verbose output'

Definition at line 5166 of file pflexible.py.

int flextest.flexread.pflexible.exit_code = main()

LateX #####.

Definition at line 5184 of file pflexible.py.

tuple flextest.flexread.pflexible.parser
Initial value:
1 = optparse.OptionParser(
2  formatter=optparse.TitledHelpFormatter(),
3  usage=globals()['__doc__'],
4  version='$Id: py.tpl 327 2008-10-05 23:38:32Z root $')

Definition at line 5161 of file pflexible.py.

flextest.flexread.pflexible.plot_footprint = plot_at_level

Definition at line 3681 of file pflexible.py.

flextest.flexread.pflexible.readgridBF = _readgrid_noFF

Definition at line 1608 of file pflexible.py.

flextest.flexread.pflexible.readheader = read_header

Definition at line 1271 of file pflexible.py.

flextest.flexread.pflexible.readheaderV6 = read_header

Definition at line 1273 of file pflexible.py.

flextest.flexread.pflexible.readheaderV8 = read_header

Definition at line 1272 of file pflexible.py.

tuple flextest.flexread.pflexible.start_time = time.time()

Definition at line 5160 of file pflexible.py.

Generated on Sun Mar 27 2016 01:16:02 for FLEXPART Testing Environment CTBTO WO8 by   doxygen 1.8.6