_get_bound_pressure_height() - Code Metrics - Inspection of "Log interpolation" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#468)

unknown

created 2017-06-27 19:59 UTC

_get_bound_pressure_height() F

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	2
Bugs	0	Features	0

Metric	Value
cc	17
c	2
b	0
f	0
dl	0
loc	88
rs	2

How to fix Long Method Complexity

# Copyright (c) 2008-2017 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of generally useful calculation tools."""

import functools
import warnings

import numpy as np
import numpy.ma as ma
from scipy.spatial import cKDTree

from . import height_to_pressure_std, pressure_to_height_std
from ..package_tools import Exporter
from ..units import check_units, units

exporter = Exporter(globals())


@exporter.export
def resample_nn_1d(a, centers):
    """Return one-dimensional nearest-neighbor indexes based on user-specified centers.

    Parameters
    ----------
    a : array-like
        1-dimensional array of numeric values from which to
        extract indexes of nearest-neighbors
    centers : array-like
        1-dimensional array of numeric values representing a subset of values to approximate

    Returns
    -------
        An array of indexes representing values closest to given array values

    """
    ix = []
    for center in centers:
        index = (np.abs(a - center)).argmin()
        if index not in ix:
            ix.append(index)
    return ix


@exporter.export
def nearest_intersection_idx(a, b):
    """Determine the index of the point just before two lines with common x values.

    Parameters
    ----------
    a : array-like
        1-dimensional array of y-values for line 1
    b : array-like
        1-dimensional array of y-values for line 2

    Returns
    -------
        An array of indexes representing the index of the values
        just before the intersection(s) of the two lines.

    """
    # Difference in the two y-value sets
    difference = a - b

    # Determine the point just before the intersection of the lines
    # Will return multiple points for multiple intersections
    sign_change_idx, = np.nonzero(np.diff(np.sign(difference)))

    return sign_change_idx


@exporter.export
def find_intersections(x, a, b, direction='all'):
    """Calculate the best estimate of intersection.

    Calculates the best estimates of the intersection of two y-value
    data sets that share a common x-value set.

    Parameters
    ----------
    x : array-like
        1-dimensional array of numeric x-values
    a : array-like
        1-dimensional array of y-values for line 1
    b : array-like
        1-dimensional array of y-values for line 2
    direction : string
        specifies direction of crossing. 'all', 'increasing' (a becoming greater than b),
        or 'decreasing' (b becoming greater than a).

    Returns
    -------
        A tuple (x, y) of array-like with the x and y coordinates of the
        intersections of the lines.

    """
    # Find the index of the points just before the intersection(s)
    nearest_idx = nearest_intersection_idx(a, b)
    next_idx = nearest_idx + 1

    # Determine the sign of the change
    sign_change = np.sign(a[next_idx] - b[next_idx])

    # x-values around each intersection
    _, x0 = _next_non_masked_element(x, nearest_idx)
    _, x1 = _next_non_masked_element(x, next_idx)

    # y-values around each intersection for the first line
    _, a0 = _next_non_masked_element(a, nearest_idx)
    _, a1 = _next_non_masked_element(a, next_idx)

    # y-values around each intersection for the second line
    _, b0 = _next_non_masked_element(b, nearest_idx)
    _, b1 = _next_non_masked_element(b, next_idx)

    # Calculate the x-intersection. This comes from finding the equations of the two lines,
    # one through (x0, a0) and (x1, a1) and the other through (x0, b0) and (x1, b1),
    # finding their intersection, and reducing with a bunch of algebra.
    delta_y0 = a0 - b0
    delta_y1 = a1 - b1
    intersect_x = (delta_y1 * x0 - delta_y0 * x1) / (delta_y1 - delta_y0)

    # Calculate the y-intersection of the lines. Just plug the x above into the equation
    # for the line through the a points. One could solve for y like x above, but this
    # causes weirder unit behavior and seems a little less good numerically.
    intersect_y = ((intersect_x - x0) / (x1 - x0)) * (a1 - a0) + a0

    # Make a mask based on the direction of sign change desired
    if direction == 'increasing':
        mask = sign_change > 0
    elif direction == 'decreasing':
        mask = sign_change < 0
    elif direction == 'all':
        return intersect_x, intersect_y
    else:
        raise ValueError('Unknown option for direction: {0}'.format(str(direction)))
    return intersect_x[mask], intersect_y[mask]


