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

Classes

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

Functions

def _readgridBF
 
def read_grid
 
def readgridV8
 
def readgridV6
 
def closest
 
def data_range
 
def main
 Below for use from command line ##########################################. More...
 

Function Documentation

def flextest.flexread.read_grid._readgridBF (   H,
  filename 
)
private 
Read grid using BinaryFile class

Definition at line 31 of file read_grid.py.

31 
32 def _readgridBF(H, filename):
33  """ Read grid using BinaryFile class"""
34  # Utility functions
35  skip = lambda n = 8 : f2.seek(n, 1)
36  getbin = lambda dtype, n = 1 : f2.read(dtype, (n,))
37 
38 
39  def getdump(n, fmt='f'):
40  """ function to get the dump values for the sparse format """
41  skip()
42  # Dfmt=[fmt]*n
43 # a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt]
44  a = f2.read(fmt, n)
45  ### DJM HACK - 2013-09-02 - without the following, 'a' ends up being
46  ### a float rather than a numpy array
47  if n==1:
48  a = np.array([a])
49 # dumplist=[a[j][0] for j in range(len(a))]
50  # dumplist=[a[j] for j in range(len(a))]
51  return a # dumplist
52 
53  def key2var(D, key):
54  cmd = "global %s; %s = D['%s'];" % (key, key, key)
55  exec(cmd)
56 
57  def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny):
58  """ function to dump sparse elements into grid, fall back method if
59  pflexcy.so module (cython) not available -- it is SLOW """
60  conc = False
61  if len(grd.shape) == 5:
62  conc = True
63  ii = 0
64  fact = 1
65  pos = 0
66  '''
67  print 'CNT_R: ' + str(cnt_r)
68  print 'DMP_R: ' + str(dmp_r)
69  print 'TYPE DMP_R: ' + str(type(dmp_r))
70  '''
71  for ir in range(cnt_r):
72 
73  if conc:
74  if dmp_r[ir] * fact > 0:
75  n = dmp_i[ii]
76  ii = ii + 1
77  fact = fact * -1.
78  else:
79  n = n + 1
80 
81  kz = n / (H.numxgrid * H.numygrid)
82  jy = (n - kz * H.numxgrid * H.numygrid) / H.numxgrid
83  ix = n - H.numxgrid * H.numygrid * kz - H.numxgrid * jy
84  grd[ix, jy, kz - 1, k, nage] = abs(dmp_r[ir])
85 
86 #
87 # print "n ==> ix,jy,kz,k,nage"
88 # print "%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage)
89 # print grd.shape
90 # print grd[0,0,0,0,0]
91 
92  else:
93  if dmp_r[ir] * fact > 0:
94  n = dmp_i[ii]
95  ii = ii + 1
96  fact = fact * -1.
97  else:
98  n = n + 1
99  jy = n / H.numxgrid
100  ix = n - H.numxgrid * jy
101  grd[ix, jy, k, nage] = abs(dmp_r[ir])
102 
103  return grd # flipud(grd.transpose())
104 
105  # # Import pflexcy.so (cython compiled version of dumpgrid)
106  try:
107 
108  from sdspflexcy import dumpdatagrid, dumpdepogrid
109  print 'using pflexcy'
110  except:
111  #print """WARNING: Using PURE Python to readgrid, execution will be slow.
112  # Try compiling the FortFlex module or the pflexcy module
113  # for your machine. For more information see the
114  # pflexible/f2py_build directory or use cython with pflexcy.pyx
115  #"""
116  dumpdatagrid = _dumpgrid
117  dumpdepogrid = _dumpgrid
118 
119  dat_cnt = 0
120  nage = 1
121  # #datagrid=np.zeros((numxgrid,numygrid,numzgrid[nspec-1],nspec,nageclass),np.float)
122  wetgrid = np.zeros((H.numxgrid, H.numygrid, H.numpointspec, 1), np.float)
123  drygrid = np.zeros((H.numxgrid, H.numygrid, H.numpointspec, 1), np.float)
124  datagrid = np.zeros((H.numxgrid, H.numygrid, H.numzgrid, H.numpointspec, nage), np.float)
125  # f = file(filename,'rb')
126  # print filename
127  f2 = BinaryFile(filename, order='fortran')
128  # read data:
129  skip(4)
130  itime = getbin('i')
131 
132  for na in range(nage):
133 
134  for ks in range(H.numpointspec):
135  # Read Wet Depostion
136  skip()
137  cnt_i = getbin('i')
138  dmp_i = getdump(cnt_i, 'i')
139  skip()
140  cnt_r = getbin('i')
141  dmp_r = getdump(cnt_r)
142  if dmp_r.any():
143  # print dmp_r, dmp_i
144  wetgrid = dumpdepogrid(dmp_i, cnt_r, dmp_r, wetgrid, ks, na, H.numxgrid, H.numygrid)
145 
146  # Read Dry Deposition
147  skip()
148  cnt_i = getbin('i')
149  dmp_i = getdump(cnt_i, 'i')
150  skip()
151  cnt_r = getbin('i')
152  dmp_r = getdump(cnt_r)
153  if dmp_r.any():
154  # print dmp_r, dmp_i
155  drygrid = dumpdepogrid(dmp_i, cnt_r, dmp_r, drygrid, ks, na, H.numxgrid, H.numygrid)
156 
157  # Read Concentrations
158  skip()
159  cnt_i = getbin('i')
160  dmp_i = getdump(cnt_i, 'i')
161  skip()
162  cnt_r = getbin('i')
163  dmp_r = getdump(cnt_r)
164  # print len(dmp_r),len(dmp_i)
165  # print cnt_r,cnt_i
166  # print dmp_i
167  # print type(dmp_i),type(cnt_r),type(dmp_r),type(datagrid),type(ks),type(na)
168  # print type(H.numxgrid),type(H.numygrid)
169  datagrid = dumpdatagrid(dmp_i, cnt_r, dmp_r, datagrid, ks, na, H.numxgrid, H.numygrid)
170 
171  # G[H['ibtime']].append(concgrid)
172  f2.close()
173 
174  return datagrid, wetgrid, drygrid, itime
flextest.flexread.read_grid._readgridBF
def _readgridBF
Definition: read_grid.py:31
flextest.flexread.read_grid.BinaryFile
HELPER FUNCTIONS ##########.
Definition: read_grid.py:742

