|
1
|
|
|
# Copyright (c) 2016,2017 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
"""Contains a collection of generally useful calculation tools.""" |
|
5
|
|
|
|
|
6
|
|
|
import functools |
|
7
|
|
|
import warnings |
|
8
|
|
|
|
|
9
|
|
|
import numpy as np |
|
10
|
|
|
import numpy.ma as ma |
|
11
|
|
|
from scipy.spatial import cKDTree |
|
12
|
|
|
|
|
13
|
|
|
from . import height_to_pressure_std, pressure_to_height_std |
|
14
|
|
|
from ..package_tools import Exporter |
|
15
|
|
|
from ..units import check_units, units |
|
16
|
|
|
|
|
17
|
|
|
exporter = Exporter(globals()) |
|
18
|
|
|
|
|
19
|
|
|
|
|
20
|
|
|
@exporter.export |
|
21
|
|
|
def resample_nn_1d(a, centers): |
|
22
|
|
|
"""Return one-dimensional nearest-neighbor indexes based on user-specified centers. |
|
23
|
|
|
|
|
24
|
|
|
Parameters |
|
25
|
|
|
---------- |
|
26
|
|
|
a : array-like |
|
27
|
|
|
1-dimensional array of numeric values from which to |
|
28
|
|
|
extract indexes of nearest-neighbors |
|
29
|
|
|
centers : array-like |
|
30
|
|
|
1-dimensional array of numeric values representing a subset of values to approximate |
|
31
|
|
|
|
|
32
|
|
|
Returns |
|
33
|
|
|
------- |
|
34
|
|
|
An array of indexes representing values closest to given array values |
|
35
|
|
|
|
|
36
|
|
|
""" |
|
37
|
|
|
ix = [] |
|
38
|
|
|
for center in centers: |
|
39
|
|
|
index = (np.abs(a - center)).argmin() |
|
40
|
|
|
if index not in ix: |
|
41
|
|
|
ix.append(index) |
|
42
|
|
|
return ix |
|
43
|
|
|
|
|
44
|
|
|
|
|
45
|
|
|
@exporter.export |
|
46
|
|
|
def nearest_intersection_idx(a, b): |
|
47
|
|
|
"""Determine the index of the point just before two lines with common x values. |
|
48
|
|
|
|
|
49
|
|
|
Parameters |
|
50
|
|
|
---------- |
|
51
|
|
|
a : array-like |
|
52
|
|
|
1-dimensional array of y-values for line 1 |
|
53
|
|
|
b : array-like |
|
54
|
|
|
1-dimensional array of y-values for line 2 |
|
55
|
|
|
|
|
56
|
|
|
Returns |
|
57
|
|
|
------- |
|
58
|
|
|
An array of indexes representing the index of the values |
|
59
|
|
|
just before the intersection(s) of the two lines. |
|
60
|
|
|
|
|
61
|
|
|
""" |
|
62
|
|
|
# Difference in the two y-value sets |
|
63
|
|
|
difference = a - b |
|
64
|
|
|
|
|
65
|
|
|
# Determine the point just before the intersection of the lines |
|
66
|
|
|
# Will return multiple points for multiple intersections |
|
67
|
|
|
sign_change_idx, = np.nonzero(np.diff(np.sign(difference))) |
|
68
|
|
|
|
|
69
|
|
|
return sign_change_idx |
|
70
|
|
|
|
|
71
|
|
|
|
|
72
|
|
|
@exporter.export |
|
73
|
|
|
@units.wraps(('=A', '=B'), ('=A', '=B', '=B')) |
|
74
|
|
|
def find_intersections(x, a, b, direction='all'): |
|
75
|
|
|
"""Calculate the best estimate of intersection. |
|
76
|
|
|
|
|
77
|
|
|
Calculates the best estimates of the intersection of two y-value |
|
78
|
|
|
data sets that share a common x-value set. |
|
79
|
|
|
|
|
80
|
|
|
Parameters |
|
81
|
|
|
---------- |
|
82
|
|
|
x : array-like |
|
83
|
|
|
1-dimensional array of numeric x-values |
|
84
|
|
|
a : array-like |
|
85
|
|
|
1-dimensional array of y-values for line 1 |
|
86
|
|
|
b : array-like |
|
87
|
|
|
1-dimensional array of y-values for line 2 |
|
88
|
|
|
direction : string, optional |
|
89
|
|
|
specifies direction of crossing. 'all', 'increasing' (a becoming greater than b), |
|
90
|
|
|
or 'decreasing' (b becoming greater than a). Defaults to 'all'. |
|
91
|
|
|
|
|
92
|
|
|
Returns |
|
93
|
|
|
------- |
|
94
|
|
|
A tuple (x, y) of array-like with the x and y coordinates of the |
|
95
|
|
|
intersections of the lines. |
|
96
|
|
|
|
|
97
|
|
|
""" |
|
98
|
|
|
# Find the index of the points just before the intersection(s) |
|
99
|
|
|
nearest_idx = nearest_intersection_idx(a, b) |
|
100
|
|
|
next_idx = nearest_idx + 1 |
|
101
|
|
|
|
|
102
|
|
|
# Determine the sign of the change |
|
103
|
|
|
sign_change = np.sign(a[next_idx] - b[next_idx]) |
|
104
|
|
|
|
|
105
|
|
|
# x-values around each intersection |
|
106
|
|
|
_, x0 = _next_non_masked_element(x, nearest_idx) |
|
107
|
|
|
_, x1 = _next_non_masked_element(x, next_idx) |
|
108
|
|
|
|
|
109
|
|
|
# y-values around each intersection for the first line |
|
110
|
|
|
_, a0 = _next_non_masked_element(a, nearest_idx) |
|
111
|
|
|
_, a1 = _next_non_masked_element(a, next_idx) |
|
112
|
|
|
|
|
113
|
|
|
# y-values around each intersection for the second line |
|
114
|
|
|
_, b0 = _next_non_masked_element(b, nearest_idx) |
|
115
|
|
|
_, b1 = _next_non_masked_element(b, next_idx) |
|
116
|
|
|
|
|
117
|
|
|
# Calculate the x-intersection. This comes from finding the equations of the two lines, |
|
118
|
|
|
# one through (x0, a0) and (x1, a1) and the other through (x0, b0) and (x1, b1), |
|
119
|
|
|
# finding their intersection, and reducing with a bunch of algebra. |
|
120
|
|
|
delta_y0 = a0 - b0 |
|
121
|
|
|
delta_y1 = a1 - b1 |
|
122
|
|
|
intersect_x = (delta_y1 * x0 - delta_y0 * x1) / (delta_y1 - delta_y0) |
|
123
|
|
|
|
|
124
|
|
|
# Calculate the y-intersection of the lines. Just plug the x above into the equation |
|
125
|
|
|
# for the line through the a points. One could solve for y like x above, but this |
|
126
|
|
|
# causes weirder unit behavior and seems a little less good numerically. |
|
127
|
|
|
intersect_y = ((intersect_x - x0) / (x1 - x0)) * (a1 - a0) + a0 |
|
128
|
|
|
|
|
129
|
|
|
# If there's no intersections, return |
|
130
|
|
|
if len(intersect_x) == 0: |
|
131
|
|
|
return intersect_x, intersect_y |
|
132
|
|
|
|
|
133
|
|
|
# Check for duplicates |
|
134
|
|
|
duplicate_mask = (np.ediff1d(intersect_x, to_end=1) != 0) |
|
135
|
|
|
|
|
136
|
|
|
# Make a mask based on the direction of sign change desired |
|
137
|
|
|
if direction == 'increasing': |
|
138
|
|
|
mask = sign_change > 0 |
|
139
|
|
|
elif direction == 'decreasing': |
|
140
|
|
|
mask = sign_change < 0 |
|
141
|
|
|
elif direction == 'all': |
|
142
|
|
|
return intersect_x[duplicate_mask], intersect_y[duplicate_mask] |
|
143
|
|
|
else: |
|
144
|
|
|
raise ValueError('Unknown option for direction: {0}'.format(str(direction))) |
|
145
|
|
|
|
|
146
|
|
|
return intersect_x[mask & duplicate_mask], intersect_y[mask & duplicate_mask] |
|
147
|
|
|
|
|
148
|
|
|
|
|
149
|
|
|
@exporter.export |
|
150
|
|
|
def interpolate_nans(x, y, kind='linear'): |
|
151
|
|
|
"""Interpolate NaN values in y. |
|
152
|
|
|
|
|
153
|
|
|
Interpolate NaN values in the y dimension. Works with unsorted x values. |
|
154
|
|
|
|
|
155
|
|
|
Parameters |
|
156
|
|
|
---------- |
|
157
|
|
|
x : array-like |
|
158
|
|
|
1-dimensional array of numeric x-values |
|
159
|
|
|
y : array-like |
|
160
|
|
|
1-dimensional array of numeric y-values |
|
161
|
|
|
kind : string |
|
162
|
|
|
specifies the kind of interpolation x coordinate - 'linear' or 'log', optional. |
|
163
|
|
|
Defaults to 'linear'. |
|
164
|
|
|
|
|
165
|
|
|
Returns |
|
166
|
|
|
------- |
|
167
|
|
|
An array of the y coordinate data with NaN values interpolated. |
|
168
|
|
|
|
|
169
|
|
|
""" |
|
170
|
|
|
x_sort_args = np.argsort(x) |
|
171
|
|
|
x = x[x_sort_args] |
|
172
|
|
|
y = y[x_sort_args] |
|
173
|
|
|
nans = np.isnan(y) |
|
174
|
|
|
if kind == 'linear': |
|
175
|
|
|
y[nans] = np.interp(x[nans], x[~nans], y[~nans]) |
|
176
|
|
|
elif kind == 'log': |
|
177
|
|
|
y[nans] = np.interp(np.log(x[nans]), np.log(x[~nans]), y[~nans]) |
|
178
|
|
|
else: |
|
179
|
|
|
raise ValueError('Unknown option for kind: {0}'.format(str(kind))) |
|
180
|
|
|
return y[x_sort_args] |
|
181
|
|
|
|
|
182
|
|
|
|
|
183
|
|
|
def _next_non_masked_element(a, idx): |
|
184
|
|
|
"""Return the next non masked element of a masked array. |
|
185
|
|
|
|
|
186
|
|
|
If an array is masked, return the next non-masked element (if the given index is masked). |
|
187
|
|
|
If no other unmasked points are after the given masked point, returns none. |
|
188
|
|
|
|
|
189
|
|
|
Parameters |
|
190
|
|
|
---------- |
|
191
|
|
|
a : array-like |
|
192
|
|
|
1-dimensional array of numeric values |
|
193
|
|
|
idx : integer |
|
194
|
|
|
index of requested element |
|
195
|
|
|
|
|
196
|
|
|
Returns |
|
197
|
|
|
------- |
|
198
|
|
|
Index of next non-masked element and next non-masked element |
|
199
|
|
|
|
|
200
|
|
|
""" |
|
201
|
|
|
try: |
|
202
|
|
|
next_idx = idx + a[idx:].mask.argmin() |
|
203
|
|
|
if ma.is_masked(a[next_idx]): |
|
204
|
|
|
return None, None |
|
205
|
|
|
else: |
|
206
|
|
|
return next_idx, a[next_idx] |
|
207
|
|
|
except (AttributeError, TypeError, IndexError): |
|
208
|
|
|
return idx, a[idx] |
|
209
|
|
|
|
|
210
|
|
|
|
|
211
|
|
|
def delete_masked_points(*arrs): |
|
212
|
|
|
"""Delete masked points from arrays. |
|
213
|
|
|
|
|
214
|
|
|
Takes arrays and removes masked points to help with calculations and plotting. |
|
215
|
|
|
|
|
216
|
|
|
Parameters |
|
217
|
|
|
---------- |
|
218
|
|
|
arrs : one or more array-like |
|
219
|
|
|
source arrays |
|
220
|
|
|
|
|
221
|
|
|
Returns |
|
222
|
|
|
------- |
|
223
|
|
|
arrs : one or more array-like |
|
224
|
|
|
arrays with masked elements removed |
|
225
|
|
|
|
|
226
|
|
|
""" |
|
227
|
|
|
if any(hasattr(a, 'mask') for a in arrs): |
|
228
|
|
|
keep = ~functools.reduce(np.logical_or, (np.ma.getmaskarray(a) for a in arrs)) |
|
229
|
|
|
return tuple(ma.asarray(a[keep]) for a in arrs) |
|
230
|
|
|
else: |
|
231
|
|
|
return arrs |
|
232
|
|
|
|
|
233
|
|
|
|
|
234
|
|
|
@exporter.export |
|
235
|
|
|
def reduce_point_density(points, radius, priority=None): |
|
236
|
|
|
r"""Return a mask to reduce the density of points in irregularly-spaced data. |
|
237
|
|
|
|
|
238
|
|
|
This function is used to down-sample a collection of scattered points (e.g. surface |
|
239
|
|
|
data), returning a mask that can be used to select the points from one or more arrays |
|
240
|
|
|
(e.g. arrays of temperature and dew point). The points selected can be controlled by |
|
241
|
|
|
providing an array of ``priority`` values (e.g. rainfall totals to ensure that |
|
242
|
|
|
stations with higher precipitation remain in the mask). |
|
243
|
|
|
|
|
244
|
|
|
Parameters |
|
245
|
|
|
---------- |
|
246
|
|
|
points : (N, K) array-like |
|
247
|
|
|
N locations of the points in K dimensional space |
|
248
|
|
|
radius : float |
|
249
|
|
|
minimum radius allowed between points |
|
250
|
|
|
priority : (N, K) array-like, optional |
|
251
|
|
|
If given, this should have the same shape as ``points``; these values will |
|
252
|
|
|
be used to control selection priority for points. |
|
253
|
|
|
|
|
254
|
|
|
Returns |
|
255
|
|
|
------- |
|
256
|
|
|
(N,) array-like of boolean values indicating whether points should be kept. This |
|
257
|
|
|
can be used directly to index numpy arrays to return only the desired points. |
|
258
|
|
|
|
|
259
|
|
|
Examples |
|
260
|
|
|
-------- |
|
261
|
|
|
>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.) |
|
262
|
|
|
array([ True, False, True], dtype=bool) |
|
263
|
|
|
>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1., |
|
264
|
|
|
... priority=np.array([0.1, 0.9, 0.3])) |
|
265
|
|
|
array([False, True, False], dtype=bool) |
|
266
|
|
|
|
|
267
|
|
|
""" |
|
268
|
|
|
# Handle 1D input |
|
269
|
|
|
if points.ndim < 2: |
|
270
|
|
|
points = points.reshape(-1, 1) |
|
271
|
|
|
|
|
272
|
|
|
# Make a kd-tree to speed searching of data. |
|
273
|
|
|
tree = cKDTree(points) |
|
274
|
|
|
|
|
275
|
|
|
# Need to use sorted indices rather than sorting the position |
|
276
|
|
|
# so that the keep mask matches *original* order. |
|
277
|
|
|
if priority is not None: |
|
278
|
|
|
# Need to sort the locations in decreasing priority. |
|
279
|
|
|
sorted_indices = np.argsort(priority)[::-1] |
|
280
|
|
|
else: |
|
281
|
|
|
# Take advantage of iterator nature of range here to avoid making big lists |
|
282
|
|
|
sorted_indices = range(len(points)) |
|
283
|
|
|
|
|
284
|
|
|
# Keep all points initially |
|
285
|
|
|
keep = np.ones(len(points), dtype=np.bool) |
|
286
|
|
|
|
|
287
|
|
|
# Loop over all the potential points |
|
288
|
|
|
for ind in sorted_indices: |
|
289
|
|
|
# Only proceed if we haven't already excluded this point |
|
290
|
|
|
if keep[ind]: |
|
291
|
|
|
# Find the neighbors and eliminate them |
|
292
|
|
|
neighbors = tree.query_ball_point(points[ind], radius) |
|
293
|
|
|
keep[neighbors] = False |
|
294
|
|
|
|
|
295
|
|
|
# We just removed ourselves, so undo that |
|
296
|
|
|
keep[ind] = True |
|
297
|
|
|
|
|
298
|
|
|
return keep |
|
299
|
|
|
|
|
300
|
|
|
|
|
301
|
|
|
def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True): |
|
302
|
|
|
"""Calculate the bounding pressure and height in a layer. |
|
303
|
|
|
|
|
304
|
|
|
Given pressure, optional heights, and a bound, return either the closest pressure/height |
|
305
|
|
|
or interpolated pressure/height. If no heights are provided, a standard atmosphere is |
|
306
|
|
|
assumed. |
|
307
|
|
|
|
|
308
|
|
|
Parameters |
|
309
|
|
|
---------- |
|
310
|
|
|
pressure : `pint.Quantity` |
|
311
|
|
|
Atmospheric pressures |
|
312
|
|
|
bound : `pint.Quantity` |
|
313
|
|
|
Bound to retrieve (in pressure or height) |
|
314
|
|
|
heights : `pint.Quantity`, optional |
|
315
|
|
|
Atmospheric heights associated with the pressure levels. Defaults to using |
|
316
|
|
|
heights calculated from ``pressure`` assuming a standard atmosphere. |
|
317
|
|
|
interpolate : boolean, optional |
|
318
|
|
|
Interpolate the bound or return the nearest. Defaults to True. |
|
319
|
|
|
|
|
320
|
|
|
Returns |
|
321
|
|
|
------- |
|
322
|
|
|
`pint.Quantity` |
|
323
|
|
|
The bound pressure and height. |
|
324
|
|
|
|
|
325
|
|
|
""" |
|
326
|
|
|
# Bound is given in pressure |
|
327
|
|
|
if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}: |
|
328
|
|
|
# If the bound is in the pressure data, we know the pressure bound exactly |
|
329
|
|
|
if bound in pressure: |
|
330
|
|
|
bound_pressure = bound |
|
331
|
|
|
# If we have heights, we know the exact height value, otherwise return standard |
|
332
|
|
|
# atmosphere height for the pressure |
|
333
|
|
|
if heights is not None: |
|
334
|
|
|
bound_height = heights[pressure == bound_pressure] |
|
335
|
|
|
else: |
|
336
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
|
337
|
|
|
# If bound is not in the data, return the nearest or interpolated values |
|
338
|
|
|
else: |
|
339
|
|
|
if interpolate: |
|
340
|
|
|
bound_pressure = bound # Use the user specified bound |
|
341
|
|
|
if heights is not None: # Interpolate heights from the height data |
|
342
|
|
|
bound_height = log_interp(bound_pressure, pressure, heights) |
|
343
|
|
|
else: # If not heights given, use the standard atmosphere |
|
344
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
|
345
|
|
|
else: # No interpolation, find the closest values |
|
346
|
|
|
idx = (np.abs(pressure - bound)).argmin() |
|
347
|
|
|
bound_pressure = pressure[idx] |
|
348
|
|
|
if heights is not None: |
|
349
|
|
|
bound_height = heights[idx] |
|
350
|
|
|
else: |
|
351
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
|
352
|
|
|
|
|
353
|
|
|
# Bound is given in height |
|
354
|
|
|
elif bound.dimensionality == {'[length]': 1.0}: |
|
355
|
|
|
# If there is height data, see if we have the bound or need to interpolate/find nearest |
|
356
|
|
|
if heights is not None: |
|
357
|
|
|
if bound in heights: # Bound is in the height data |
|
358
|
|
|
bound_height = bound |
|
359
|
|
|
bound_pressure = pressure[heights == bound] |
|
360
|
|
|
else: # Bound is not in the data |
|
361
|
|
|
if interpolate: |
|
362
|
|
|
bound_height = bound |
|
363
|
|
|
|
|
364
|
|
|
# Need to cast back to the input type since interp (up to at least numpy |
|
365
|
|
|
# 1.13 always returns float64. This can cause upstream users problems, |
|
366
|
|
|
# resulting in something like np.append() to upcast. |
|
367
|
|
|
bound_pressure = np.interp(np.atleast_1d(bound), heights, |
|
368
|
|
|
pressure).astype(bound.dtype) * pressure.units |
|
369
|
|
|
else: |
|
370
|
|
|
idx = (np.abs(heights - bound)).argmin() |
|
371
|
|
|
bound_pressure = pressure[idx] |
|
372
|
|
|
bound_height = heights[idx] |
|
373
|
|
|
else: # Don't have heights, so assume a standard atmosphere |
|
374
|
|
|
bound_height = bound |
|
375
|
|
|
bound_pressure = height_to_pressure_std(bound) |
|
376
|
|
|
# If interpolation is on, this is all we need, if not, we need to go back and |
|
377
|
|
|
# find the pressure closest to this and refigure the bounds |
|
378
|
|
|
if not interpolate: |
|
379
|
|
|
idx = (np.abs(pressure - bound_pressure)).argmin() |
|
380
|
|
|
bound_pressure = pressure[idx] |
|
381
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
|
382
|
|
|
|
|
383
|
|
|
# Bound has invalid units |
|
384
|
|
|
else: |
|
385
|
|
|
raise ValueError('Bound must be specified in units of length or pressure.') |
|
386
|
|
|
|
|
387
|
|
|
# If the bound is out of the range of the data, we shouldn't extrapolate |
|
388
|
|
|
if (bound_pressure < np.min(pressure)) or (bound_pressure > np.max(pressure)): |
|
389
|
|
|
raise ValueError('Specified bound is outside pressure range.') |
|
390
|
|
|
if heights is not None: |
|
391
|
|
|
if (bound_height > np.max(heights)) or (bound_height < np.min(heights)): |
|
392
|
|
|
raise ValueError('Specified bound is outside height range.') |
|
393
|
|
|
|
|
394
|
|
|
return bound_pressure, bound_height |
|
395
|
|
|
|
|
396
|
|
|
|
|
397
|
|
|
@exporter.export |
|
398
|
|
|
@check_units('[length]') |
|
399
|
|
|
def get_layer_heights(heights, depth, *args, **kwargs): |
|
400
|
|
|
"""Return an atmospheric layer from upper air data with the requested bottom and depth. |
|
401
|
|
|
|
|
402
|
|
|
This function will subset an upper air dataset to contain only the specified layer using |
|
403
|
|
|
the heights only. |
|
404
|
|
|
|
|
405
|
|
|
Parameters |
|
406
|
|
|
---------- |
|
407
|
|
|
heights : array-like |
|
408
|
|
|
Atmospheric heights |
|
409
|
|
|
depth : `pint.Quantity` |
|
410
|
|
|
The thickness of the layer |
|
411
|
|
|
*args : array-like |
|
412
|
|
|
Atmospheric variable(s) measured at the given pressures |
|
413
|
|
|
bottom : `pint.Quantity`, optional |
|
414
|
|
|
The bottom of the layer |
|
415
|
|
|
interpolate : bool, optional |
|
416
|
|
|
Interpolate the top and bottom points if they are not in the given data. Defaults |
|
417
|
|
|
to True. |
|
418
|
|
|
with_agl : bool, optional |
|
419
|
|
|
Returns the heights as above ground level by subtracting the minimum height in the |
|
420
|
|
|
provided heights. Defaults to False. |
|
421
|
|
|
|
|
422
|
|
|
Returns |
|
423
|
|
|
------- |
|
424
|
|
|
`pint.Quantity, pint.Quantity` |
|
425
|
|
|
The height and data variables of the layer |
|
426
|
|
|
|
|
427
|
|
|
""" |
|
428
|
|
|
bottom = kwargs.pop('bottom', None) |
|
429
|
|
|
interpolate = kwargs.pop('interpolate', True) |
|
430
|
|
|
with_agl = kwargs.pop('with_agl', False) |
|
431
|
|
|
|
|
432
|
|
|
# Make sure pressure and datavars are the same length |
|
433
|
|
|
for datavar in args: |
|
434
|
|
|
if len(heights) != len(datavar): |
|
435
|
|
|
raise ValueError('Height and data variables must have the same length.') |
|
436
|
|
|
|
|
437
|
|
|
# If we want things in AGL, subtract the minimum height from all height values |
|
438
|
|
|
if with_agl: |
|
439
|
|
|
sfc_height = np.min(heights) |
|
440
|
|
|
heights -= sfc_height |
|
441
|
|
|
|
|
442
|
|
|
# If the bottom is not specified, make it the surface |
|
443
|
|
|
if bottom is None: |
|
444
|
|
|
bottom = heights[0] |
|
445
|
|
|
|
|
446
|
|
|
# Make heights and arguments base units |
|
447
|
|
|
heights = heights.to_base_units() |
|
448
|
|
|
bottom = bottom.to_base_units() |
|
449
|
|
|
|
|
450
|
|
|
# Calculate the top of the layer |
|
451
|
|
|
top = bottom + depth |
|
452
|
|
|
|
|
453
|
|
|
ret = [] # returned data variables in layer |
|
454
|
|
|
|
|
455
|
|
|
# Ensure heights are sorted in ascending order |
|
456
|
|
|
sort_inds = np.argsort(heights) |
|
457
|
|
|
heights = heights[sort_inds] |
|
458
|
|
|
|
|
459
|
|
|
# Mask based on top and bottom |
|
460
|
|
|
inds = (heights >= bottom) & (heights <= top) |
|
461
|
|
|
heights_interp = heights[inds] |
|
462
|
|
|
|
|
463
|
|
|
# Interpolate heights at bounds if necessary and sort |
|
464
|
|
|
if interpolate: |
|
465
|
|
|
# If we don't have the bottom or top requested, append them |
|
466
|
|
|
if top not in heights_interp: |
|
467
|
|
|
heights_interp = np.sort(np.append(heights_interp, top)) * heights.units |
|
468
|
|
|
if bottom not in heights_interp: |
|
469
|
|
|
heights_interp = np.sort(np.append(heights_interp, bottom)) * heights.units |
|
470
|
|
|
|
|
471
|
|
|
ret.append(heights_interp) |
|
472
|
|
|
|
|
473
|
|
|
for datavar in args: |
|
474
|
|
|
# Ensure that things are sorted in ascending order |
|
475
|
|
|
datavar = datavar[sort_inds] |
|
476
|
|
|
|
|
477
|
|
|
if interpolate: |
|
478
|
|
|
# Interpolate for the possibly missing bottom/top values |
|
479
|
|
|
datavar_interp = interp(heights_interp, heights, datavar) |
|
480
|
|
|
datavar = datavar_interp |
|
481
|
|
|
else: |
|
482
|
|
|
datavar = datavar[inds] |
|
483
|
|
|
|
|
484
|
|
|
ret.append(datavar) |
|
485
|
|
|
return ret |
|
486
|
|
|
|
|
487
|
|
|
|
|
488
|
|
|
@exporter.export |
|
489
|
|
|
@check_units('[pressure]') |
|
490
|
|
|
def get_layer(pressure, *args, **kwargs): |
|
491
|
|
|
r"""Return an atmospheric layer from upper air data with the requested bottom and depth. |
|
492
|
|
|
|
|
493
|
|
|
This function will subset an upper air dataset to contain only the specified layer. The |
|
494
|
|
|
bottom of the layer can be specified with a pressure or height above the surface |
|
495
|
|
|
pressure. The bottom defaults to the surface pressure. The depth of the layer can be |
|
496
|
|
|
specified in terms of pressure or height above the bottom of the layer. If the top and |
|
497
|
|
|
bottom of the layer are not in the data, they are interpolated by default. |
|
498
|
|
|
|
|
499
|
|
|
Parameters |
|
500
|
|
|
---------- |
|
501
|
|
|
pressure : array-like |
|
502
|
|
|
Atmospheric pressure profile |
|
503
|
|
|
*args : array-like |
|
504
|
|
|
Atmospheric variable(s) measured at the given pressures |
|
505
|
|
|
heights: array-like, optional |
|
506
|
|
|
Atmospheric heights corresponding to the given pressures. Defaults to using |
|
507
|
|
|
heights calculated from ``p`` assuming a standard atmosphere. |
|
508
|
|
|
bottom : `pint.Quantity`, optional |
|
509
|
|
|
The bottom of the layer as a pressure or height above the surface pressure. Defaults |
|
510
|
|
|
to the lowest pressure or height given. |
|
511
|
|
|
depth : `pint.Quantity`, optional |
|
512
|
|
|
The thickness of the layer as a pressure or height above the bottom of the layer. |
|
513
|
|
|
Defaults to 100 hPa. |
|
514
|
|
|
interpolate : bool, optional |
|
515
|
|
|
Interpolate the top and bottom points if they are not in the given data. Defaults |
|
516
|
|
|
to True. |
|
517
|
|
|
|
|
518
|
|
|
Returns |
|
519
|
|
|
------- |
|
520
|
|
|
`pint.Quantity, pint.Quantity` |
|
521
|
|
|
The pressure and data variables of the layer |
|
522
|
|
|
|
|
523
|
|
|
""" |
|
524
|
|
|
# Pop off keyword arguments |
|
525
|
|
|
heights = kwargs.pop('heights', None) |
|
526
|
|
|
bottom = kwargs.pop('bottom', None) |
|
527
|
|
|
depth = kwargs.pop('depth', 100 * units.hPa) |
|
528
|
|
|
interpolate = kwargs.pop('interpolate', True) |
|
529
|
|
|
|
|
530
|
|
|
# If we get the depth kwarg, but it's None, set it to the default as well |
|
531
|
|
|
if depth is None: |
|
532
|
|
|
depth = 100 * units.hPa |
|
533
|
|
|
|
|
534
|
|
|
# Make sure pressure and datavars are the same length |
|
535
|
|
|
for datavar in args: |
|
536
|
|
|
if len(pressure) != len(datavar): |
|
537
|
|
|
raise ValueError('Pressure and data variables must have the same length.') |
|
538
|
|
|
|
|
539
|
|
|
# If the bottom is not specified, make it the surface pressure |
|
540
|
|
|
if bottom is None: |
|
541
|
|
|
bottom = pressure[0] |
|
542
|
|
|
|
|
543
|
|
|
bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom, |
|
544
|
|
|
heights=heights, |
|
545
|
|
|
interpolate=interpolate) |
|
546
|
|
|
|
|
547
|
|
|
# Calculate the top if whatever units depth is in |
|
548
|
|
|
if depth.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}: |
|
549
|
|
|
top = bottom_pressure - depth |
|
550
|
|
|
elif depth.dimensionality == {'[length]': 1}: |
|
551
|
|
|
top = bottom_height + depth |
|
552
|
|
|
else: |
|
553
|
|
|
raise ValueError('Depth must be specified in units of length or pressure') |
|
554
|
|
|
|
|
555
|
|
|
top_pressure, _ = _get_bound_pressure_height(pressure, top, heights=heights, |
|
556
|
|
|
interpolate=interpolate) |
|
557
|
|
|
|
|
558
|
|
|
ret = [] # returned data variables in layer |
|
559
|
|
|
|
|
560
|
|
|
# Ensure pressures are sorted in ascending order |
|
561
|
|
|
sort_inds = np.argsort(pressure) |
|
562
|
|
|
pressure = pressure[sort_inds] |
|
563
|
|
|
|
|
564
|
|
|
# Mask based on top and bottom pressure |
|
565
|
|
|
inds = (pressure <= bottom_pressure) & (pressure >= top_pressure) |
|
566
|
|
|
p_interp = pressure[inds] |
|
567
|
|
|
|
|
568
|
|
|
# Interpolate pressures at bounds if necessary and sort |
|
569
|
|
|
if interpolate: |
|
570
|
|
|
# If we don't have the bottom or top requested, append them |
|
571
|
|
|
if top_pressure not in p_interp: |
|
572
|
|
|
p_interp = np.sort(np.append(p_interp, top_pressure)) * pressure.units |
|
573
|
|
|
if bottom_pressure not in p_interp: |
|
574
|
|
|
p_interp = np.sort(np.append(p_interp, bottom_pressure)) * pressure.units |
|
575
|
|
|
|
|
576
|
|
|
ret.append(p_interp[::-1]) |
|
577
|
|
|
|
|
578
|
|
|
for datavar in args: |
|
579
|
|
|
# Ensure that things are sorted in ascending order |
|
580
|
|
|
datavar = datavar[sort_inds] |
|
581
|
|
|
|
|
582
|
|
|
if interpolate: |
|
583
|
|
|
# Interpolate for the possibly missing bottom/top values |
|
584
|
|
|
datavar_interp = log_interp(p_interp, pressure, datavar) |
|
585
|
|
|
datavar = datavar_interp |
|
586
|
|
|
else: |
|
587
|
|
|
datavar = datavar[inds] |
|
588
|
|
|
|
|
589
|
|
|
ret.append(datavar[::-1]) |
|
590
|
|
|
return ret |
|
591
|
|
|
|
|
592
|
|
|
|
|
593
|
|
|
@exporter.export |
|
594
|
|
|
@units.wraps(None, ('=A', '=A')) |
|
595
|
|
|
def interp(x, xp, *args, **kwargs): |
|
596
|
|
|
r"""Interpolates data with any shape over a specified axis. |
|
597
|
|
|
|
|
598
|
|
|
Interpolation over a specified axis for arrays of any shape. |
|
599
|
|
|
|
|
600
|
|
|
Parameters |
|
601
|
|
|
---------- |
|
602
|
|
|
x : array-like |
|
603
|
|
|
1-D array of desired interpolated values. |
|
604
|
|
|
|
|
605
|
|
|
xp : array-like |
|
606
|
|
|
The x-coordinates of the data points. |
|
607
|
|
|
|
|
608
|
|
|
args : array-like |
|
609
|
|
|
The data to be interpolated. Can be multiple arguments, all must be the same shape as |
|
610
|
|
|
xp. |
|
611
|
|
|
|
|
612
|
|
|
axis : int, optional |
|
613
|
|
|
The axis to interpolate over. Defaults to 0. |
|
614
|
|
|
|
|
615
|
|
|
fill_value: float, optional |
|
616
|
|
|
Specify handling of interpolation points out of data bounds. If None, will return |
|
617
|
|
|
ValueError if points are out of bounds. Defaults to nan. |
|
618
|
|
|
|
|
619
|
|
|
Returns |
|
620
|
|
|
------- |
|
621
|
|
|
array-like |
|
622
|
|
|
Interpolated values for each point with coordinates sorted in ascending order. |
|
623
|
|
|
|
|
624
|
|
|
Examples |
|
625
|
|
|
-------- |
|
626
|
|
|
>>> x = np.array([1., 2., 3., 4.]) |
|
627
|
|
|
>>> y = np.array([1., 2., 3., 4.]) |
|
628
|
|
|
>>> x_interp = np.array([2.5, 3.5]) |
|
629
|
|
|
>>> metpy.calc.interp(x_interp, x, y) |
|
630
|
|
|
array([ 2.5, 3.5]) |
|
631
|
|
|
|
|
632
|
|
|
Notes |
|
633
|
|
|
----- |
|
634
|
|
|
xp and args must be the same shape. |
|
635
|
|
|
|
|
636
|
|
|
""" |
|
637
|
|
|
# Pull out keyword args |
|
638
|
|
|
fill_value = kwargs.pop('fill_value', np.nan) |
|
639
|
|
|
axis = kwargs.pop('axis', 0) |
|
640
|
|
|
|
|
641
|
|
|
# Make x an array |
|
642
|
|
|
x = np.asanyarray(x).reshape(-1) |
|
643
|
|
|
|
|
644
|
|
|
# Save number of dimensions in xp |
|
645
|
|
|
ndim = xp.ndim |
|
646
|
|
|
|
|
647
|
|
|
# Sort input data |
|
648
|
|
|
sort_args = np.argsort(xp, axis=axis) |
|
649
|
|
|
sort_x = np.argsort(x) |
|
650
|
|
|
|
|
651
|
|
|
# indices for sorting |
|
652
|
|
|
sorter = broadcast_indices(xp, sort_args, ndim, axis) |
|
653
|
|
|
|
|
654
|
|
|
# sort xp |
|
655
|
|
|
xp = xp[sorter] |
|
656
|
|
|
# Ensure pressure in increasing order |
|
657
|
|
|
variables = [arr[sorter] for arr in args] |
|
658
|
|
|
|
|
659
|
|
|
# Make x broadcast with xp |
|
660
|
|
|
x_array = x[sort_x] |
|
661
|
|
|
expand = [np.newaxis] * ndim |
|
662
|
|
|
expand[axis] = slice(None) |
|
663
|
|
|
x_array = x_array[expand] |
|
664
|
|
|
|
|
665
|
|
|
# Calculate value above interpolated value |
|
666
|
|
|
minv = np.apply_along_axis(np.searchsorted, axis, xp, x[sort_x]) |
|
667
|
|
|
minv2 = np.copy(minv) |
|
668
|
|
|
|
|
669
|
|
|
# If fill_value is none and data is out of bounds, raise value error |
|
670
|
|
|
if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None: |
|
671
|
|
|
raise ValueError('Interpolation point out of data bounds encountered') |
|
672
|
|
|
|
|
673
|
|
|
# Warn if interpolated values are outside data bounds, will make these the values |
|
674
|
|
|
# at end of data range. |
|
675
|
|
|
if np.max(minv) == xp.shape[axis]: |
|
676
|
|
|
warnings.warn('Interpolation point out of data bounds encountered') |
|
677
|
|
|
minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1 |
|
678
|
|
|
if np.min(minv) == 0: |
|
679
|
|
|
minv2[minv == 0] = 1 |
|
680
|
|
|
|
|
681
|
|
|
# Get indices for broadcasting arrays |
|
682
|
|
|
above = broadcast_indices(xp, minv2, ndim, axis) |
|
683
|
|
|
below = broadcast_indices(xp, minv2 - 1, ndim, axis) |
|
684
|
|
|
|
|
685
|
|
|
if np.any(x_array < xp[below]): |
|
686
|
|
|
warnings.warn('Interpolation point out of data bounds encountered') |
|
687
|
|
|
|
|
688
|
|
|
# Create empty output list |
|
689
|
|
|
ret = [] |
|
690
|
|
|
|
|
691
|
|
|
# Calculate interpolation for each variable |
|
692
|
|
|
for var in variables: |
|
693
|
|
|
var_interp = var[below] + ((x_array - xp[below]) / |
|
694
|
|
|
(xp[above] - xp[below])) * (var[above] - |
|
695
|
|
|
var[below]) |
|
696
|
|
|
|
|
697
|
|
|
# Set points out of bounds to fill value. |
|
698
|
|
|
var_interp[minv == xp.shape[axis]] = fill_value |
|
699
|
|
|
var_interp[x_array < xp[below]] = fill_value |
|
700
|
|
|
|
|
701
|
|
|
# Check for input points in decreasing order and return output to match. |
|
702
|
|
|
if x[0] > x[-1]: |
|
703
|
|
|
var_interp = np.swapaxes(np.swapaxes(var_interp, 0, axis)[::-1], 0, axis) |
|
704
|
|
|
# Output to list |
|
705
|
|
|
ret.append(var_interp) |
|
706
|
|
|
if len(ret) == 1: |
|
707
|
|
|
return ret[0] |
|
708
|
|
|
else: |
|
709
|
|
|
return ret |
|
710
|
|
|
|
|
711
|
|
|
|
|
712
|
|
|
def broadcast_indices(x, minv, ndim, axis): |
|
713
|
|
|
"""Calculate index values to properly broadcast index array within data array. |
|
714
|
|
|
|
|
715
|
|
|
See usage in interp. |
|
716
|
|
|
""" |
|
717
|
|
|
ret = [] |
|
718
|
|
|
for dim in range(ndim): |
|
719
|
|
|
if dim == axis: |
|
720
|
|
|
ret.append(minv) |
|
721
|
|
|
else: |
|
722
|
|
|
broadcast_slice = [np.newaxis] * ndim |
|
723
|
|
|
broadcast_slice[dim] = slice(None) |
|
724
|
|
|
dim_inds = np.arange(x.shape[dim]) |
|
725
|
|
|
ret.append(dim_inds[broadcast_slice]) |
|
726
|
|
|
return ret |
|
727
|
|
|
|
|
728
|
|
|
|
|
729
|
|
|
@exporter.export |
|
730
|
|
|
@units.wraps(None, ('=A', '=A')) |
|
731
|
|
|
def log_interp(x, xp, *args, **kwargs): |
|
732
|
|
|
r"""Interpolates data with logarithmic x-scale over a specified axis. |
|
733
|
|
|
|
|
734
|
|
|
Interpolation on a logarithmic x-scale for interpolation values in pressure coordintates. |
|
735
|
|
|
|
|
736
|
|
|
Parameters |
|
737
|
|
|
---------- |
|
738
|
|
|
x : array-like |
|
739
|
|
|
1-D array of desired interpolated values. |
|
740
|
|
|
|
|
741
|
|
|
xp : array-like |
|
742
|
|
|
The x-coordinates of the data points. |
|
743
|
|
|
|
|
744
|
|
|
args : array-like |
|
745
|
|
|
The data to be interpolated. Can be multiple arguments, all must be the same shape as |
|
746
|
|
|
xp. |
|
747
|
|
|
|
|
748
|
|
|
axis : int, optional |
|
749
|
|
|
The axis to interpolate over. Defaults to 0. |
|
750
|
|
|
|
|
751
|
|
|
fill_value: float, optional |
|
752
|
|
|
Specify handling of interpolation points out of data bounds. If None, will return |
|
753
|
|
|
ValueError if points are out of bounds. Defaults to nan. |
|
754
|
|
|
|
|
755
|
|
|
Returns |
|
756
|
|
|
------- |
|
757
|
|
|
array-like |
|
758
|
|
|
Interpolated values for each point with coordinates sorted in ascending order. |
|
759
|
|
|
|
|
760
|
|
|
Examples |
|
761
|
|
|
-------- |
|
762
|
|
|
>>> x_log = np.array([1e3, 1e4, 1e5, 1e6]) |
|
763
|
|
|
>>> y_log = np.log(x_log) * 2 + 3 |
|
764
|
|
|
>>> x_interp = np.array([5e3, 5e4, 5e5]) |
|
765
|
|
|
>>> metpy.calc.log_interp(x_interp, x_log, y_log) |
|
766
|
|
|
array([ 20.03438638, 24.63955657, 29.24472675]) |
|
767
|
|
|
|
|
768
|
|
|
Notes |
|
769
|
|
|
----- |
|
770
|
|
|
xp and args must be the same shape. |
|
771
|
|
|
|
|
772
|
|
|
""" |
|
773
|
|
|
# Pull out kwargs |
|
774
|
|
|
fill_value = kwargs.pop('fill_value', np.nan) |
|
775
|
|
|
axis = kwargs.pop('axis', 0) |
|
776
|
|
|
|
|
777
|
|
|
# Log x and xp |
|
778
|
|
|
log_x = np.log(x) |
|
779
|
|
|
log_xp = np.log(xp) |
|
780
|
|
|
return interp(log_x, log_xp, *args, axis=axis, fill_value=fill_value) |
|
781
|
|
|
|
|
782
|
|
|
|
|
783
|
|
|
def _greater_or_close(a, value, **kwargs): |
|
784
|
|
|
r"""Compare values for greater or close to boolean masks. |
|
785
|
|
|
|
|
786
|
|
|
Returns a boolean mask for values greater than or equal to a target within a specified |
|
787
|
|
|
absolute or relative tolerance (as in :func:`numpy.isclose`). |
|
788
|
|
|
|
|
789
|
|
|
Parameters |
|
790
|
|
|
---------- |
|
791
|
|
|
a : array-like |
|
792
|
|
|
Array of values to be compared |
|
793
|
|
|
value : float |
|
794
|
|
|
Comparison value |
|
795
|
|
|
|
|
796
|
|
|
Returns |
|
797
|
|
|
------- |
|
798
|
|
|
array-like |
|
799
|
|
|
Boolean array where values are greater than or nearly equal to value. |
|
800
|
|
|
|
|
801
|
|
|
""" |
|
802
|
|
|
return np.greater(a, value) | np.isclose(a, value, **kwargs) |
|
803
|
|
|
|
|
804
|
|
|
|
|
805
|
|
|
def _less_or_close(a, value, **kwargs): |
|
806
|
|
|
r"""Compare values for less or close to boolean masks. |
|
807
|
|
|
|
|
808
|
|
|
Returns a boolean mask for values less than or equal to a target within a specified |
|
809
|
|
|
absolute or relative tolerance (as in :func:`numpy.isclose`). |
|
810
|
|
|
|
|
811
|
|
|
Parameters |
|
812
|
|
|
---------- |
|
813
|
|
|
a : array-like |
|
814
|
|
|
Array of values to be compared |
|
815
|
|
|
value : float |
|
816
|
|
|
Comparison value |
|
817
|
|
|
|
|
818
|
|
|
Returns |
|
819
|
|
|
------- |
|
820
|
|
|
array-like |
|
821
|
|
|
Boolean array where values are less than or nearly equal to value. |
|
822
|
|
|
|
|
823
|
|
|
""" |
|
824
|
|
|
return np.less(a, value) | np.isclose(a, value, **kwargs) |
|
825
|
|
|
|