@exporter.export
def interpolate_nans(x, y, kind='linear'):
    """Interpolate NaN values in y.

    Interpolate NaN values in the y dimension. Works with unsorted x values.

    Parameters
    ----------
    x : array-like
        1-dimensional array of numeric x-values
    y : array-like
        1-dimensional array of numeric y-values
    kind : string
        specifies the kind of interpolation x coordinate - 'linear' or 'log'

    Returns
    -------
        An array of the y coordinate data with NaN values interpolated.

    """
    x_sort_args = np.argsort(x)
    x = x[x_sort_args]
    y = y[x_sort_args]
    nans = np.isnan(y)
    if kind == 'linear':
        y[nans] = np.interp(x[nans], x[~nans], y[~nans])
    elif kind == 'log':
        y[nans] = np.interp(np.log(x[nans]), np.log(x[~nans]), y[~nans])
    else:
        raise ValueError('Unknown option for kind: {0}'.format(str(kind)))
    return y[x_sort_args]


def _next_non_masked_element(a, idx):
    """Return the next non masked element of a masked array.

    If an array is masked, return the next non-masked element (if the given index is masked).
    If no other unmasked points are after the given masked point, returns none.

    Parameters
    ----------
    a : array-like
        1-dimensional array of numeric values
    idx : integer
        index of requested element

    Returns
    -------
        Index of next non-masked element and next non-masked element

    """
    try:
        next_idx = idx + a[idx:].mask.argmin()
        if ma.is_masked(a[next_idx]):
            return None, None
        else:
            return next_idx, a[next_idx]
    except (AttributeError, TypeError, IndexError):
        return idx, a[idx]


def delete_masked_points(*arrs):
    """Delete masked points from arrays.

    Takes arrays and removes masked points to help with calculations and plotting.

    Parameters
    ----------
    arrs : one or more array-like
        source arrays

    Returns
    -------
    arrs : one or more array-like
        arrays with masked elements removed

    """
    if any(hasattr(a, 'mask') for a in arrs):
        keep = ~functools.reduce(np.logical_or, (np.ma.getmaskarray(a) for a in arrs))
        return tuple(ma.asarray(a[keep]) for a in arrs)
    else:
        return arrs


@exporter.export
def reduce_point_density(points, radius, priority=None):
    r"""Return a mask to reduce the density of points in irregularly-spaced data.

    This function is used to down-sample a collection of scattered points (e.g. surface
    data), returning a mask that can be used to select the points from one or more arrays
    (e.g. arrays of temperature and dew point). The points selected can be controlled by
    providing an array of ``priority`` values (e.g. rainfall totals to ensure that
    stations with higher precipitation remain in the mask).

    Parameters
    ----------
    points : (N, K) array-like
        N locations of the points in K dimensional space
    radius : float
        minimum radius allowed between points
    priority : (N, K) array-like, optional
        If given, this should have the same shape as ``points``; these values will
        be used to control selection priority for points.

    Returns
    -------
        (N,) array-like of boolean values indicating whether points should be kept. This
        can be used directly to index numpy arrays to return only the desired points.

    Examples
    --------
    >>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.)
    array([ True, False,  True], dtype=bool)
    >>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.,
    ... priority=np.array([0.1, 0.9, 0.3]))
    array([False,  True, False], dtype=bool)

    """
    # Handle 1D input
    if points.ndim < 2:
        points = points.reshape(-1, 1)

    # Make a kd-tree to speed searching of data.
    tree = cKDTree(points)

    # Need to use sorted indices rather than sorting the position
    # so that the keep mask matches *original* order.
    if priority is not None:
        # Need to sort the locations in decreasing priority.
        sorted_indices = np.argsort(priority)[::-1]
    else:
        # Take advantage of iterator nature of range here to avoid making big lists
        sorted_indices = range(len(points))

    # Keep all points initially
    keep = np.ones(len(points), dtype=np.bool)

    # Loop over all the potential points
    for ind in sorted_indices:
        # Only proceed if we haven't already excluded this point
        if keep[ind]:
            # Find the neighbors and eliminate them
            neighbors = tree.query_ball_point(points[ind], radius)
            keep[neighbors] = False

            # We just removed ourselves, so undo that
            keep[ind] = True

    return keep