Here is the caller graph for this function:

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

Definition at line 909 of file read_grid.py.

909 
910 def closest(num, numlist):
911  """ returns the index of the *closest* value in a list """
912  # check if we're using datetimes
913  dates = False
914  if isinstance(num, datetime.datetime):
915  dates = True
916  if dates:
917  num = date2num(num)
918  assert isinstance(numlist[0], datetime.datetime), \
919  "num is date, numlist must be a list of dates"
920  numlist = date2num(numlist)
921 
922  return (np.abs(numlist - num)).argmin()
923 
924 
flextest.flexread.read_grid.closest
def closest
Definition: read_grid.py:909
def flextest.flexread.read_grid.data_range (   data,
  min = 'median' 
) 
return a data range for flexpart data

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

Definition at line 952 of file read_grid.py.

952 
953 def data_range(data, min='median'):
954  """
955  return a data range for flexpart data
956 
957  optional keyword min = ['median', 'mean', 'min']
958  """
959  dmax = np.nanmax(data)
960  if np.isnan(dmax):
961  dmax = 1e5
962 
963  if min == 'mean':
964  dmin = np.mean(data[data.nonzero()])
965  elif min == 'median':
966  dmin = np.median(data[data.nonzero()])
967  else:
968  dmin = np.nanmin(data[data.nonzero()])
969 
970  if np.isnan(dmin):
971  dmin = 1e-5
972 
973  return [dmin, dmax]
flextest.flexread.read_grid.data_range
def data_range
Definition: read_grid.py:952
def flextest.flexread.read_grid.main ( )

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

Definition at line 976 of file read_grid.py.

976 
977 def main ():
978 
979  global options, args
980  here = os.path.curdir
981  dataDir = options.pathname
982  H = read_header(dataDir, **{'readp':1}) # readheader and release pointspiu
983  G = readgridV8(H, **{'unit':'time'})
flextest.flexread.read_grid.readgridV8
def readgridV8
Definition: read_grid.py:247
flextest.flexread.read_grid.main
def main
Below for use from command line ##########################################.
Definition: read_grid.py:976

Here is the call graph for this function:

def flextest.flexread.read_grid.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 175 of file read_grid.py.

