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

Classes

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

Functions

def gridarea
 Functions for reading FLEXPART output #####. More...
 
def _get_header_version
 
def read_header
 
def closest
 

Variables

 readheader = read_header
 
 readheaderV8 = read_header
 
 readheaderV6 = read_header
 

Function Documentation

def flextest.flexread.read_header._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 84 of file read_header.py.

84 
85 def _get_header_version(bf):
86  """
87  Open and read the binary file (bf) header only to the point of
88  the Flexpart version string, return the string, seek to the
89  start of the file.
90 
91  """
92  try:
93  bf = BinaryFile(bf)
94  except:
95  bf = bf
96 
97  ret = bf.tell()
98  bf.seek(12) # start of version string
99  version = bf.read('13S')
100  bf.seek(ret)
101 
102  return version
flextest.flexread.read_header.BinaryFile
HELPER FUNCTIONS ##########.
Definition: read_header.py:569
flextest.flexread.read_header._get_header_version
def _get_header_version
Definition: read_header.py:84

Here is the caller graph for this function:

def flextest.flexread.read_header.closest (   num,
  numlist 
) 
returns the index of the *closest* value in a list

Definition at line 733 of file read_header.py.

733 
734 def closest(num, numlist):
735  """ returns the index of the *closest* value in a list """
736  # check if we're using datetimes
737  dates = False
738  if isinstance(num, datetime.datetime):
739  dates = True
740  if dates:
741  num = date2num(num)
742  assert isinstance(numlist[0], datetime.datetime), \
743  "num is date, numlist must be a list of dates"
744  numlist = date2num(numlist)
745 
746  return (np.abs(numlist - num)).argmin()
747 
748 
flextest.flexread.read_header.closest
def closest
Definition: read_header.py:733
def flextest.flexread.read_header.gridarea (   H)

Functions for reading FLEXPART output #####. 

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 36 of file read_header.py.

36 
37 def gridarea(H):
38  """returns an array of area corresponding to each nx,ny,nz
39 
40  Usage::
41 
42  > area = gridarea(H)
43 
44 
45  Returns
46  OUT = array area corresponding to nx,ny,nz
47 
48  Arguments
49  H = :class:`Header` object from readheader function.
50 
51  """
52 
53 
54  pih = pi / 180.
55  r_earth = 6.371e6
56  cosfunc = lambda y : cos(y * pih) * r_earth
57 
58  nx = H['numxgrid']
59  ny = H['numygrid']
60  outlat0 = H['outlat0']
61  dyout = H['dyout']
62  dxout = H['dxout']
63  area = np.zeros((nx, ny))
64 
65  for iy in range(ny):
66  ylata = outlat0 + (float(iy) + 0.5) * dyout # NEED TO Check this, iy since arrays are 0-index
67  ylatp = ylata + 0.5 * dyout
68  ylatm = ylata - 0.5 * dyout
69  if (ylatm < 0 and ylatp > 0): hzone = dyout * r_earth * pih
70  else:
71  # cosfact = cosfunc(ylata)
72  cosfactp = cosfunc(ylatp)
73  cosfactm = cosfunc(ylatm)
74  if cosfactp < cosfactm:
75  hzone = sqrt(r_earth ** 2 - cosfactp ** 2) - sqrt(r_earth ** 2 - cosfactm ** 2)
76  else:
77  hzone = sqrt(r_earth ** 2 - cosfactm ** 2) - sqrt(r_earth ** 2 - cosfactp ** 2)
78 
79  gridarea = 2.*pi * r_earth * hzone * dxout / 360.
80  for ix in range(nx):
81  area[ix, iy] = gridarea
82 
83  return area
flextest.flexread.read_header.gridarea
def gridarea
Functions for reading FLEXPART output #####.
Definition: read_header.py:36

Here is the caller graph for this function:

def flextest.flexread.read_header.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 103 of file read_header.py.