@exporter.export
@units.wraps(None, ('=A', '=A', None))
def log_interp(x, xp, *args, **kwargs):
    r"""Interpolates data with logarithmic x-scale over a specified axis.

    Interpolation on a logarithmic x-scale for interpolation values in pressure coordintates.

    Parameters
    ----------
    x : array-like
        1-D array of desired interpolated values.

    xp : array-like
        The x-coordinates of the data points.

    args : array-like
        The y-coordinates of the data points, same shape as xp.

    axis : int, optional
        The axis to interpolate over. Defaults to 0.

    fill_value: float, optional
        Specify handling of interpolation points out of data bounds. If None, will return
        ValueError if points are out of bounds. Defaults to nan.

    Returns
    -------
    array-like
        Interpolated values for each point with coordinates sorted in ascending order.

    Examples
     --------
     >>> x_log = np.array([1e3, 1e4, 1e5, 1e6])
     >>> y_log = np.log(x_log) * 2 + 3
     >>> x_interp = np.array([5e3, 5e4, 5e5])
     >>> metpy.calc.log_interp(x_interp, x_log, y_log)
     array([ 20.03438638,  24.63955657,  29.24472675])

    Notes
    -----
    xp and args must be the same shape.

    """
    # Pull out keyword args
    fill_value = kwargs.pop('fill_value', np.nan)
    axis = kwargs.pop('axis', 0)

    # Save number of dimensions in xp
    ndim = xp.ndim

    # Sort input data
    sort_args = np.argsort(xp, axis=axis)
    sort_x = np.argsort(x)

    # indices for sorting
    sorter = broadcast_indices(xp, sort_args, ndim, axis)

    # Ensure pressure in increasing order
    variables = [arr[sorter] for arr in args]

    # Make x broadcast with xp
    x_array = x[sort_x]
    expand = [np.newaxis] * ndim
    expand[axis] = slice(None)
    x_array = x_array[expand]

    # Log x and xp
    log_x_array = np.log(x_array)
    log_xp = np.log(xp[sorter])

    # Calculate value above interpolated value
    minv = np.apply_along_axis(np.searchsorted, axis, xp[sorter], x[sort_x])
    minv2 = np.copy(minv)

    # If extrap is none and data is out of bounds, raise value error
    if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None:
        raise ValueError('Interpolation point out of data bounds encountered')

    # Warn if interpolated values are outside data bounds, will make these the values
    # at end of data range.
    if np.max(minv) == xp.shape[axis]:
        warnings.warn('Interpolation point out of data bounds encountered')
        minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1
    if np.max(minv) == 0:
        warnings.warn('Interpolation point out of data bounds encountered')
        minv2[minv == 0] = 1

    # Get indicies for broadcasting arrays
    above = broadcast_indices(xp, minv2, ndim, axis)
    below = broadcast_indices(xp, minv2 - 1, ndim, axis)

    # Create empty output list
    ret = []

    # Calculate interpolation for each variable
    for var in variables:
        var_interp = var[below] + ((log_x_array - log_xp[below]) /
                                   (log_xp[above] - log_xp[below])) * (var[above] -
                                                                       var[below])

        # set to nan if fill_value = nan
        var_interp[minv == xp.shape[axis]] = fill_value
        var_interp[minv == 0] = fill_value

        # Output to list
        ret.append(var_interp)
    if len(ret) == 1:
        return ret[0]
    else:
        return ret


def broadcast_indices(x, minv, ndim, axis):
    """Calculate index values to properly broadcast index array within data array.

<<<<<<< HEAD
    See usage in log_interp.
    """
    ret = []
    for dim in range(ndim):
        if dim == axis:
            ret.append(minv)
        else:
            broadcast_slice = [np.newaxis] * ndim
            broadcast_slice[dim] = slice(None)
            dim_inds = np.arange(x.shape[dim])
            ret.append(dim_inds[broadcast_slice])
    return ret


def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True):
    """Calculate the bounding pressure and height in a layer.

    Given pressure, optional heights, and a bound, return either the closest pressure/height
    or interpolated pressure/height. If no heights are provided, a standard atmosphere is
    assumed.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Atmospheric pressures
    bound : `pint.Quantity`
        Bound to retrieve (in pressure or height)
    heights : `pint.Quantity`
        Atmospheric heights associated with the pressure levels
    interpolate : boolean
        Interpolate the bound or return the nearest

    Returns
    -------
    `pint.Quantity`
        The bound pressure and height.

    """
    # Bound is given in pressure
    if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
        # If the bound is in the pressure data, we know the pressure bound exactly
        if bound in pressure:
            bound_pressure = bound
            # If we have heights, we know the exact height value, otherwise return standard
            # atmosphere height for the pressure
            if heights is not None:
                bound_height = heights[pressure == bound_pressure]
            else:
                bound_height = pressure_to_height_std(bound_pressure)
        # If bound is not in the data, return the nearest or interpolated values
        else:
            if interpolate:
                bound_pressure = bound  # Use the user specified bound
                if heights is not None:  # Interpolate heights from the height data
                    bound_height = log_interp(bound_pressure, pressure, heights)
                else:  # If not heights given, use the standard atmosphere
                    bound_height = pressure_to_height_std(bound_pressure)
            else:  # No interpolation, find the closest values
                idx = (np.abs(pressure - bound)).argmin()
                bound_pressure = pressure[idx]
                if heights is not None:
                    bound_height = heights[idx]
                else:
                    bound_height = pressure_to_height_std(bound_pressure)

    # Bound is given in height
    elif bound.dimensionality == {'[length]': 1.0}:
        # If there is height data, see if we have the bound or need to interpolate/find nearest
        if heights is not None:
            if bound in heights:  # Bound is in the height data
                bound_height = bound
                bound_pressure = pressure[heights == bound]
            else:  # Bound is not in the data
                if interpolate:
                    bound_height = bound
                    bound_pressure = height_to_pressure_std(bound_height)
                else:
                    idx = (np.abs(heights - bound)).argmin()
                    bound_pressure = pressure[idx]
                    bound_height = heights[idx]
        else:  # Don't have heights, so assume a standard atmosphere
            bound_height = bound
            bound_pressure = height_to_pressure_std(bound)
            # If interpolation is on, this is all we need, if not, we need to go back and
            # find the pressure closest to this and refigure the bounds
            if not interpolate:
                idx = (np.abs(pressure - bound_pressure)).argmin()
                bound_pressure = pressure[idx]
                bound_height = pressure_to_height_std(bound_pressure)

    # Bound has invalid units
    else:
        raise ValueError('Bound must be specified in units of length or pressure.')

    # If the bound is out of the range of the data, we shouldn't extrapolate
    if (bound_pressure < np.min(pressure)) or (bound_pressure > np.max(pressure)):
        raise ValueError('Specified bound is outside pressure range.')
    if heights is not None:
        if (bound_height > np.max(heights)) or (bound_height < np.min(heights)):
            raise ValueError('Specified bound is outside height range.')

    return bound_pressure, bound_height