175 
176 def read_grid(H, **kwargs):
177  """
178  Accepts a header object as input, returns dictionary of Grid values
179  keyed by species and datestring from H['available_dates']. A grid instance from this
180  dictionary may be passed to the get_slabs function to return a dictionary of slabs
181  suitable for plotting with the plotting routines. See below for more information.
182 
183  **DEPENDENCY**
184  Requires FortFlex.so module compiled using f2py. See FortFlex.f for more details.
185 
186  Usage::
187 
188  > FLEXDATA = read_grid(H,**kwargs)
189 
190  Returns:
191 
192  A grid dictionary key by tuples with (species,available_dates), where species is
193  and integer. See grid.keys()
194 
195  FLEXDATA[(s,datestring)]['grid']
196  FLEXDATA[(s,datestring)]['itime']
197  FLEXDATA[(s,datestring)]['shape']
198  FLEXDATA[(s,datestring)]['max']
199  FLEXDATA[(s,datestring)]['min']
200  FLEXDATA[(s,datestring)]['timestamp']
201  FLEXDATA[(s,datestring)]['species']
202  FLEXDATA[(s,datestring)]['gridfile']
203  FLEXDATA[(s,datestring)]['rel_i']
204  FLEXDATA[(s,datestring)]['spec_i']
205 
206  Arguments
207 
208  .. tabularcolumns:: |l|L|
209 
210  ============= ========================================
211  keyword Description [default]
212  ============= ========================================
213  date which yyyymmddhhmmss from available_dates
214  or use (time_ret)
215  time_ret index to time
216  unit 'conc', 'pptv', ['time'], 'footprint'
217  nspec_ret numspecies
218  pspec_ret index to ???
219  age_ret index to ageclass
220  nested obtained from H['nested']
221  BinaryFile Use BinaryFile vs. FortFlex [False]
222  getwet True, [False]
223  getdry True, [False]
224  scaledepo A float value to scale deposition [1.0]
225  scaleconc A float value to scale conc [1.0]
226  decayconstant A float for decay const. [9e6]
227  calcfoot Will cause footprint to be calculated
228  [False]
229  verbose more output
230  version A string 'V8' or 'V6', ['V8']
231  ============= ========================================
232 
233 
234  .. note::
235  most arguments are able to be extracted fro the header "H"
236 
237  """
238 
239  if H.version == 'V9':
240  return readgridV9(H, **kwargs)
241  if H.version == 'V8':
242  return readgridV8(H, **kwargs)
243  if H.version == 'V6':
244  return readgridV6(H, **kwargs)
245  else:
246  raise IOError("No version attribute defined for Header.")
flextest.flexread.read_grid.readgridV6
def readgridV6
Definition: read_grid.py:498
flextest.flexread.read_grid.readgridV8
def readgridV8
Definition: read_grid.py:247
flextest.flexread.read_grid.read_grid
def read_grid
Definition: read_grid.py:175

Here is the call graph for this function:

def flextest.flexread.read_grid.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 498 of file read_grid.py.

