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
|
|
|
# Make sure pressure is monotonically decreasing |
327
|
|
|
sort_inds = np.argsort(pressure)[::-1] |
328
|
|
|
pressure = pressure[sort_inds] |
329
|
|
|
if heights is not None: |
330
|
|
|
heights = heights[sort_inds] |
331
|
|
|
|
332
|
|
|
# Bound is given in pressure |
333
|
|
|
if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}: |
334
|
|
|
# If the bound is in the pressure data, we know the pressure bound exactly |
335
|
|
|
if bound in pressure: |
336
|
|
|
bound_pressure = bound |
337
|
|
|
# If we have heights, we know the exact height value, otherwise return standard |
338
|
|
|
# atmosphere height for the pressure |
339
|
|
|
if heights is not None: |
340
|
|
|
bound_height = heights[pressure == bound_pressure] |
341
|
|
|
else: |
342
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
343
|
|
|
# If bound is not in the data, return the nearest or interpolated values |
344
|
|
|
else: |
345
|
|
|
if interpolate: |
346
|
|
|
bound_pressure = bound # Use the user specified bound |
347
|
|
|
if heights is not None: # Interpolate heights from the height data |
348
|
|
|
bound_height = log_interp(bound_pressure, pressure, heights) |
349
|
|
|
else: # If not heights given, use the standard atmosphere |
350
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
351
|
|
|
else: # No interpolation, find the closest values |
352
|
|
|
idx = (np.abs(pressure - bound)).argmin() |
353
|
|
|
bound_pressure = pressure[idx] |
354
|
|
|
if heights is not None: |
355
|
|
|
bound_height = heights[idx] |
356
|
|
|
else: |
357
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
358
|
|
|
|
359
|
|
|
# Bound is given in height |
360
|
|
|
elif bound.dimensionality == {'[length]': 1.0}: |
361
|
|
|
# If there is height data, see if we have the bound or need to interpolate/find nearest |
362
|
|
|
if heights is not None: |
363
|
|
|
if bound in heights: # Bound is in the height data |
364
|
|
|
bound_height = bound |
365
|
|
|
bound_pressure = pressure[heights == bound] |
366
|
|
|
else: # Bound is not in the data |
367
|
|
|
if interpolate: |
368
|
|
|
bound_height = bound |
369
|
|
|
|
370
|
|
|
# Need to cast back to the input type since interp (up to at least numpy |
371
|
|
|
# 1.13 always returns float64. This can cause upstream users problems, |
372
|
|
|
# resulting in something like np.append() to upcast. |
373
|
|
|
bound_pressure = np.interp(np.atleast_1d(bound), heights, |
374
|
|
|
pressure).astype(bound.dtype) * pressure.units |
375
|
|
|
else: |
376
|
|
|
idx = (np.abs(heights - bound)).argmin() |
377
|
|
|
bound_pressure = pressure[idx] |
378
|
|
|
bound_height = heights[idx] |
379
|
|
|
else: # Don't have heights, so assume a standard atmosphere |
380
|
|
|
bound_height = bound |
381
|
|
|
bound_pressure = height_to_pressure_std(bound) |
382
|
|
|
# If interpolation is on, this is all we need, if not, we need to go back and |
383
|
|
|
# find the pressure closest to this and refigure the bounds |
384
|
|
|
if not interpolate: |
385
|
|
|
idx = (np.abs(pressure - bound_pressure)).argmin() |
386
|
|
|
bound_pressure = pressure[idx] |
387
|
|
|
bound_height = pressure_to_height_std(bound_pressure) |
388
|
|
|
|
389
|
|
|
# Bound has invalid units |
390
|
|
|
else: |
391
|
|
|
raise ValueError('Bound must be specified in units of length or pressure.') |
392
|
|
|
|
393
|
|
|
# If the bound is out of the range of the data, we shouldn't extrapolate |
394
|
|
|
if not (_greater_or_close(bound_pressure, np.nanmin(pressure) * pressure.units) and |
395
|
|
|
_less_or_close(bound_pressure, np.nanmax(pressure) * pressure.units)): |
396
|
|
|
raise ValueError('Specified bound is outside pressure range.') |
397
|
|
|
if heights is not None: |
398
|
|
|
if not (_less_or_close(bound_height, np.nanmax(heights) * heights.units) and |
399
|
|
|
_greater_or_close(bound_height, np.nanmin(heights) * heights.units)): |
400
|
|
|
raise ValueError('Specified bound is outside height range.') |
401
|
|
|
|
402
|
|
|
return bound_pressure, bound_height |
403
|
|
|
|
404
|
|
|
|
405
|
|
|
@exporter.export |
406
|
|
|
@check_units('[length]') |
407
|
|
|
def get_layer_heights(heights, depth, *args, **kwargs): |
408
|
|
|
"""Return an atmospheric layer from upper air data with the requested bottom and depth. |
409
|
|
|
|
410
|
|
|
This function will subset an upper air dataset to contain only the specified layer using |
411
|
|
|
the heights only. |
412
|
|
|
|
413
|
|
|
Parameters |
414
|
|
|
---------- |
415
|
|
|
heights : array-like |
416
|
|
|
Atmospheric heights |
417
|
|
|
depth : `pint.Quantity` |
418
|
|
|
The thickness of the layer |
419
|
|
|
*args : array-like |
420
|
|
|
Atmospheric variable(s) measured at the given pressures |
421
|
|
|
bottom : `pint.Quantity`, optional |
422
|
|
|
The bottom of the layer |
423
|
|
|
interpolate : bool, optional |
424
|
|
|
Interpolate the top and bottom points if they are not in the given data. Defaults |
425
|
|
|
to True. |
426
|
|
|
with_agl : bool, optional |
427
|
|
|
Returns the heights as above ground level by subtracting the minimum height in the |
428
|
|
|
provided heights. Defaults to False. |
429
|
|
|
|
430
|
|
|
Returns |
431
|
|
|
------- |
432
|
|
|
`pint.Quantity, pint.Quantity` |
433
|
|
|
The height and data variables of the layer |
434
|
|
|
|
435
|
|
|
""" |
436
|
|
|
bottom = kwargs.pop('bottom', None) |
437
|
|
|
interpolate = kwargs.pop('interpolate', True) |
438
|
|
|
with_agl = kwargs.pop('with_agl', False) |
439
|
|
|
|
440
|
|
|
# Make sure pressure and datavars are the same length |
441
|
|
|
for datavar in args: |
442
|
|
|
if len(heights) != len(datavar): |
443
|
|
|
raise ValueError('Height and data variables must have the same length.') |
444
|
|
|
|
445
|
|
|
# If we want things in AGL, subtract the minimum height from all height values |
446
|
|
|
if with_agl: |
447
|
|
|
sfc_height = np.min(heights) |
448
|
|
|
heights -= sfc_height |
449
|
|
|
|
450
|
|
|
# If the bottom is not specified, make it the surface |
451
|
|
|
if bottom is None: |
452
|
|
|
bottom = heights[0] |
453
|
|
|
|
454
|
|
|
# Make heights and arguments base units |
455
|
|
|
heights = heights.to_base_units() |
456
|
|
|
bottom = bottom.to_base_units() |
457
|
|
|
|
458
|
|
|
# Calculate the top of the layer |
459
|
|
|
top = bottom + depth |
460
|
|
|
|
461
|
|
|
ret = [] # returned data variables in layer |
462
|
|
|
|
463
|
|
|
# Ensure heights are sorted in ascending order |
464
|
|
|
sort_inds = np.argsort(heights) |
465
|
|
|
heights = heights[sort_inds] |
466
|
|
|
|
467
|
|
|
# Mask based on top and bottom |
468
|
|
|
inds = _greater_or_close(heights, bottom) & _less_or_close(heights, top) |
469
|
|
|
heights_interp = heights[inds] |
470
|
|
|
|
471
|
|
|
# Interpolate heights at bounds if necessary and sort |
472
|
|
|
if interpolate: |
473
|
|
|
# If we don't have the bottom or top requested, append them |
474
|
|
|
if top not in heights_interp: |
475
|
|
|
heights_interp = np.sort(np.append(heights_interp, top)) * heights.units |
476
|
|
|
if bottom not in heights_interp: |
477
|
|
|
heights_interp = np.sort(np.append(heights_interp, bottom)) * heights.units |
478
|
|
|
|
479
|
|
|
ret.append(heights_interp) |
480
|
|
|
|
481
|
|
|
for datavar in args: |
482
|
|
|
# Ensure that things are sorted in ascending order |
483
|
|
|
datavar = datavar[sort_inds] |
484
|
|
|
|
485
|
|
|
if interpolate: |
486
|
|
|
# Interpolate for the possibly missing bottom/top values |
487
|
|
|
datavar_interp = interp(heights_interp, heights, datavar) |
488
|
|
|
datavar = datavar_interp |
489
|
|
|
else: |
490
|
|
|
datavar = datavar[inds] |
491
|
|
|
|
492
|
|
|
ret.append(datavar) |
493
|
|
|
return ret |
494
|
|
|
|
495
|
|
|
|
496
|
|
|
@exporter.export |
497
|
|
|
@check_units('[pressure]') |
498
|
|
|
def get_layer(pressure, *args, **kwargs): |
499
|
|
|
r"""Return an atmospheric layer from upper air data with the requested bottom and depth. |
500
|
|
|
|
501
|
|
|
This function will subset an upper air dataset to contain only the specified layer. The |
502
|
|
|
bottom of the layer can be specified with a pressure or height above the surface |
503
|
|
|
pressure. The bottom defaults to the surface pressure. The depth of the layer can be |
504
|
|
|
specified in terms of pressure or height above the bottom of the layer. If the top and |
505
|
|
|
bottom of the layer are not in the data, they are interpolated by default. |
506
|
|
|
|
507
|
|
|
Parameters |
508
|
|
|
---------- |
509
|
|
|
pressure : array-like |
510
|
|
|
Atmospheric pressure profile |
511
|
|
|
*args : array-like |
512
|
|
|
Atmospheric variable(s) measured at the given pressures |
513
|
|
|
heights: array-like, optional |
514
|
|
|
Atmospheric heights corresponding to the given pressures. Defaults to using |
515
|
|
|
heights calculated from ``p`` assuming a standard atmosphere. |
516
|
|
|
bottom : `pint.Quantity`, optional |
517
|
|
|
The bottom of the layer as a pressure or height above the surface pressure. Defaults |
518
|
|
|
to the lowest pressure or height given. |
519
|
|
|
depth : `pint.Quantity`, optional |
520
|
|
|
The thickness of the layer as a pressure or height above the bottom of the layer. |
521
|
|
|
Defaults to 100 hPa. |
522
|
|
|
interpolate : bool, optional |
523
|
|
|
Interpolate the top and bottom points if they are not in the given data. Defaults |
524
|
|
|
to True. |
525
|
|
|
|
526
|
|
|
Returns |
527
|
|
|
------- |
528
|
|
|
`pint.Quantity, pint.Quantity` |
529
|
|
|
The pressure and data variables of the layer |
530
|
|
|
|
531
|
|
|
""" |
532
|
|
|
# Pop off keyword arguments |
533
|
|
|
heights = kwargs.pop('heights', None) |
534
|
|
|
bottom = kwargs.pop('bottom', None) |
535
|
|
|
depth = kwargs.pop('depth', 100 * units.hPa) |
536
|
|
|
interpolate = kwargs.pop('interpolate', True) |
537
|
|
|
|
538
|
|
|
# If we get the depth kwarg, but it's None, set it to the default as well |
539
|
|
|
if depth is None: |
540
|
|
|
depth = 100 * units.hPa |
541
|
|
|
|
542
|
|
|
# Make sure pressure and datavars are the same length |
543
|
|
|
for datavar in args: |
544
|
|
|
if len(pressure) != len(datavar): |
545
|
|
|
raise ValueError('Pressure and data variables must have the same length.') |
546
|
|
|
|
547
|
|
|
# If the bottom is not specified, make it the surface pressure |
548
|
|
|
if bottom is None: |
549
|
|
|
bottom = np.nanmax(pressure) * pressure.units |
550
|
|
|
|
551
|
|
|
bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom, |
552
|
|
|
heights=heights, |
553
|
|
|
interpolate=interpolate) |
554
|
|
|
|
555
|
|
|
# Calculate the top if whatever units depth is in |
556
|
|
|
if depth.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}: |
557
|
|
|
top = bottom_pressure - depth |
558
|
|
|
elif depth.dimensionality == {'[length]': 1}: |
559
|
|
|
top = bottom_height + depth |
560
|
|
|
else: |
561
|
|
|
raise ValueError('Depth must be specified in units of length or pressure') |
562
|
|
|
|
563
|
|
|
top_pressure, _ = _get_bound_pressure_height(pressure, top, heights=heights, |
564
|
|
|
interpolate=interpolate) |
565
|
|
|
|
566
|
|
|
ret = [] # returned data variables in layer |
567
|
|
|
|
568
|
|
|
# Ensure pressures are sorted in ascending order |
569
|
|
|
sort_inds = np.argsort(pressure) |
570
|
|
|
pressure = pressure[sort_inds] |
571
|
|
|
|
572
|
|
|
# Mask based on top and bottom pressure |
573
|
|
|
inds = (_less_or_close(pressure, bottom_pressure) & |
574
|
|
|
_greater_or_close(pressure, top_pressure)) |
575
|
|
|
p_interp = pressure[inds] |
576
|
|
|
|
577
|
|
|
# Interpolate pressures at bounds if necessary and sort |
578
|
|
|
if interpolate: |
579
|
|
|
# If we don't have the bottom or top requested, append them |
580
|
|
|
if not np.any(np.isclose(top_pressure, p_interp)): |
581
|
|
|
p_interp = np.sort(np.append(p_interp, top_pressure)) * pressure.units |
582
|
|
|
if not np.any(np.isclose(bottom_pressure, p_interp)): |
583
|
|
|
p_interp = np.sort(np.append(p_interp, bottom_pressure)) * pressure.units |
584
|
|
|
|
585
|
|
|
ret.append(p_interp[::-1]) |
586
|
|
|
|
587
|
|
|
for datavar in args: |
588
|
|
|
# Ensure that things are sorted in ascending order |
589
|
|
|
datavar = datavar[sort_inds] |
590
|
|
|
|
591
|
|
|
if interpolate: |
592
|
|
|
# Interpolate for the possibly missing bottom/top values |
593
|
|
|
datavar_interp = log_interp(p_interp, pressure, datavar) |
594
|
|
|
datavar = datavar_interp |
595
|
|
|
else: |
596
|
|
|
datavar = datavar[inds] |
597
|
|
|
|
598
|
|
|
ret.append(datavar[::-1]) |
599
|
|
|
return ret |
600
|
|
|
|
601
|
|
|
|
602
|
|
|
@exporter.export |
603
|
|
|
@units.wraps(None, ('=A', '=A')) |
604
|
|
|
def interp(x, xp, *args, **kwargs): |
605
|
|
|
r"""Interpolates data with any shape over a specified axis. |
606
|
|
|
|
607
|
|
|
Interpolation over a specified axis for arrays of any shape. |
608
|
|
|
|
609
|
|
|
Parameters |
610
|
|
|
---------- |
611
|
|
|
x : array-like |
612
|
|
|
1-D array of desired interpolated values. |
613
|
|
|
|
614
|
|
|
xp : array-like |
615
|
|
|
The x-coordinates of the data points. |
616
|
|
|
|
617
|
|
|
args : array-like |
618
|
|
|
The data to be interpolated. Can be multiple arguments, all must be the same shape as |
619
|
|
|
xp. |
620
|
|
|
|
621
|
|
|
axis : int, optional |
622
|
|
|
The axis to interpolate over. Defaults to 0. |
623
|
|
|
|
624
|
|
|
fill_value: float, optional |
625
|
|
|
Specify handling of interpolation points out of data bounds. If None, will return |
626
|
|
|
ValueError if points are out of bounds. Defaults to nan. |
627
|
|
|
|
628
|
|
|
Returns |
629
|
|
|
------- |
630
|
|
|
array-like |
631
|
|
|
Interpolated values for each point with coordinates sorted in ascending order. |
632
|
|
|
|
633
|
|
|
Examples |
634
|
|
|
-------- |
635
|
|
|
>>> x = np.array([1., 2., 3., 4.]) |
636
|
|
|
>>> y = np.array([1., 2., 3., 4.]) |
637
|
|
|
>>> x_interp = np.array([2.5, 3.5]) |
638
|
|
|
>>> metpy.calc.interp(x_interp, x, y) |
639
|
|
|
array([ 2.5, 3.5]) |
640
|
|
|
|
641
|
|
|
Notes |
642
|
|
|
----- |
643
|
|
|
xp and args must be the same shape. |
644
|
|
|
|
645
|
|
|
""" |
646
|
|
|
# Pull out keyword args |
647
|
|
|
fill_value = kwargs.pop('fill_value', np.nan) |
648
|
|
|
axis = kwargs.pop('axis', 0) |
649
|
|
|
|
650
|
|
|
# Make x an array |
651
|
|
|
x = np.asanyarray(x).reshape(-1) |
652
|
|
|
|
653
|
|
|
# Save number of dimensions in xp |
654
|
|
|
ndim = xp.ndim |
655
|
|
|
|
656
|
|
|
# Sort input data |
657
|
|
|
sort_args = np.argsort(xp, axis=axis) |
658
|
|
|
sort_x = np.argsort(x) |
659
|
|
|
|
660
|
|
|
# indices for sorting |
661
|
|
|
sorter = broadcast_indices(xp, sort_args, ndim, axis) |
662
|
|
|
|
663
|
|
|
# sort xp |
664
|
|
|
xp = xp[sorter] |
665
|
|
|
# Ensure pressure in increasing order |
666
|
|
|
variables = [arr[sorter] for arr in args] |
667
|
|
|
|
668
|
|
|
# Make x broadcast with xp |
669
|
|
|
x_array = x[sort_x] |
670
|
|
|
expand = [np.newaxis] * ndim |
671
|
|
|
expand[axis] = slice(None) |
672
|
|
|
x_array = x_array[expand] |
673
|
|
|
|
674
|
|
|
# Calculate value above interpolated value |
675
|
|
|
minv = np.apply_along_axis(np.searchsorted, axis, xp, x[sort_x]) |
676
|
|
|
minv2 = np.copy(minv) |
677
|
|
|
|
678
|
|
|
# If fill_value is none and data is out of bounds, raise value error |
679
|
|
|
if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None: |
680
|
|
|
raise ValueError('Interpolation point out of data bounds encountered') |
681
|
|
|
|
682
|
|
|
# Warn if interpolated values are outside data bounds, will make these the values |
683
|
|
|
# at end of data range. |
684
|
|
|
if np.max(minv) == xp.shape[axis]: |
685
|
|
|
warnings.warn('Interpolation point out of data bounds encountered') |
686
|
|
|
minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1 |
687
|
|
|
if np.min(minv) == 0: |
688
|
|
|
minv2[minv == 0] = 1 |
689
|
|
|
|
690
|
|
|
# Get indices for broadcasting arrays |
691
|
|
|
above = broadcast_indices(xp, minv2, ndim, axis) |
692
|
|
|
below = broadcast_indices(xp, minv2 - 1, ndim, axis) |
693
|
|
|
|
694
|
|
|
if np.any(x_array < xp[below]): |
695
|
|
|
warnings.warn('Interpolation point out of data bounds encountered') |
696
|
|
|
|
697
|
|
|
# Create empty output list |
698
|
|
|
ret = [] |
699
|
|
|
|
700
|
|
|
# Calculate interpolation for each variable |
701
|
|
|
for var in variables: |
702
|
|
|
var_interp = var[below] + ((x_array - xp[below]) / |
703
|
|
|
(xp[above] - xp[below])) * (var[above] - |
704
|
|
|
var[below]) |
705
|
|
|
|
706
|
|
|
# Set points out of bounds to fill value. |
707
|
|
|
var_interp[minv == xp.shape[axis]] = fill_value |
708
|
|
|
var_interp[x_array < xp[below]] = fill_value |
709
|
|
|
|
710
|
|
|
# Check for input points in decreasing order and return output to match. |
711
|
|
|
if x[0] > x[-1]: |
712
|
|
|
var_interp = np.swapaxes(np.swapaxes(var_interp, 0, axis)[::-1], 0, axis) |
713
|
|
|
# Output to list |
714
|
|
|
ret.append(var_interp) |
715
|
|
|
if len(ret) == 1: |
716
|
|
|
return ret[0] |
717
|
|
|
else: |
718
|
|
|
return ret |
719
|
|
|
|
720
|
|
|
|
721
|
|
|
def broadcast_indices(x, minv, ndim, axis): |
722
|
|
|
"""Calculate index values to properly broadcast index array within data array. |
723
|
|
|
|
724
|
|
|
See usage in interp. |
725
|
|
|
""" |
726
|
|
|
ret = [] |
727
|
|
|
for dim in range(ndim): |
728
|
|
|
if dim == axis: |
729
|
|
|
ret.append(minv) |
730
|
|
|
else: |
731
|
|
|
broadcast_slice = [np.newaxis] * ndim |
732
|
|
|
broadcast_slice[dim] = slice(None) |
733
|
|
|
dim_inds = np.arange(x.shape[dim]) |
734
|
|
|
ret.append(dim_inds[broadcast_slice]) |
735
|
|
|
return ret |
736
|
|
|
|
737
|
|
|
|
738
|
|
|
@exporter.export |
739
|
|
|
@units.wraps(None, ('=A', '=A')) |
740
|
|
|
def log_interp(x, xp, *args, **kwargs): |
741
|
|
|
r"""Interpolates data with logarithmic x-scale over a specified axis. |
742
|
|
|
|
743
|
|
|
Interpolation on a logarithmic x-scale for interpolation values in pressure coordintates. |
744
|
|
|
|
745
|
|
|
Parameters |
746
|
|
|
---------- |
747
|
|
|
x : array-like |
748
|
|
|
1-D array of desired interpolated values. |
749
|
|
|
|
750
|
|
|
xp : array-like |
751
|
|
|
The x-coordinates of the data points. |
752
|
|
|
|
753
|
|
|
args : array-like |
754
|
|
|
The data to be interpolated. Can be multiple arguments, all must be the same shape as |
755
|
|
|
xp. |
756
|
|
|
|
757
|
|
|
axis : int, optional |
758
|
|
|
The axis to interpolate over. Defaults to 0. |
759
|
|
|
|
760
|
|
|
fill_value: float, optional |
761
|
|
|
Specify handling of interpolation points out of data bounds. If None, will return |
762
|
|
|
ValueError if points are out of bounds. Defaults to nan. |
763
|
|
|
|
764
|
|
|
Returns |
765
|
|
|
------- |
766
|
|
|
array-like |
767
|
|
|
Interpolated values for each point with coordinates sorted in ascending order. |
768
|
|
|
|
769
|
|
|
Examples |
770
|
|
|
-------- |
771
|
|
|
>>> x_log = np.array([1e3, 1e4, 1e5, 1e6]) |
772
|
|
|
>>> y_log = np.log(x_log) * 2 + 3 |
773
|
|
|
>>> x_interp = np.array([5e3, 5e4, 5e5]) |
774
|
|
|
>>> metpy.calc.log_interp(x_interp, x_log, y_log) |
775
|
|
|
array([ 20.03438638, 24.63955657, 29.24472675]) |
776
|
|
|
|
777
|
|
|
Notes |
778
|
|
|
----- |
779
|
|
|
xp and args must be the same shape. |
780
|
|
|
|
781
|
|
|
""" |
782
|
|
|
# Pull out kwargs |
783
|
|
|
fill_value = kwargs.pop('fill_value', np.nan) |
784
|
|
|
axis = kwargs.pop('axis', 0) |
785
|
|
|
|
786
|
|
|
# Log x and xp |
787
|
|
|
log_x = np.log(x) |
788
|
|
|
log_xp = np.log(xp) |
789
|
|
|
return interp(log_x, log_xp, *args, axis=axis, fill_value=fill_value) |
790
|
|
|
|
791
|
|
|
|
792
|
|
|
def _greater_or_close(a, value, **kwargs): |
793
|
|
|
r"""Compare values for greater or close to boolean masks. |
794
|
|
|
|
795
|
|
|
Returns a boolean mask for values greater than or equal to a target within a specified |
796
|
|
|
absolute or relative tolerance (as in :func:`numpy.isclose`). |
797
|
|
|
|
798
|
|
|
Parameters |
799
|
|
|
---------- |
800
|
|
|
a : array-like |
801
|
|
|
Array of values to be compared |
802
|
|
|
value : float |
803
|
|
|
Comparison value |
804
|
|
|
|
805
|
|
|
Returns |
806
|
|
|
------- |
807
|
|
|
array-like |
808
|
|
|
Boolean array where values are greater than or nearly equal to value. |
809
|
|
|
|
810
|
|
|
""" |
811
|
|
|
return np.greater(a, value) | np.isclose(a, value, **kwargs) |
812
|
|
|
|
813
|
|
|
|
814
|
|
|
def _less_or_close(a, value, **kwargs): |
815
|
|
|
r"""Compare values for less or close to boolean masks. |
816
|
|
|
|
817
|
|
|
Returns a boolean mask for values less than or equal to a target within a specified |
818
|
|
|
absolute or relative tolerance (as in :func:`numpy.isclose`). |
819
|
|
|
|
820
|
|
|
Parameters |
821
|
|
|
---------- |
822
|
|
|
a : array-like |
823
|
|
|
Array of values to be compared |
824
|
|
|
value : float |
825
|
|
|
Comparison value |
826
|
|
|
|
827
|
|
|
Returns |
828
|
|
|
------- |
829
|
|
|
array-like |
830
|
|
|
Boolean array where values are less than or nearly equal to value. |
831
|
|
|
|
832
|
|
|
""" |
833
|
|
|
return np.less(a, value) | np.isclose(a, value, **kwargs) |
834
|
|
|
|