103 
104 def read_header(pathname, **kwargs):
105  """
106  The readheader function returns a special class (Structure) which behaves
107  like a dictionary. It contains all the metadata from the simulation which
108  is contained in the "header" or "header_nest" binary files from the model
109  output.
110 
111  .. warning::
112  It is recommended to use the :class:`Header` class: H = pf.Header(path)
113 
114  This version is using the BinaryFile class rather than FortFlex.
115 
116  Usage::
117 
118  > H = read_header(pathname) #Don't include header filename
119 
120  Returns a dictionary
121 
122  H = dictionary like object with all the run metadata.
123  Arguments
124 
125  .. tabularcolumns:: |l|L|
126 
127  ============= ========================================
128  keyword Description [default]
129  ============= ========================================
130  pathname FLEXPART run output directory
131  readp read release points 0=no, [1]=y
132  nested nested output [False] or True
133  headerfile provide custom name for header file
134  datefile provide a custom name for the date file
135  verbose print information while loading header
136  ============= ========================================
137 
138  .. note::
139  **This function is in development**
140 
141  Please report any bugs found.
142 
143  .. TODO::
144 
145  probably a lot of things, among which ...
146 
147  [] choose skip/getbin or direct seek/read
148  [] define output keys in docstring
149 
150 
151  .. note::
152  The user is no longer required to indicate which version of FLEXPART
153  the header is from. Checks are made, and the header is read accordingly.
154  Please report any problems...
155 
156  """
157 
158  OPS = Structure()
159  OPS.readp = True
160  OPS.nested = False
161  OPS.ltopo = 1 # 1 for AGL, 0 for ASL
162  OPS.verbose = False
163  OPS.headerfile = None
164  OPS.datefile = None
165 
166  # # add keyword overides and options to header
167  for k in kwargs.keys():
168  if k not in OPS.keys():
169  print("WARNING: {0} not a valid input option.".format(k))
170 
171  # BW compat fixes
172  if 'nest' in kwargs.keys():
173  raise IOError("nest is no longer a valid keyword, see docs. \n Now use nested=True or nested=False")
174 
175 
176  if 'nested' in kwargs.keys():
177  if kwargs['nested'] is 1:
178  print("WARNING, use of nested=1, deprecated converting to nested=True")
179  kwargs['nested'] = True
180 
181  if 'nested' in kwargs.keys():
182  if kwargs['nested'] is 0:
183  print("WARNING, use of nested=0, deprecated converting to nested=False")
184  kwargs['nested'] = False
185 
186  OPS.update(kwargs)
187 
188  if OPS.verbose:
189  print "Reading Header with:\n"
190 
191  for o in OPS:
192  print "%s ==> %s" % (o, OPS[o])
193 
194  # Define utility functions for reading binary file
195  skip = lambda n = 8 : bf.seek(n, 1)
196  getbin = lambda dtype, n = 1 : bf.read(dtype, (n,))
197 
198 
199  h = Structure()
200 
201 
202  if OPS.headerfile:
203  filename = os.path.join(pathname, OPS.headerfile);
204 
205  elif OPS.nested is True:
206  filename = os.path.join(pathname, 'header_nest')
207  h['nested'] = True;
208  else:
209  filename = os.path.join(pathname, 'header')
210  h['nested'] = False;
211 
212  # Open header file in binary format
213  if not os.path.exists(filename):
214  raise IOError("No such file: {0}".format(filename))
215  else:
216  try:
217  bf = BinaryFile(filename, order="fortran")
218  except:
219  raise IOError("Error opening: {0} with BinaryFile class".format(filename))
220 
221  # Get available_dates from dates file in same directory as header
222  if OPS.datefile:
223  datefile = os.path.join(pathname, OPS.datefile)
224  else:
225  datefile = os.path.join(pathname, 'dates')
226 
227  if not os.path.exists(datefile):
228  raise IOError("No DATEFILE: {0}".format(datefile))
229  else:
230  try:
231  fd = file(datefile, 'r').readlines()
232  except:
233  raise IOError("Could not read datefile: {0}".format(datefile))
234 
235  # get rid of any duplicate dates (a fix for the forecast system)
236  fd = sorted(list(set(fd)))
237  h['available_dates'] = [d.strip('\n') for d in fd]
238 
239 
240  # which version format is header file:
241  version = _get_header_version(bf)
242 
243 
244  # required containers
245  junk = [] # for catching unused output
246  h['nz_list'] = []
247  h['species'] = []
248  h['wetdep'] = []
249  h['drydep'] = []
250  h['ireleasestart'] = []
251  h['ireleaseend'] = []
252  h['compoint'] = []
253 
254  # Define Header format and create Dictionary Keys
255  I = {0:'_0', 1:'ibdate', 2:'ibtime', 3:'flexpart', \
256  4:'_1', 5:'loutstep', 6:'loutaver', 7:'loutsample', \
257  8:'_2', 9:'outlon0', 10:'outlat0', 11:'numxgrid', \
258  12:'numygrid', 13:'dxout', 14:'dyout', 15:'_3', 16:'numzgrid', \
259  }
260  # format for binary reading first part of the header file
261  Dfmt = ['i', 'i', 'i', '13S', '2i', 'i', 'i', 'i', '2i', 'f', 'f', 'i', 'i', 'f', 'f', '2i', 'i']
262  if bf:
263  a = [bf.read(fmt) for fmt in Dfmt]
264  for j in range(len(a)):
265  h[I[j]] = a[j]
266 
267  h['outheight'] = np.array([bf.read('f') for i in range(h['numzgrid'])])
268  junk.append(bf.read('2i'))
269  h['jjjjmmdd'] = bf.read('i')
270  h['hhmmss'] = bf.read('i')
271  junk.append(bf.read('2i'))
272  h['nspec'] = bf.read('i') / 3 # why!?
273  h['numpointspec'] = bf.read('i')
274  junk.append(bf.read('2i'))
275 
276  # Read in the species names and levels for each nspec
277  # temp dictionaries
278  for i in range(h['nspec']):
279  if 'V6' in version:
280 
281  h['wetdep'].append(''.join([bf.read('c') for i in range(10)]).strip())
282  junk.append(bf.read('2i'))
283  junk.append(bf.read('i'))
284  h['drydep'].append(''.join([bf.read('c') for i in range(10)]).strip())
285  junk.append(bf.read('2i'))
286  h['nz_list'].append(getbin('i'))
287  h['species'].append(''.join([getbin('c') for i in range(10)]).strip())
288 
289  else:
290 
291  junk.append(bf.read('i'))
292  h['wetdep'].append(''.join([bf.read('c') for i in range(10)]).strip())
293  junk.append(bf.read('2i'))
294  junk.append(bf.read('i'))
295  h['drydep'].append(''.join([bf.read('c') for i in range(10)]).strip())
296  junk.append(bf.read('2i'))
297  h['nz_list'].append(bf.read('i'))
298  h['species'].append(''.join([bf.read('c') for i in range(10)]).strip())
299  junk.append(bf.read('2i'))
300 
301  if 'V6' in version:
302  bf.seek(8, 1)
303  # pdb.set_trace()
304  h['numpoint'] = bf.read('i')
305 
306  # read release info if requested
307  # if OPS.readp pass has changed, we cycle through once,
308  # then break the loop if OPS.readp is false. This is done
309  # in order to get some of the information into the header.
310  before_readp = bf.tell()
311 
312  # initialise release fields
313  I = {2:'kindz', 3:'xp1', 4:'yp1', 5:'xp2', \
314  6:'yp2', 7:'zpoint1', 8:'zpoint2', 9:'npart', 10:'mpart'}
315 
316  for k, v in I.iteritems():
317  h[v] = np.zeros(h['numpoint']) # create zero-filled lists in H dict
318 
319  h['xmass'] = np.zeros((h['numpoint'], h['nspec']))
320 
321 
322  if 'V6' in version:
323  skip()
324  _readV6(bf, h)
325  else:
326  junk.append(bf.read('i'))
327  for i in range(h['numpoint']):
328  junk.append(bf.read('i'))
329 
330  h['ireleasestart'].append(bf.read('i'))
331  h['ireleaseend'].append(bf.read('i'))
332  h['kindz'][i] = bf.read('h') # This is an int16, might need to to change something
333  junk.append(bf.read('2i'))
334 
335  h['xp1'][i] = bf.read('f')
336  h['yp1'][i] = bf.read('f')
337  h['xp2'][i] = bf.read('f')
338  h['yp2'][i] = bf.read('f')
339  h['zpoint1'][i] = bf.read('f')
340  h['zpoint2'][i] = bf.read('f')
341 
342  junk.append(bf.read('2i'))
343  h['npart'][i] = bf.read('i')
344  h['mpart'][i] = bf.read('i')
345 
346  junk.append(bf.read('i'))
347  # initialise release fields
348 
349  l = bf.read('i') # get compoint length?
350  gt = bf.tell() + l # create 'goto' point
351  sp = ''
352  while re.search("\w", bf.read('c')): # collect the characters for the compoint
353  bf.seek(-1, 1)
354  sp = sp + bf.read('c')
355 
356  bf.seek(gt) # skip ahead to gt point
357 
358  h['compoint'].append(sp) # species names in dictionary for each nspec
359  # h['compoint'].append(''.join([bf.read('c') for i in range(45)]))
360 
361  junk.append(bf.read('i'))
362  # now loop for nspec to get xmass
363  for v in range(h['nspec']):
364  Dfmt = ['i', 'f', '2i', 'f', '2i', 'f', 'i']
365  a = [bf.read(fmt) for fmt in Dfmt]
366  h['xmass'][i, v] = a[1]
367 
368  if OPS.readp is False:
369  """
370  We get the first set of release points here in order
371  to get some information, but skip reading the rest
372  """
373  bf.seek(before_readp)
374  if 'V6' in version:
375  bf.seek(157 * h['numpoint'], 1);
376  else:
377  bf.seek(119 * h['numpoint'] + (h['nspec'] * 36) * h['numpoint'] + 4, 1)
378  break
379 
380 
381 
382  if 'V6' in version:
383  h['method'] = bf.read('i')
384 
385  else:
386  junk.append(bf.read('i'))
387  junk.append(bf.read('i'))
388 
389  h['lsubgrid'] = bf.read('i')
390  h['lconvection'] = bf.read('i')
391  h['ind_source'] = bf.read('i')
392  h['ind_receptor'] = bf.read('i')
393 
394  if 'V6' in version:
395  pass
396  else:
397  junk.append(bf.read('2i'))
398 
399  h['nageclass'] = bf.read('i')
400 
401  Lage_fmt = ['i'] * h.nageclass
402  h['lage'] = [bf.read(fmt) for fmt in Lage_fmt]
403 
404  # Orography
405  nx = h['numxgrid']
406  ny = h['numygrid']
407  Dfmt = ['f'] * nx
408  h['oro'] = np.zeros((nx, ny), np.float)
409  junk.append(bf.read('2i'))
410 
411  for ix in range(nx):
412  h['oro'][ix] = bf.read('f', ny)
413  bf.seek(8, 1)
414 
415  # Why was this? / deprecated.
416  # if h['loutstep'] < 0:
417  # h['nspec'] = h['numpoint']
418 
419  bf.close()
420 
421 
422 
423 
424  ############# ADDITIONS ###########
425  # attributes to header that can be added after reading below here
426  h['pathname'] = pathname
427  h['decayconstant'] = 0
428  h.path = pathname
429 
430 
431  # Calculate Height (outheight + topography)
432  # There is an offset issue here related to the 0-indexing. Be careful.
433  oro = h['oro'] # z is a numpy array
434  nx = h['numxgrid']
435  ny = h['numygrid']
436  nz = h['numzgrid']
437 
438  Heightnn = np.zeros((nx, ny, nz), np.float)
439  for ix in range(nx):
440  if OPS.ltopo == 1:
441  Heightnn[ix, :, 0] = oro[ix, :]
442  else:
443  Heightnn[ix, :, 0] = np.zeros(ny)
444 
445  for iz in range(nz):
446  if OPS.ltopo == 1:
447  Heightnn[ix, :, iz] = [h['outheight'][iz] + oro[ix, y] for y in range(ny)]
448  else:
449  Heightnn[ix, :, iz] = h['outheight'][iz]
450 
451  h['Area'] = gridarea(h)
452  h['Heightnn'] = Heightnn
453  h['nx'] = nx
454  h['ny'] = ny
455  h['nz'] = nz
456 
457  # Convert ireleasestart and ireleaseend to datetimes, fix issue #10
458  start_day = datetime.datetime.strptime(str(h['ibdate']), '%Y%m%d')
459  H, M, S , = [int(str(h['ibtime']).zfill(6)[i:i + 2]) for i in range(0, 6, 2)]
460  start_time = datetime.timedelta(hours=H, minutes=M, seconds=S)
461  h['simulationstart'] = start_day + start_time
462 
463  if OPS.readp:
464  releasestart, releaseend = [], []
465  for i in range(h.numpointspec):
466  releasestart.append(h.simulationstart + \
467  datetime.timedelta(seconds=int(h.ireleasestart[i])))
468  releaseend.append(h.simulationstart + \
469  datetime.timedelta(seconds=int(h.ireleaseend[i])))
470  h.releasestart = releasestart
471  h.releaseend = releaseend[:h.numpointspec]
472  h.releasetimes = [b - ((b - a) / 2) for a, b in zip(h.releasestart, h.releaseend)]
473 
474  # Add datetime objects for dates
475  available_dates_dt = []
476  for i in h.available_dates:
477  available_dates_dt.append(datetime.datetime(
478  int(i[:4]), int(i[4:6]), int(i[6:8]), int(i[8:10]), int(i[10:12]), int(i[12:])))
479  h.available_dates_dt = available_dates_dt
480  h.first_date = available_dates_dt[0]
481  h.last_date = available_dates_dt[-1]
482  h.ageclasses = np.array([act - h.simulationstart for act in h.available_dates_dt])
483  h.numageclasses = len(h.ageclasses)
484 
485  # Add other helpful attributes
486  h.nxmax = h.numxgrid
487  h.nymax = h.numygrid
488  h.nzmax = h.numzgrid
489  h.maxspec = h.nspec
490  h.maxpoint = h.numpoint
491  h.maxageclass = h.numageclasses
492 
493  h.area = h.Area # fix an annoyance
494 
495  if OPS.readp:
496  h.xpoint = h.xp1
497  h.ypoint = h.yp1
498 
499  # Add release unit derived from kindz
500  if 'kindz' not in h.keys():
501  h.kindz = [0]
502  h.alt_unit = 'unkn.'
503  if 3 in h.kindz:
504  h.alt_unit = 'hPa'
505  elif 2 in h.kindz:
506  h.alt_unit = 'm.a.s.l.'
507  elif 1 in h.kindz:
508  h.alt_unit = 'm.a.g.l.'
509 
510  #
511  if h.loutstep > 0:
512  h.direction = 'forward'
513  h.unit = 'conc' # could be pptv
514  h.plot_unit = 'ppb'
515  else:
516  h.direction = 'backward'
517  h.unit = 'time'
518  h.plot_unit = 'ns / kg' # Not sure about this
519 
520  # Units based on Table 1, ACP 2005
521  if h.direction == 'forward':
522  if h.ind_source == 1:
523  if h.ind_receptor == 1:
524  h.output_unit = 'ng m-3'
525  if h.ind_receptor == 2:
526  h.output_unit = 'pptm'
527  if h.ind_source == 2:
528  if h.ind_receptor == 1:
529  h.output_unit = 'ng m-3'
530  if h.ind_receptor == 2:
531  h.output_unit = 'pptm'
532  if h.direction == 'backward':
533  if h.ind_source == 1:
534  if h.ind_receptor == 1:
535  h.output_unit = 's'
536  if h.ind_receptor == 2:
537  h.output_unit = 's m^3 kg-1'
538  if h.ind_source == 2:
539  if h.ind_receptor == 1:
540  h.output_unit = 's kg m-3'
541  if h.ind_receptor == 2:
542  h.output_unit = 's'
543 
544 
545 
546  # Add layer thickness
547  layerthickness = [h.outheight[0]]
548  for i, lh in enumerate(h.outheight[1:]):
549  layerthickness.append(lh - h.outheight[i])
550  h.layerthickness = layerthickness
551 
552  h.options = OPS
553  h.fp_version = h.flexpart
554  h.junk = junk # the extra bits that were read... more for debugging
555 
556 
557  #print('Read {0} Header: {1}'.format(h.fp_version, filename))
558 
559  return h
560 
# BW Compatability
flextest.flexread.read_header.BinaryFile
HELPER FUNCTIONS ##########.
Definition: read_header.py:569
flextest.flexread.read_header._get_header_version
def _get_header_version
Definition: read_header.py:84
flextest.flexread.read_header.gridarea
def gridarea
Functions for reading FLEXPART output #####.
Definition: read_header.py:36
flextest.flexread.read_header.Structure
Definition: read_header.py:749
flextest.flexread.read_header.read_header
def read_header
Definition: read_header.py:103

Here is the call graph for this function:

Variable Documentation

flextest.flexread.read_header.readheader = read_header

Definition at line 561 of file read_header.py.

flextest.flexread.read_header.readheaderV6 = read_header

Definition at line 563 of file read_header.py.

flextest.flexread.read_header.readheaderV8 = read_header

Definition at line 562 of file read_header.py.

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