Completed
Push — master ( 41e9c7...323290 )
by Ryan
19s
created

get_layer_heights()   F

Complexity

Conditions 10

Size

Total Lines 89

Duplication

Lines 0
Ratio 0 %

Importance

Changes 1
Bugs 0 Features 0
Metric Value
cc 10
c 1
b 0
f 0
dl 0
loc 89
rs 3.4615

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_layer_heights() 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
    # 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