@exporter.export
@check_units('[pressure]')
def get_layer(p, *args, **kwargs):
    r"""Return an atmospheric layer from upper air data with the requested bottom and depth.

    This function will subset an upper air dataset to contain only the specified layer. The
    bottom of the layer can be specified with a pressure or height above the surface
    pressure. The bottom defaults to the surface pressure. The depth of the layer can be
    specified in terms of pressure or height above the bottom of the layer. If the top and
    bottom of the layer are not in the data, they are interpolated by default.

    Parameters
    ----------
    p : array-like
        Atmospheric pressure profile
    *args : array-like
        Atmospheric variable(s) measured at the given pressures
    heights: array-like
        Atmospheric heights corresponding to the given pressures
    bottom : `pint.Quantity`
        The bottom of the layer as a pressure or height above the surface pressure
    depth : `pint.Quantity`
        The thickness of the layer as a pressure or height above the bottom of the layer
    interpolate : bool
        Interpolate the top and bottom points if they are not in the given data

    Returns
    -------
    `pint.Quantity, pint.Quantity`
        The pressure and data variables of the layer

    """
    # Pop off keyword arguments
    heights = kwargs.pop('heights', None)
    bottom = kwargs.pop('bottom', None)
    depth = kwargs.pop('depth', 100 * units.hPa)
    interpolate = kwargs.pop('interpolate', True)

    # Make sure pressure and datavars are the same length
    for datavar in args:
        if len(p) != len(datavar):
            raise ValueError('Pressure and data variables must have the same length.')

    # If the bottom is not specified, make it the surface pressure
    if bottom is None:
        bottom = p[0]

    bottom_pressure, bottom_height = _get_bound_pressure_height(p, bottom, heights=heights,
                                                                interpolate=interpolate)

    # Calculate the top if whatever units depth is in
    if depth.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
        top = bottom_pressure - depth
    elif depth.dimensionality == {'[length]': 1}:
        top = bottom_height + depth
    else:
        raise ValueError('Depth must be specified in units of length or pressure')

    top_pressure, _ = _get_bound_pressure_height(p, top, heights=heights,
                                                 interpolate=interpolate)

    ret = []  # returned data variables in layer

    # Ensure pressures are sorted in ascending order
    sort_inds = np.argsort(p)
    p = p[sort_inds]

    # Mask based on top and bottom pressure
    inds = (p <= bottom_pressure) & (p >= top_pressure)
    p_interp = p[inds]

    # Interpolate pressures at bounds if necessary and sort
    if interpolate:
        # If we don't have the bottom or top requested, append them
        if top_pressure not in p_interp:
            p_interp = np.sort(np.append(p_interp, top_pressure)) * p.units
        if bottom_pressure not in p_interp:
            p_interp = np.sort(np.append(p_interp, bottom_pressure)) * p.units

    ret.append(p_interp[::-1])

    for datavar in args:
        # Ensure that things are sorted in ascending order
        datavar = datavar[sort_inds]

        if interpolate:
            # Interpolate for the possibly missing bottom/top values
            datavar_interp = log_interp(p_interp, p, datavar)
            datavar = datavar_interp
        else:
            datavar = datavar[inds]

        ret.append(datavar[::-1])

    return ret


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			def find_intersections(x, a, b, direction='all'):
74			"""Calculate the best estimate of intersection.
75
76			Calculates the best estimates of the intersection of two y-value
77			data sets that share a common x-value set.
78
79			Parameters
80			----------
81			x : array-like
82			1-dimensional array of numeric x-values
83			a : array-like
84			1-dimensional array of y-values for line 1
85			b : array-like
86			1-dimensional array of y-values for line 2
87			direction : string
88			specifies direction of crossing. 'all', 'increasing' (a becoming greater than b),
89			or 'decreasing' (b becoming greater than a).
90
91			Returns
92			-------
93			A tuple (x, y) of array-like with the x and y coordinates of the
94			intersections of the lines.
95
96			"""
97			# Find the index of the points just before the intersection(s)
98			nearest_idx = nearest_intersection_idx(a, b)
99			next_idx = nearest_idx + 1
100
101			# Determine the sign of the change
102			sign_change = np.sign(a[next_idx] - b[next_idx])
103
104			# x-values around each intersection
105			_, x0 = _next_non_masked_element(x, nearest_idx)
106			_, x1 = _next_non_masked_element(x, next_idx)
107
108			# y-values around each intersection for the first line
109			_, a0 = _next_non_masked_element(a, nearest_idx)
110			_, a1 = _next_non_masked_element(a, next_idx)
111
112			# y-values around each intersection for the second line
113			_, b0 = _next_non_masked_element(b, nearest_idx)
114			_, b1 = _next_non_masked_element(b, next_idx)
115
116			# Calculate the x-intersection. This comes from finding the equations of the two lines,
117			# one through (x0, a0) and (x1, a1) and the other through (x0, b0) and (x1, b1),
118			# finding their intersection, and reducing with a bunch of algebra.
119			delta_y0 = a0 - b0
120			delta_y1 = a1 - b1
121			intersect_x = (delta_y1 * x0 - delta_y0 * x1) / (delta_y1 - delta_y0)
122
123			# Calculate the y-intersection of the lines. Just plug the x above into the equation
124			# for the line through the a points. One could solve for y like x above, but this
125			# causes weirder unit behavior and seems a little less good numerically.
126			intersect_y = ((intersect_x - x0) / (x1 - x0)) * (a1 - a0) + a0
127
128			# Make a mask based on the direction of sign change desired
129			if direction == 'increasing':
130			mask = sign_change > 0
131			elif direction == 'decreasing':
132			mask = sign_change < 0
133			elif direction == 'all':
134			return intersect_x, intersect_y
135			else:
136			raise ValueError('Unknown option for direction: {0}'.format(str(direction)))
137			return intersect_x[mask], intersect_y[mask]
138
139
140			@exporter.export
141			def interpolate_nans(x, y, kind='linear'):
142			"""Interpolate NaN values in y.
143
144			Interpolate NaN values in the y dimension. Works with unsorted x values.
145
146			Parameters
147			----------
148			x : array-like
149			1-dimensional array of numeric x-values
150			y : array-like
151			1-dimensional array of numeric y-values
152			kind : string
153			specifies the kind of interpolation x coordinate - 'linear' or 'log'
154
155			Returns
156			-------
157			An array of the y coordinate data with NaN values interpolated.
158
159			"""
160			x_sort_args = np.argsort(x)
161			x = x[x_sort_args]
162			y = y[x_sort_args]
163			nans = np.isnan(y)
164			if kind == 'linear':
165			y[nans] = np.interp(x[nans], x[~nans], y[~nans])
166			elif kind == 'log':
167			y[nans] = np.interp(np.log(x[nans]), np.log(x[~nans]), y[~nans])
168			else:
169			raise ValueError('Unknown option for kind: {0}'.format(str(kind)))
170			return y[x_sort_args]
171
172
173			def _next_non_masked_element(a, idx):
174			"""Return the next non masked element of a masked array.
175
176			If an array is masked, return the next non-masked element (if the given index is masked).
177			If no other unmasked points are after the given masked point, returns none.
178
179			Parameters
180			----------
181			a : array-like
182			1-dimensional array of numeric values
183			idx : integer
184			index of requested element
185
186			Returns
187			-------
188			Index of next non-masked element and next non-masked element
189
190			"""
191			try:
192			next_idx = idx + a[idx:].mask.argmin()
193			if ma.is_masked(a[next_idx]):
194			return None, None
195			else:
196			return next_idx, a[next_idx]
197			except (AttributeError, TypeError, IndexError):
198			return idx, a[idx]
199
200
201			def delete_masked_points(*arrs):
202			"""Delete masked points from arrays.
203
204			Takes arrays and removes masked points to help with calculations and plotting.
205
206			Parameters
207			----------
208			arrs : one or more array-like
209			source arrays
210
211			Returns
212			-------
213			arrs : one or more array-like
214			arrays with masked elements removed
215
216			"""
217			if any(hasattr(a, 'mask') for a in arrs):
218			keep = ~functools.reduce(np.logical_or, (np.ma.getmaskarray(a) for a in arrs))
219			return tuple(ma.asarray(a[keep]) for a in arrs)
220			else:
221			return arrs
222
223
224			@exporter.export
225			def reduce_point_density(points, radius, priority=None):
226			r"""Return a mask to reduce the density of points in irregularly-spaced data.
227
228			This function is used to down-sample a collection of scattered points (e.g. surface
229			data), returning a mask that can be used to select the points from one or more arrays
230			(e.g. arrays of temperature and dew point). The points selected can be controlled by
231			providing an array of ``priority`` values (e.g. rainfall totals to ensure that
232			stations with higher precipitation remain in the mask).
233
234			Parameters
235			----------
236			points : (N, K) array-like
237			N locations of the points in K dimensional space
238			radius : float
239			minimum radius allowed between points
240			priority : (N, K) array-like, optional
241			If given, this should have the same shape as ``points``; these values will
242			be used to control selection priority for points.
243
244			Returns
245			-------
246			(N,) array-like of boolean values indicating whether points should be kept. This
247			can be used directly to index numpy arrays to return only the desired points.
248
249			Examples
250			--------
251			>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.)
252			array([ True, False, True], dtype=bool)
253			>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.,
254			... priority=np.array([0.1, 0.9, 0.3]))
255			array([False, True, False], dtype=bool)
256
257			"""
258			# Handle 1D input
259			if points.ndim < 2:
260			points = points.reshape(-1, 1)
261
262			# Make a kd-tree to speed searching of data.
263			tree = cKDTree(points)
264
265			# Need to use sorted indices rather than sorting the position
266			# so that the keep mask matches original order.
267			if priority is not None:
268			# Need to sort the locations in decreasing priority.
269			sorted_indices = np.argsort(priority)[::-1]
270			else:
271			# Take advantage of iterator nature of range here to avoid making big lists
272			sorted_indices = range(len(points))
273
274			# Keep all points initially
275			keep = np.ones(len(points), dtype=np.bool)
276
277			# Loop over all the potential points
278			for ind in sorted_indices:
279			# Only proceed if we haven't already excluded this point
280			if keep[ind]:
281			# Find the neighbors and eliminate them
282			neighbors = tree.query_ball_point(points[ind], radius)
283			keep[neighbors] = False
284
285			# We just removed ourselves, so undo that
286			keep[ind] = True
287
288			return keep
289
290
291			@exporter.export
292			@units.wraps(None, ('=A', '=A', None))
293			def log_interp(x, xp, args, *kwargs):
294			r"""Interpolates data with logarithmic x-scale over a specified axis.
295
296			Interpolation on a logarithmic x-scale for interpolation values in pressure coordintates.
297
298			Parameters
299			----------
300			x : array-like
301			1-D array of desired interpolated values.
302
303			xp : array-like
304			The x-coordinates of the data points.
305
306			args : array-like
307			The y-coordinates of the data points, same shape as xp.
308
309			axis : int, optional
310			The axis to interpolate over. Defaults to 0.
311
312			fill_value: float, optional
313			Specify handling of interpolation points out of data bounds. If None, will return
314			ValueError if points are out of bounds. Defaults to nan.
315
316			Returns
317			-------
318			array-like
319			Interpolated values for each point with coordinates sorted in ascending order.
320
321			Examples
322			--------
323			>>> x_log = np.array([1e3, 1e4, 1e5, 1e6])
324			>>> y_log = np.log(x_log) * 2 + 3
325			>>> x_interp = np.array([5e3, 5e4, 5e5])
326			>>> metpy.calc.log_interp(x_interp, x_log, y_log)
327			array([ 20.03438638, 24.63955657, 29.24472675])
328
329			Notes
330			-----
331			xp and args must be the same shape.
332
333			"""
334			# Pull out keyword args
335			fill_value = kwargs.pop('fill_value', np.nan)
336			axis = kwargs.pop('axis', 0)
337
338			# Save number of dimensions in xp
339			ndim = xp.ndim
340
341			# Sort input data
342			sort_args = np.argsort(xp, axis=axis)
343			sort_x = np.argsort(x)
344
345			# indices for sorting
346			sorter = broadcast_indices(xp, sort_args, ndim, axis)
347
348			# Ensure pressure in increasing order
349			variables = [arr[sorter] for arr in args]
350
351			# Make x broadcast with xp
352			x_array = x[sort_x]
353			expand = [np.newaxis] * ndim
354			expand[axis] = slice(None)
355			x_array = x_array[expand]
356
357			# Log x and xp
358			log_x_array = np.log(x_array)
359			log_xp = np.log(xp[sorter])
360
361			# Calculate value above interpolated value
362			minv = np.apply_along_axis(np.searchsorted, axis, xp[sorter], x[sort_x])
363			minv2 = np.copy(minv)
364
365			# If extrap is none and data is out of bounds, raise value error
366			if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None:
367			raise ValueError('Interpolation point out of data bounds encountered')
368
369			# Warn if interpolated values are outside data bounds, will make these the values
370			# at end of data range.
371			if np.max(minv) == xp.shape[axis]:
372			warnings.warn('Interpolation point out of data bounds encountered')
373			minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1
374			if np.max(minv) == 0:
375			warnings.warn('Interpolation point out of data bounds encountered')
376			minv2[minv == 0] = 1
377
378			# Get indicies for broadcasting arrays
379			above = broadcast_indices(xp, minv2, ndim, axis)
380			below = broadcast_indices(xp, minv2 - 1, ndim, axis)
381
382			# Create empty output list
383			ret = []
384
385			# Calculate interpolation for each variable
386			for var in variables:
387			var_interp = var[below] + ((log_x_array - log_xp[below]) /
388			(log_xp[above] - log_xp[below])) * (var[above] -
389			var[below])
390
391			# set to nan if fill_value = nan
392			var_interp[minv == xp.shape[axis]] = fill_value
393			var_interp[minv == 0] = fill_value
394
395			# Output to list
396			ret.append(var_interp)
397			if len(ret) == 1:
398			return ret[0]
399			else:
400			return ret
401
402
403			def broadcast_indices(x, minv, ndim, axis):
404			"""Calculate index values to properly broadcast index array within data array.
405
406			<<<<<<< HEAD
407			See usage in log_interp.
408			"""
409			ret = []
410			for dim in range(ndim):
411			if dim == axis:
412			ret.append(minv)
413			else:
414			broadcast_slice = [np.newaxis] * ndim
415			broadcast_slice[dim] = slice(None)
416			dim_inds = np.arange(x.shape[dim])
417			ret.append(dim_inds[broadcast_slice])
418			return ret
419
420
421			def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True):
422			"""Calculate the bounding pressure and height in a layer.
423
424			Given pressure, optional heights, and a bound, return either the closest pressure/height
425			or interpolated pressure/height. If no heights are provided, a standard atmosphere is
426			assumed.
427
428			Parameters
429			----------
430			pressure : `pint.Quantity`
431			Atmospheric pressures
432			bound : `pint.Quantity`
433			Bound to retrieve (in pressure or height)
434			heights : `pint.Quantity`
435			Atmospheric heights associated with the pressure levels
436			interpolate : boolean
437			Interpolate the bound or return the nearest
438
439			Returns
440			-------
441			`pint.Quantity`
442			The bound pressure and height.
443
444			"""
445			# Bound is given in pressure
446			if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
447			# If the bound is in the pressure data, we know the pressure bound exactly
448			if bound in pressure:
449			bound_pressure = bound
450			# If we have heights, we know the exact height value, otherwise return standard
451			# atmosphere height for the pressure
452			if heights is not None:
453			bound_height = heights[pressure == bound_pressure]
454			else:
455			bound_height = pressure_to_height_std(bound_pressure)
456			# If bound is not in the data, return the nearest or interpolated values
457			else:
458			if interpolate:
459			bound_pressure = bound # Use the user specified bound
460			if heights is not None: # Interpolate heights from the height data
461			bound_height = log_interp(bound_pressure, pressure, heights)
462			else: # If not heights given, use the standard atmosphere
463			bound_height = pressure_to_height_std(bound_pressure)
464			else: # No interpolation, find the closest values
465			idx = (np.abs(pressure - bound)).argmin()
466			bound_pressure = pressure[idx]
467			if heights is not None:
468			bound_height = heights[idx]
469			else:
470			bound_height = pressure_to_height_std(bound_pressure)
471
472			# Bound is given in height
473			elif bound.dimensionality == {'[length]': 1.0}:
474			# If there is height data, see if we have the bound or need to interpolate/find nearest
475			if heights is not None:
476			if bound in heights: # Bound is in the height data
477			bound_height = bound
478			bound_pressure = pressure[heights == bound]
479			else: # Bound is not in the data
480			if interpolate:
481			bound_height = bound
482			bound_pressure = height_to_pressure_std(bound_height)
483			else:
484			idx = (np.abs(heights - bound)).argmin()
485			bound_pressure = pressure[idx]
486			bound_height = heights[idx]
487			else: # Don't have heights, so assume a standard atmosphere
488			bound_height = bound
489			bound_pressure = height_to_pressure_std(bound)
490			# If interpolation is on, this is all we need, if not, we need to go back and
491			# find the pressure closest to this and refigure the bounds
492			if not interpolate:
493			idx = (np.abs(pressure - bound_pressure)).argmin()
494			bound_pressure = pressure[idx]
495			bound_height = pressure_to_height_std(bound_pressure)
496
497			# Bound has invalid units
498			else:
499			raise ValueError('Bound must be specified in units of length or pressure.')
500
501			# If the bound is out of the range of the data, we shouldn't extrapolate
502			if (bound_pressure < np.min(pressure)) or (bound_pressure > np.max(pressure)):
503			raise ValueError('Specified bound is outside pressure range.')
504			if heights is not None:
505			if (bound_height > np.max(heights)) or (bound_height < np.min(heights)):
506			raise ValueError('Specified bound is outside height range.')
507
508			return bound_pressure, bound_height
509
510
511			@exporter.export
512			@check_units('[pressure]')
513			def get_layer(p, args, *kwargs):
514			r"""Return an atmospheric layer from upper air data with the requested bottom and depth.
515
516			This function will subset an upper air dataset to contain only the specified layer. The
517			bottom of the layer can be specified with a pressure or height above the surface
518			pressure. The bottom defaults to the surface pressure. The depth of the layer can be
519			specified in terms of pressure or height above the bottom of the layer. If the top and
520			bottom of the layer are not in the data, they are interpolated by default.
521
522			Parameters
523			----------
524			p : array-like
525			Atmospheric pressure profile
526			*args : array-like
527			Atmospheric variable(s) measured at the given pressures
528			heights: array-like
529			Atmospheric heights corresponding to the given pressures
530			bottom : `pint.Quantity`
531			The bottom of the layer as a pressure or height above the surface pressure
532			depth : `pint.Quantity`
533			The thickness of the layer as a pressure or height above the bottom of the layer
534			interpolate : bool
535			Interpolate the top and bottom points if they are not in the given data
536
537			Returns
538			-------
539			`pint.Quantity, pint.Quantity`
540			The pressure and data variables of the layer
541
542			"""
543			# Pop off keyword arguments
544			heights = kwargs.pop('heights', None)
545			bottom = kwargs.pop('bottom', None)
546			depth = kwargs.pop('depth', 100 * units.hPa)
547			interpolate = kwargs.pop('interpolate', True)
548
549			# Make sure pressure and datavars are the same length
550			for datavar in args:
551			if len(p) != len(datavar):
552			raise ValueError('Pressure and data variables must have the same length.')
553
554			# If the bottom is not specified, make it the surface pressure
555			if bottom is None:
556			bottom = p[0]
557
558			bottom_pressure, bottom_height = _get_bound_pressure_height(p, bottom, heights=heights,
559			interpolate=interpolate)
560
561			# Calculate the top if whatever units depth is in
562			if depth.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
563			top = bottom_pressure - depth
564			elif depth.dimensionality == {'[length]': 1}:
565			top = bottom_height + depth
566			else:
567			raise ValueError('Depth must be specified in units of length or pressure')
568
569			top_pressure, _ = _get_bound_pressure_height(p, top, heights=heights,
570			interpolate=interpolate)
571
572			ret = [] # returned data variables in layer
573
574			# Ensure pressures are sorted in ascending order
575			sort_inds = np.argsort(p)
576			p = p[sort_inds]
577
578			# Mask based on top and bottom pressure
579			inds = (p <= bottom_pressure) & (p >= top_pressure)
580			p_interp = p[inds]
581
582			# Interpolate pressures at bounds if necessary and sort
583			if interpolate:
584			# If we don't have the bottom or top requested, append them
585			if top_pressure not in p_interp:
586			p_interp = np.sort(np.append(p_interp, top_pressure)) * p.units
587			if bottom_pressure not in p_interp:
588			p_interp = np.sort(np.append(p_interp, bottom_pressure)) * p.units
589
590			ret.append(p_interp[::-1])
591
592			for datavar in args:
593			# Ensure that things are sorted in ascending order
594			datavar = datavar[sort_inds]
595
596			if interpolate:
597			# Interpolate for the possibly missing bottom/top values
598			datavar_interp = log_interp(p_interp, p, datavar)
599			datavar = datavar_interp
600			else:
601			datavar = datavar[inds]
602
603			ret.append(datavar[::-1])
604
605			return ret
606

Unidata / MetPy

Pull Request — master (#468)

_get_bound_pressure_height() F

Complexity

Size

Duplication

Importance

How to fix Long Method Complexity

Long Method

Complexity

Duplication Side-by-Side

Filter issues like