Completed
Pull Request — master (#577)
by
unknown
49s
created

get_layer_heights()   F

Complexity

Conditions 10

Size

Total Lines 87

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 87
rs 3.5294

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