_get_bound_pressure_height()   F
last analyzed

Complexity

Conditions 18

Size

Total Lines 102

Duplication

Lines 0
Ratio 0 %

Importance

Changes 4
Bugs 0 Features 0
Metric Value
cc 18
c 4
b 0
f 0
dl 0
loc 102
rs 2

How to fix   Long Method    Complexity   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

Complexity

Complex classes like _get_bound_pressure_height() often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.

Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.

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