498 
499 def readgridV6(H, **kwargs):
500 
501  """
502  Accepts a header object as input, and selects appropriate readgrid function
503  to use for reading in data from the flexpart binary Fortran files.
504 
505  See the :func:`read_grid` for more information on keyword arguments
506 
507  This is the 'V6' version of the function.
508 
509  """
510  # # OPS is the options Structure, sets defaults, then update w/ kwargs
511  OPS = Structure()
512  OPS.unit = 'time'
513  OPS.nspec_ret = 0
514  OPS.pspec_ret = 0
515  OPS.age_ret = 0
516  OPS.time_ret = 0
517  OPS.scaledepo = 1.0
518  OPS.scaleconc = 1.0
519  OPS.decayconstant = 9e6
520  OPS.date = None
521  OPS.calcfoot = False
522  OPS.verbose = False
523  OPS.BinaryFile = False
524  # # add keyword overides and options to header
525  OPS.update(kwargs)
526  # H.update(OPS)
527 
528  # # set up the return dictionary (FLEXDATA updates fd, fd is returned)
529  FLEXDATA = {}
530  fd = Structure()
531  fd.options = Structure()
532  # # add the requests to the fd object to be returned
533  fd.options.update(OPS)
534 
535  # # What direction is the run?
536  unit = OPS.unit
537  if H['loutstep'] > 0:
538  forward = True
539  if unit == 'time':
540  # # default forward unit
541  unit = 'conc'
542  else:
543  forward = False
544 
545  # # What species to return?
546  nspec_ret = OPS.nspec_ret
547  if isinstance(nspec_ret, int):
548  nspec_ret = [nspec_ret]
549  assert iter(nspec_ret), "nspec_ret must be iterable."
550 
551  # # get times to return
552  get_dates = None
553  if OPS.time_ret is not None:
554  get_dates = []
555  time_ret = OPS.time_ret
556  if isinstance(time_ret, int) == True:
557  time_ret = [time_ret]
558 
559  if time_ret[0] < 0:
560  if forward == False:
561  # # get all dates for calculating footprint.
562  time_ret = np.arange(len(H.available_dates))
563  else:
564  raise ValueError("Must enter a positive time_ret for forward runs")
565 
566  for t in time_ret:
567  get_dates.append(H.available_dates[t])
568 
569 
570  # # define what dates to extract if user has explicitly defined a 'date'
571  if OPS.date != None:
572  date = OPS.date
573  if time_ret is not None:
574  Warning("overwriting time_ret variable, date was requested")
575  get_dates = []
576  if not isinstance(date, list):
577  date = date.strip().split(',')
578  for d in date:
579  try:
580  get_dates.append(H.available_dates[H.available_dates.index(d)])
581  time_ret = None
582  except:
583  _shout("Cannot find date: %s in H['available_dates']\n" % d)
584 
585  if get_dates is None:
586  raise ValueError("Must provide either time_ret or date value.")
587  else:
588  # # assign grid dates for indexing fd
589  fd.grid_dates = get_dates[:]
590 
591  #print 'getting grid for: ', get_dates
592  # Some predifinitions
593  fail = 0
594  # set filename prefix
595  prefix = ['grid_conc_', 'grid_pptv_', \
596  'grid_time_', 'footprint_', 'footprint_total', \
597  'grid_conc_nest_', 'grid_pptv_nest_', \
598  'grid_time_nest_', 'footprint_nest_', 'footprint_total_nest'
599  ]
600 
601  units = ['conc', 'pptv', 'time', 'footprint', 'footprint_total']
602  unit_i = units.index(unit)
603 
604  useFortFlex = False
605  if not useFortFlex:
606  readgrid = _readgridBF
607  OPS.BinaryFile = True
608 
609  # cheat: use key2var function to get values from header dict, H
610  # need to change this to using H.xxxxx
611  # headervars = ['nested','nspec','numxgrid','numygrid','numzgrid','nageclass',\
612  # 'dates','pathname','decayconstant','numpoint','numpointspec',
613  # 'area','Heightnn','lage']
614  # for k in headervars:
615  # key2var(H,k)
616 
617 
618 
619  # reserve output fields
620  #print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint
621 
622  # -------------------------------------------------
623 
624  # if forward:
625  # numpointspec = 1
626  # else:
627  # numpointspec = numpoint
628  # numpointspec = H['numpointspec']
629  # if unit_i == 4:
630  # grid=np.empty((numxgrid,numygrid,1,nspec_ret,pspec_ret,age_ret,len(get_dates)))
631  # else:
632  # grid=np.empty((numxgrid,numygrid,numzgrid,nspec_ret,pspec_ret,age_ret,len(get_dates)))
633  # zplot=np.empty((numxgrid,numygrid,numzgrid,numpointspec))
634 
635  #--------------------------------------------------
636  # Loop over all times, given in field H['available_dates']
637  #--------------------------------------------------
638 
639  for date_i in range(len(get_dates)):
640  datestring = get_dates[date_i]
641  #print datestring
642  FLEXDATA[datestring] = {}
643  for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1):
644  total_footprint = False
645  FLEXDATA[(s, datestring)] = Structure()
646  # spec_fid = '_'+str(s+1).zfill(3)
647 
648  if unit_i != 4:
649  filename = os.path.join(H['pathname'], \
650  prefix[(unit_i) + (H.nested * 5)] + datestring)
651  H.zdims = H.numzgrid
652 
653  else:
654  # grid total footprint
655  print "Total footprint"
656  total_footprint = True
657  filename = os.path.join(H['pathname'], \
658  prefix[(unit_i) + (H.nested * 5)])
659  H.zdims = 1
660 
661  if os.path.exists(filename):
662  H.filename = filename
663  # print 'reading: ' + filename
664  if OPS.verbose:
665  print 'with values:'
666  inputvars = ['filename', 'numxgrid', 'numygrid',
667  'zdims', 'numpoint', 'nageclass', \
668  'scaledepo', 'scaleconc',
669  'decayconstant', 'numpointspec']
670  for v in inputvars:
671  print v, " ==> ", H[v]
672 
673 
674  if OPS.BinaryFile:
675  # print 'Using BinaryFile'
676  gridT, wetgrid, drygrid, itime = _readgridBF(H, filename)
677  else:
678  gridT, wetgrid, drygrid, itime = readgrid(filename, \
679  H.numxgrid, H.numygrid,
680  H.zdims, H.numpointspec, H.nageclass, \
681  OPS.scaledepo, OPS.scaleconc, H.decayconstant)
682 
683 
684  if forward:
685  zplot = gridT[:, :, :, :, 0]
686  else:
687  zplot = gridT[:, :, :, :, 0]
688 
689  # if total_footprint:
690  # zplot = np.squeeze(gridT)
691 
692  if OPS.calcfoot:
693  zplot = sumgrid(zplot, gridT, \
694  H.area, H.Heightnn)
695 
696 
697  # # get the total column and prep the grid
698  if H.direction == 'forward':
699  # not trying to do anything here... must be done
700  # after retrieving the grid
701  # D = get_slabs(H,np.squeeze(zplot))
702  # rel_i = H.available_dates.index(datestring)
703  D = zplot
704  rel_i = 'k'
705  else:
706  D = zplot
707  rel_i = 'k'
708 
709  # # NOTE:
710  # # If you're changing things here, you might want to change
711  # # them in fill_backward as well, yes I know... something is
712  # # poorly designed ;(
713  FLEXDATA[(s, datestring)]['grid'] = D # zplot
714  FLEXDATA[(s, datestring)]['itime'] = itime
715  FLEXDATA[(s, datestring)]['shape'] = zplot.shape
716  FLEXDATA[(s, datestring)]['max'] = zplot.max()
717  FLEXDATA[(s, datestring)]['min'] = zplot.min()
718  FLEXDATA[(s, datestring)]['timestamp'] = datetime.datetime.strptime(datestring, '%Y%m%d%H%M%S')
719  FLEXDATA[(s, datestring)]['species'] = H['species'][s]
720  FLEXDATA[(s, datestring)]['gridfile'] = filename
721  FLEXDATA[(s, datestring)]['rel_i'] = rel_i
722  FLEXDATA[(s, datestring)]['spec_i'] = s
723 
724 
725  else:
726  _shout('***ERROR: file %s not found! \n' % filename)
727  fail = 1
728  fd.set_with_dict(FLEXDATA)
729 
730  try:
731  # just for testing, set the first available grid as a shortcut
732  # this will be removed.
733  qind = (0, fd.grid_dates[0])
734  fd.grid = fd[qind][fd[qind].keys()[0]].grid
735  except:
736  pass
737 
738 
739  return fd
740 
flextest.flexread.read_grid.readgridV6
def readgridV6
Definition: read_grid.py:498
flextest.flexread.read_grid._readgridBF
def _readgridBF
Definition: read_grid.py:31
flextest.flexread.read_grid.Structure
Definition: read_grid.py:925

