FLEXPART Testing Environment CTBTO WO8
 All Classes Namespaces Files Functions Variables Pages
FlexpartOutput.py
Go to the documentation of this file.
1 # -*- coding: utf-8 -*-
2 """
3 FlexpartOutput.py
4 """
5 
6 import logging
7 import numpy as np
8 import os
9 import sys
10 
11 #import pflexible as pf
12 import read_header as rh
13 import read_grid as rg
14 
15 class FlexpartOutput(object):
16 
17  """High level interface to encapsulate and provide access to
18  the output of a specified FLEXPART simulation"""
19 
20  def __init__(self, output_dir=None, nest=False):
21 
22  """Constructor
23  output_dir: full path to directory of FLEXPART output
24  nest: Denotes whether this object is for the mother grid or a nest
25  """
26 
27 
28 
29  self._output_dir = output_dir
30 
31  if self._output_dir and os.path.isdir(self._output_dir):
32  try:
33  self._Header = rh.Header( self._output_dir, nested=nest )
34  except:
35  raise IOError('Problem getting output data header')
36  else:
37  raise IOError('Bad output data dir: ' + str(self._output_dir))
38 
39  #print self._Header
40  def get_horiz_slice(self, timestamp=None, level=1, species=1,
41  release=1, age_class=1):
42 
43  """Extracts 2D horizontal slice. Note that level, species and
44  release index values for the user start at 1, but they are all
45  decremented internal to this method so they start at 0.
46  Returned 2D grid is a NumPy array
47 
48  Note that if timestamp is not defined, then a timestamp "close to
49  the middle" will be selected.
50 
51  If there is any kind of error, we return None
52  """
53 
54  # if timestamp is not defined, choose a "mid way" timestamp for default
55  if not timestamp:
56  ntimes = len(self._Header['available_dates'])
57  tstamp_idx = ntimes/2
58  timestamp = self._Header['available_dates'][tstamp_idx]
59 
60  # Get the full grid for specified timestamp and species
61  so_far_so_good = True
62  try:
63 
64  G = rg.read_grid(self._Header, date=timestamp,
65  nspec_ret=species-1,
66  pspec_ret=release-1,
67  age_ret=age_class-1,
68  )
69 
70  except:
71  logging.error('Problem getting G from read_grid()')
72  so_far_so_good = False
73 
74  if so_far_so_good:
75  try:
76  # Get the grid (x,y,z,release)
77  the_grid = G[ (species-1, timestamp) ]
78  except:
79  logging.error('Problem getting the_grid from G')
80  so_far_so_good = False
81 
82  if so_far_so_good:
83  try:
84  # Get the 2D slice horizontal slice (x,y)
85  the_slice = the_grid.grid[:, :, level-1, release-1]
86 
87  except:
88  logging.error('Problem getting the_slice from the_grid')
89  so_far_so_good = False
90 
91  if so_far_so_good:
92  return the_slice
93  else:
94  return None
95 
96 
97  def get_deposition(self, timestamp=None, species=1,
98  release=1, age_class=1,
99  depo_type='dry'):
100 
101  """Extracts 2D horizontal slice of deposition.
102  Note that species and
103  release index values for the user start at 1, but they are all
104  decremented internal to this method so they start at 0.
105  Returned 2D grid is a NumPy array
106 
107  Note that if timestamp is not defined, then a timestamp "close to
108  the middle" will be selected.
109 
110  If there is any kind of error, we return None
111  """
112 
113  # if timestamp is not defined, choose a "mid way" timestamp for default
114  if not timestamp:
115  ntimes = len(self._Header['available_dates'])
116  tstamp_idx = ntimes/2
117  timestamp = self._Header['available_dates'][tstamp_idx]
118 
119  # Get the full grid for specified timestamp and species
120  so_far_so_good = True
121  try:
122  dry_flag=False; wet_flag=False
123  if depo_type=='dry':
124  dry_flag = True
125  if depo_type=='wet':
126  wet_flag = True
127 
128  G = rg.read_grid(self._Header, date=timestamp,
129  nspec_ret=species-1,
130  pspec_ret=release-1,
131  age_ret=age_class-1,
132  getdry=dry_flag,
133  getwet=wet_flag
134  )
135 
136  #print "G"
137  #print G
138  #print G.shape
139  #sys.exit()
140 
141  except:
142  logging.error('Problem getting G from read_grid()')
143  so_far_so_good = False
144 
145  if so_far_so_good:
146  try:
147  # Get the grid (x,y,z,release)
148  the_grid = G[ (species-1, timestamp) ]
149  except:
150  logging.error('Problem getting the_grid from G')
151  so_far_so_good = False
152 
153  if so_far_so_good:
154  try:
155  # Get the 2D slice horizontal slice (x,y)
156 
157  # 2016-03-23 -- DJM modification - apparently John put a
158  # "squeeze()" call on these arrays when they came out in
159  # read_grid(), so they can come out as either 2D or 3D. I
160  # need consistency, so if they're 2D, I'm expanding back out
161  # to 3D
162  #
163  # WARNING - DJM - 2016-03-23 - I'm not sure why I used
164  # "release" as an index in one of the following rather than
165  # species. I "think" it was logical, and I just went by
166  # the dimensions I had available in a test case, but I'm
167  # not 100% positive that this ain't going to get us into
168  # trouble one day.
169  if depo_type=='dry':
170  #print the_grid.dry
171  #print 'dry shape: ', the_grid.dry.shape
172  if the_grid.dry.ndim == 2:
173  dry_grid = np.expand_dims(the_grid.dry, axis=2)
174  else:
175  dry_grid = the_grid.dry
176  the_depo = dry_grid[:, :, release-1]
177  elif depo_type=='wet':
178  #print 'wet shape: ', the_grid.wet.shape
179  if the_grid.wet.ndim == 2:
180  wet_grid = np.expand_dims(the_grid.wet, axis=2)
181  else:
182  wet_grid = the_grid.wet
183  the_depo = wet_grid[:, :, release-1]
184  else:
185  the_depo = None
186  except:
187  logging.error('Problem getting the_depo from the_grid')
188  so_far_so_good = False
189 
190  if so_far_so_good:
191  return the_depo
192  else:
193  return None
194 
195 
196 
197  def get_volume(self, timestamp=None, level_list=None, species=1,
198  release=1, age_class=1):
199 
200  """Extracts 3D volume. if level_list is not defined, we return the full
201  volume for the timestamp. If level_list is defined,then we extract the
202  specified levels and create a volume from them
203 
204  Note that level, species and
205  release index values for the user start at 1, but they are all
206  decremented internal to this method so they start at 0.
207  Returned 2D grid is a NumPy array
208 
209  If there is any kind of error, we return None
210  """
211 
212  # if timestamp is not defined, choose a "mid way" timestamp for default
213  if not timestamp:
214  ntimes = len(self._Header['available_dates'])
215  tstamp_idx = ntimes/2
216  timestamp = self._Header['available_dates'][tstamp_idx]
217 
218  # Get the full grid for specified timestamp and species
219  so_far_so_good = True
220  try:
221  G = rg.read_grid(self._Header, date=timestamp,
222  nspec_ret=species-1,
223  pspec_ret=release-1,
224  age_ret=age_class-1,
225  )
226 
227  except:
228  logging.error('Problem getting G from read_grid()')
229  so_far_so_good = False
230 
231  if so_far_so_good:
232  try:
233  # Get the grid (x,y,z,release)
234  the_grid = G[ (species-1, timestamp) ]
235  except:
236  logging.error('Problem getting the_grid from G')
237  so_far_so_good = False
238 
239  if so_far_so_good:
240  if level_list:
241  # Decrement indices to 0 start
242  level_indices = list(map(lambda x: x-1, level_list))
243  try:
244  # Get the volume (x,y,z)
245  the_volume = the_grid.grid[:, :, level_indices, release-1]
246  except:
247  logging.error('Problem getting the_volume from the_grid')
248  so_far_so_good = False
249  else:
250  # Extract the full volume
251  # Make sure all level numbers are less than levels dimension
252  # of grid
253 
254  try:
255  # Get the volume (x,y,z)
256  the_volume = the_grid.grid[:, :, :, release-1]
257  except:
258  logging.error('Problem getting the_volume from the_grid')
259  so_far_so_good = False
260 
261  if so_far_so_good:
262  return the_volume
263  else:
264  return None
265 
266 
267  def get_horiz_timeseries(self, timestamp_list=None,
268  level=1, species=1,
269  release=1, age_class=1):
270 
271  """Extracts timeseries of 2D horizontal slices. If timestamp_list is
272  defined correctly, it extracts just the timestamps (in YYYYMMDDHHMMss
273  format) listed. Otherwise, it extracts all timestamps. The timeseries
274  is returned in a 3D NumPy array, indexed by (t, x, y)"""
275 
276  # Get the dimensions to create a 3D (t, x, y) array
277  if not timestamp_list:
278  timestamp_list = self._Header['available_dates']
279  Nx = self._Header.nx
280  Ny = self._Header.ny
281 
282  slice_timeseries = np.ndarray( (len(timestamp_list), Nx, Ny),
283  dtype=np.float64 )
284 
285 
286  #print Nx, Ny
287 
288  time_idx = 0
289  for ts in timestamp_list:
290  next_slice = self.get_horiz_slice(timestamp=ts,
291  level=level,
292  species=species,
293  release=release,
294  age_class=age_class)
295 
296 
297 
298  #print next_slice
299  slice_timeseries[time_idx, :, :] = next_slice
300  time_idx += 1
301 
302 
303  return slice_timeseries
304 
305 
306  def get_deposition_timeseries(self, timestamp_list=None,
307  species=1,
308  release=1, age_class=1,
309  depo_type='dry'):
310 
311  """Extracts timeseries of 2D deposition. If timestamp_list is
312  defined correctly, it extracts just the timestamps (in YYYYMMDDHHMMss
313  format) listed. Otherwise, it extracts all timestamps. The timeseries
314  is returned in a 3D NumPy array, indexed by (t, x, y)"""
315 
316  # Get the dimensions to create a 3D (t, x, y) array
317  if not timestamp_list:
318  timestamp_list = self._Header['available_dates']
319  Nx = self._Header.nx
320  Ny = self._Header.ny
321 
322  depo_timeseries = np.ndarray( (len(timestamp_list), Nx, Ny),
323  dtype=np.float64 )
324 
325 
326  #print Nx, Ny
327 
328  time_idx = 0
329  for ts in timestamp_list:
330  next_depo = self.get_deposition(timestamp=ts,
331  species=species,
332  release=release,
333  age_class=age_class,
334  depo_type=depo_type
335  )
336 
337 
338  #print 'next_depo: ', next_depo
339  #print next_depo.shape
340 
341  #print next_slice
342  depo_timeseries[time_idx, :, :] = next_depo
343  time_idx += 1
344 
345 
346  return depo_timeseries
347 
348 
349  def get_volume_timeseries(self, timestamp_list=None,
350  level_list=None, species=1,
351  release=1, age_class=1):
352 
353  """Extracts timeseries of 2D horizontal slices. If timestamp_list is
354  defined correctly, it extracts just the timestamps (in YYYYMMDDHHMMss
355  format) listed. Otherwise, it extracts all timestamps. The timeseries
356  is returned in a 3D NumPy array, indexed by (t, x, y)"""
357 
358  # Get the dimensions to create a 3D (t, x, y) array
359  if not timestamp_list:
360  timestamp_list = self._Header['available_dates']
361  Nx = self._Header.nx
362  Ny = self._Header.ny
363 
364  if not level_list:
365  Nz = self._Header.nz
366  else:
367  Nz = len( level_list )
368 
369  volume_timeseries = np.ndarray( (len(timestamp_list), Nx, Ny, Nz),
370  dtype=np.float64 )
371 
372 
373  #print Nx, Ny, Nz
374 
375  time_idx = 0
376  for ts in timestamp_list:
377  next_volume = self.get_volume(timestamp=ts,
378  level_list=level_list,
379  species=species,
380  release=release,
381  age_class=age_class)
382 
383 
384 
385  #print next_volume
386  volume_timeseries[time_idx, :, :] = next_volume
387  time_idx += 1
388 
389 
390  return volume_timeseries
391 
392 
393 
394  def get_timestamp_list(self):
395 
396  """
397  Returns the list of timestamps as Python list
398  """
399  timestamp_list = self._Header['available_dates']
400 
401  return timestamp_list
402 
403 
404 
405 
406 if __name__=="_main__":
407  pass
Generated on Sun Mar 27 2016 01:16:01 for FLEXPART Testing Environment CTBTO WO8 by   doxygen 1.8.6