Paparazzi UAS  v5.15_devel-81-gd13dafb
Paparazzi is a free software Unmanned Aircraft System.
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
mesonh_atmosphere.py
Go to the documentation of this file.
1 """
2  Reading mesoNH netCDF files and providing an interface to read chunks
3 of data and interpolate points.
4 
5  TODO : test with large netCDF files, add mapping of variables as an argument.
6 """
7 
8 from __future__ import division, print_function, absolute_import
9 from netCDF4 import MFDataset, Dataset
10 import numpy as np
11 from scipy.interpolate import RegularGridInterpolator
12 
13 
14 def find_lt(a, x):
15  'Find rightmost value less than x'
16  return np.searchsorted(a, x, side='left')-1
17 
18 
19 def find_le(a, x):
20  'Find rightmost value less than or equal to x'
21  return np.searchsorted(a, x, side='right')-1
22 
23 
24 def find_gt(a, x):
25  'Find leftmost value greater than x'
26  return np.searchsorted(a, x, side='right')
27 
28 
29 def find_ge(a, x):
30  'Find leftmost item greater than or equal to x'
31  return np.searchsorted(a, x, side='left')
32 
33 WX = 'UT' # Wind x,y and z composants
34 WY = 'VT'
35 WZ = 'WT'
36 PTEMP = 'THT' # potential temperature
37 RVAP = 'RVT' # vapor mixing ratio
38 LWATER = 'RCT' # liquid water
39 ZSCALE = 'VLEV' # vertical level (height) scale
40 XSCALE = 'S_N_direction' # south -> north axis scale
41 YSCALE = 'W_E_direction' # west -> east axis scale
42 TIME = 'time' # time dimension
43 X = 2
44 Y = 3
45 Z = 1
46 T = 0
47 
48 
50  """
51  Interpolation in space-time from mesoNH data over a grid.
52 
53  The functions valued on the grid are assumed to be periodic in the
54  x and y axis (south->north and west->east) :
55  f(x)=f(x+xmax-xmin) where [xmin,xmax] are the bounds on x
56  of the mesoNH grid.
57 
58  Clipping is done on the z (height) and time axis :
59  f(z)=f(min(max(z,zmin),zmax)) where [zmin,zmax] are the bound on z
60  of the mesoNH grid.
61 
62  Two interpolation methods are currently supported : nearest neighbour
63  and linear.
64 
65  Uses RegularGridInterpolator from scipy package.
66 
67  """
68  data = []
69 
70  interpolator = None
71  boundsmax = None
72  boundsmin = None
73 
74  def __init__(self, files, tstep, tinit=0):
75  if np.isscalar(files):
76  self.data = Dataset(files)
77  else:
78  self.data = MFDataset(files)
79  pz = self.data.variables[ZSCALE][:, 0, 0]
80  px = self.data.variables[XSCALE]
81  py = self.data.variables[YSCALE]
82  pt = np.array([tinit + i * tstep for i in
83  range(self.data.variables[TIME].shape[0])])
84  self.boundsmax = [pt[-1], pz[-1], px[-1], py[-1]]
85  self.boundsmin = [pt[0], pz[0], px[0], py[0]]
86 
87  self.grid_coordinates = (pt, pz, px, py)
88  self.grid_shape = (len(pt), len(pz), len(px), len(py))
89 
90  def get_points(self, points, var, method='nearest'):
91  """ Get value of variable on points
92 
93  Arguments:
94  points: a ndarray containing the point coordinates on the last
95  dimension
96  var: the name of the variable in the mesoNH file(s)
97  method: 'nearest' and 'linear' interpolation are currently supported
98  """
99  points = np.array(points)
100  p = self._apply_bounds(points)
101  caxes = tuple(range(p.ndim-1))
102  bounds = zip(p.min(axis=caxes), p.max(axis=caxes))
103  interpolator = self._get_interpolator(bounds, var, method)
104  return interpolator(p, method).squeeze()
105 
106  def _get_interpolator(self, bounds, var, method="nearest"):
107  slice_indexes, coordinates = self._slicyfy(bounds)
108  values = self._get_var_values(var, slice_indexes)
109  ip = RegularGridInterpolator(coordinates, values, method)
110  return ip
111 
112  def _slicyfy(self, bounds):
113  slice_indexes = ()
114  coordinates = ()
115  for d, b in enumerate(bounds):
116  dslice = slice(find_le(self.grid_coordinates[d], b[0]),
117  find_gt(self.grid_coordinates[d], b[1])+1)
118  slice_indexes += dslice,
119  coordinates += self.grid_coordinates[d][dslice],
120 
121  return slice_indexes, coordinates
122 
123  def _get_var_values(self, var, idx=Ellipsis):
124  return self.data[var][idx]
125 
126  def _apply_bounds(self, point):
127  x = point[..., X]
128  y = point[..., Y]
129  z = point[..., Z]
130  t = point[..., T]
131 
132  p = np.ndarray(point.shape)
133  p[..., T] = np.clip(t, self.boundsmin[T], self.boundsmax[T])
134  p[..., Z] = np.clip(z, self.boundsmin[Z], self.boundsmax[Z])
135  p[..., X] = (x-self.boundsmin[X]) % (self.boundsmax[X]-self.boundsmin[X])\
136  + self.boundsmin[X]
137  p[..., Y] = (y-self.boundsmin[Y]) % (self.boundsmax[Y]-self.boundsmin[Y])\
138  + self.boundsmin[Y]
139 
140  return p
141 
142  def get_wind(self, points, method='nearest'):
143  """Convenience method for getting 3D wind. See get_points."""
144  return np.array((self.get_points(points, WX, method),
145  self.get_points(points, WY, method),
146  self.get_points(points, WZ, method)
147  ))