Here is the call graph for this function:

Here is the caller graph for this function:

def flextest.flexread.read_grid.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 247 of file read_grid.py.

247 
248 def readgridV8(H, **kwargs):
249 
250  """
251  Accepts a header object as input, and selects appropriate readgrid function
252  to use for reading in data from the flexpart binary Fortran files.
253 
254  See the :func:`read_grid` for more information on keyword arguments
255 
256  This is the 'V8' version of the function.
257 
258  """
259  # # OPS is the options Structure, sets defaults, then update w/ kwargs
260  OPS = Structure()
261  OPS.unit = H.unit
262  OPS.getwet = False
263  OPS.getdry = False
264  OPS.nspec_ret = 0
265  OPS.npspec_int = False # allows to select an index of npsec when calling readgrid
266  OPS.pspec_ret = 0
267  OPS.age_ret = 0
268  OPS.time_ret = 0
269  OPS.scaledepo = 1.0
270  OPS.scaleconc = 1.0
271  OPS.decayconstant = 9e6
272  OPS.date = None
273  OPS.calcfoot = False
274  OPS.verbose = False
275  OPS.BinaryFile = False
276  OPS.version = 'V8'
277  # # add keyword overides and options to header
278  OPS.update(kwargs)
279  # H.update(OPS)
280 
281  # # set up the return dictionary (FLEXDATA updates fd, fd is returned)
282  FLEXDATA = {}
283  fd = Structure()
284  fd.options = Structure()
285 
286  # # What direction is the run?
287  unit = OPS.unit
288 
289  if H['loutstep'] > 0:
290  forward = True
291  if unit == 'time':
292  # # default forward unit
293  unit = 'conc'
294  OPS.unit = unit
295  else:
296  forward = False
297 
298  # # What species to return?
299  nspec_ret = OPS.nspec_ret
300  if isinstance(nspec_ret, int):
301  nspec_ret = [nspec_ret]
302  assert iter(nspec_ret), "nspec_ret must be iterable."
303 
304  # # get times to return
305  get_dates = None
306  if OPS.time_ret is not None:
307  get_dates = []
308  time_ret = OPS.time_ret
309  if isinstance(time_ret, int) == True:
310  time_ret = [time_ret]
311 
312  if time_ret[0] < 0:
313  #if forward == False:
314  # # get all dates for calculating footprint.
315  time_ret = np.arange(len(H.available_dates))
316  #else:
317  # raise ValueError("Must enter a positive time_ret for forward runs")
318 
319  for t in time_ret:
320  get_dates.append(H.available_dates[t])
321 
322 
323  # # define what dates to extract if user has explicitly defined a 'date'
324  if OPS.date != None:
325  date = OPS.date
326  if time_ret is not None:
327  Warning("overwriting time_ret variable, date was requested")
328  get_dates = []
329  if not isinstance(date, list):
330  date = date.strip().split(',')
331  for d in date:
332  try:
333  get_dates.append(H.available_dates[H.available_dates.index(d)])
334  time_ret = None
335  except:
336  _shout("Cannot find date: %s in H['available_dates']\n" % d)
337 
338  if get_dates is None:
339  raise ValueError("Must provide either time_ret or date value.")
340  else:
341  # # assign grid dates for indexing fd
342  fd.grid_dates = get_dates[:]
343 
344  #print 'getting grid for: ', get_dates
345  # Some predifinitions
346  fail = 0
347  # set filename prefix
348  prefix = ['grid_conc_', 'grid_pptv_', \
349  'grid_time_', 'footprint_', 'footprint_total', \
350  'grid_conc_nest_', 'grid_pptv_nest_', \
351  'grid_time_nest_', 'footprint_nest_', 'footprint_total_nest'
352  ]
353 
354  units = ['conc', 'pptv', 'time', 'footprint', 'footprint_total']
355  unit_i = units.index(unit)
356 
357  useFortFlex = False
358  if not useFortFlex:
359  readgrid = _readgridBF
360  OPS.BinaryFile = True
361 
362 
363  # reserve output fields
364  #print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint
365 
366  # -------------------------------------------------
367 
368  # # add the requests to the fd object to be returned
369  OPS.unit = unit
370  fd.options.update(OPS)
371 
372  #--------------------------------------------------
373  # Loop over all times, given in field H['available_dates']
374  #--------------------------------------------------
375 
376  for date_i in range(len(get_dates)):
377  datestring = get_dates[date_i]
378  #print datestring
379  for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1):A
380 
381  FLEXDATA[(s, datestring)] = Structure()
382  spec_fid = '_' + str(s + 1).zfill(3)
383 
384  if unit_i != 4:
385  filename = os.path.join(H['pathname'], \
386  prefix[(unit_i) + (H.nested * 5)] + datestring + spec_fid)
387  H.zdims = H.numzgrid
388 
389  else:
390  # grid total footprint
391  print "Total footprint"
392  filename = os.path.join(H['pathname'], \
393  prefix[(unit_i) + (H.nested * 5)] + spec_fid)
394  H.zdims = 1
395 
396  if os.path.exists(filename):
397  H.filename = filename
398  # print 'reading: ' + filename
399  if OPS.verbose:
400  print 'with values:'
401  inputvars = ['filename', 'numxgrid', 'numygrid',
402  'zdims', 'numpoint', 'nageclass', \
403  'scaledepo', 'scaleconc',
404  'decayconstant', 'numpointspec']
405  for v in inputvars:
406  print v, " ==> ", H[v]
407 
408 
409  if OPS.BinaryFile:
410  #print("Reading {0} with BinaryFile".format(filename))
411  gridT, wetgrid, drygrid, itime = _readgridBF(H, filename)
412  else:
413  # # Quick fix for Sabine's Ship releases, added nspec_int so that only one
414  # # field of the nspec dimension is actually read
415  if OPS.npspec_int is not False:
416  npspec_int = OPS.npspec_int
417  numpointspec = 1
418  else:
419  npspec_int = 0
420  numpointspec = H.numpointspec
421 
422  gridT, wetgrid, drygrid, itime = readgrid(filename, \
423  H.numxgrid, H.numygrid,
424  H.zdims, numpointspec, H.nageclass, \
425  OPS.scaledepo, OPS.scaleconc, H.decayconstant, npspec_int)
426 
427  if OPS.getwet:
428  wet = wetgrid.squeeze()
429  if OPS.getdry:
430  dry = drygrid.squeeze()
431  if forward:
432  zplot = gridT[:, :, :, :, 0]
433  else:
434  zplot = gridT[:, :, :, :, 0]
435 
436  if OPS.calcfoot:
437 
438  zplot = sumgrid(zplot, gridT, \
439  H.area, H.Heightnn)
440 
441 
442  # # get the total column and prep the grid
443  if H.direction == 'forward':
444  # not trying to do anything here... must be done
445  # after retrieving the grid
446  # D = get_slabs(H,np.squeeze(zplot))
447  rel_i = 0 #H.available_dates.index(datestring)
448  D = zplot
449 
450  else:
451  D = zplot
452  rel_i = 'k'
453 
454  # # NOTE:
455  # # If you're changing things here, you might want to change
456  # # them in fill_backward as well, yes I know... something is
457  # # poorly designed ;(
458 
459  FLEXDATA[(s, datestring)]['grid'] = D # zplot
460 
461  FLEXDATA[(s, datestring)]['itime'] = itime
462 
463  FLEXDATA[(s, datestring)]['shape'] = zplot.shape
464 
465  FLEXDATA[(s, datestring)]['max'] = zplot.max()
466 
467  FLEXDATA[(s, datestring)]['min'] = zplot.min()
468  FLEXDATA[(s, datestring)]['timestamp'] = datetime.datetime.strptime(datestring, '%Y%m%d%H%M%S')
469  FLEXDATA[(s, datestring)]['species'] = H['species'][s]
470  FLEXDATA[(s, datestring)]['gridfile'] = filename
471  FLEXDATA[(s, datestring)]['rel_i'] = rel_i
472  FLEXDATA[(s, datestring)]['spec_i'] = s
473  if OPS.getwet:
474  FLEXDATA[(s, datestring)]['wet'] = wet
475  else:
476  FLEXDATA[(s, datestring)]['wet'] = None
477 
478  if OPS.getdry:
479  FLEXDATA[(s, datestring)]['dry'] = dry
480  else:
481  FLEXDATA[(s, datestring)]['dry'] = None
482 
483  else:
484  _shout('***ERROR: file %s not found! \n' % filename)
485  fail = 1
486 
487  fd.set_with_dict(FLEXDATA)
488  try:
489  # just for testing, set the first available grid as a shortcut
490  # this will be removed.
491  qind = (nspec_ret[0], fd.grid_dates[0])
492  fd.grid = fd[qind][fd[qind].keys()[0]].grid
493  except:
494  pass
495 
496 
497  return fd
flextest.flexread.read_grid._readgridBF
def _readgridBF
Definition: read_grid.py:31
flextest.flexread.read_grid.readgridV8
def readgridV8
Definition: read_grid.py:247
flextest.flexread.read_grid.Structure
Definition: read_grid.py:925

Here is the call graph for this function:

Here is the caller graph for this function:

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