delete_masked_points() - Code Metrics - Inspection of "Skew T Masked Data" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#378)

unknown

created 2017-04-03 15:36 UTC

delete_masked_points() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	5
dl	0
loc	21
rs	8.439
c	0
b	0
f	0

# 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 numpy as np
import numpy.ma as ma
from scipy.spatial import cKDTree

from ..package_tools import Exporter

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 is 'linear':
        y[nans] = np.interp(x[nans], x[~nans], y[~nans])
    elif kind is '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


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
8			import numpy as np
9			import numpy.ma as ma
10			from scipy.spatial import cKDTree
11
12			from ..package_tools import Exporter
13
14			exporter = Exporter(globals())
15
16
17			@exporter.export
18			def resample_nn_1d(a, centers):
19			"""Return one-dimensional nearest-neighbor indexes based on user-specified centers.
20
21			Parameters
22			----------
23			a : array-like
24			1-dimensional array of numeric values from which to
25			extract indexes of nearest-neighbors
26			centers : array-like
27			1-dimensional array of numeric values representing a subset of values to approximate
28
29			Returns
30			-------
31			An array of indexes representing values closest to given array values
32			"""
33			ix = []
34			for center in centers:
35			index = (np.abs(a - center)).argmin()
36			if index not in ix:
37			ix.append(index)
38			return ix
39
40
41			@exporter.export
42			def nearest_intersection_idx(a, b):
43			"""Determine the index of the point just before two lines with common x values.
44
45			Parameters
46			----------
47			a : array-like
48			1-dimensional array of y-values for line 1
49			b : array-like
50			1-dimensional array of y-values for line 2
51
52			Returns
53			-------
54			An array of indexes representing the index of the values
55			just before the intersection(s) of the two lines.
56			"""
57			# Difference in the two y-value sets
58			difference = a - b
59
60			# Determine the point just before the intersection of the lines
61			# Will return multiple points for multiple intersections
62			sign_change_idx, = np.nonzero(np.diff(np.sign(difference)))
63
64			return sign_change_idx
65
66
67			@exporter.export
68			def find_intersections(x, a, b, direction='all'):
69			"""Calculate the best estimate of intersection.
70
71			Calculates the best estimates of the intersection of two y-value
72			data sets that share a common x-value set.
73
74			Parameters
75			----------
76			x : array-like
77			1-dimensional array of numeric x-values
78			a : array-like
79			1-dimensional array of y-values for line 1
80			b : array-like
81			1-dimensional array of y-values for line 2
82			direction : string
83			specifies direction of crossing. 'all', 'increasing' (a becoming greater than b),
84			or 'decreasing' (b becoming greater than a).
85
86			Returns
87			-------
88			A tuple (x, y) of array-like with the x and y coordinates of the
89			intersections of the lines.
90			"""
91			# Find the index of the points just before the intersection(s)
92			nearest_idx = nearest_intersection_idx(a, b)
93			next_idx = nearest_idx + 1
94
95			# Determine the sign of the change
96			sign_change = np.sign(a[next_idx] - b[next_idx])
97
98			# x-values around each intersection
99			_, x0 = _next_non_masked_element(x, nearest_idx)
100			_, x1 = _next_non_masked_element(x, next_idx)
101
102			# y-values around each intersection for the first line
103			_, a0 = _next_non_masked_element(a, nearest_idx)
104			_, a1 = _next_non_masked_element(a, next_idx)
105
106			# y-values around each intersection for the second line
107			_, b0 = _next_non_masked_element(b, nearest_idx)
108			_, b1 = _next_non_masked_element(b, next_idx)
109
110			# Calculate the x-intersection. This comes from finding the equations of the two lines,
111			# one through (x0, a0) and (x1, a1) and the other through (x0, b0) and (x1, b1),
112			# finding their intersection, and reducing with a bunch of algebra.
113			delta_y0 = a0 - b0
114			delta_y1 = a1 - b1
115			intersect_x = (delta_y1 * x0 - delta_y0 * x1) / (delta_y1 - delta_y0)
116
117			# Calculate the y-intersection of the lines. Just plug the x above into the equation
118			# for the line through the a points. One could solve for y like x above, but this
119			# causes weirder unit behavior and seems a little less good numerically.
120			intersect_y = ((intersect_x - x0) / (x1 - x0)) * (a1 - a0) + a0
121
122			# Make a mask based on the direction of sign change desired
123			if direction == 'increasing':
124			mask = sign_change > 0
125			elif direction == 'decreasing':
126			mask = sign_change < 0
127			elif direction == 'all':
128			return intersect_x, intersect_y
129			else:
130			raise ValueError('Unknown option for direction: {0}'.format(str(direction)))
131			return intersect_x[mask], intersect_y[mask]
132
133
134			@exporter.export
135			def interpolate_nans(x, y, kind='linear'):
136			"""Interpolate NaN values in y.
137
138			Interpolate NaN values in the y dimension. Works with unsorted x values.
139
140			Parameters
141			----------
142			x : array-like
143			1-dimensional array of numeric x-values
144			y : array-like
145			1-dimensional array of numeric y-values
146			kind : string
147			specifies the kind of interpolation x coordinate - 'linear' or 'log'
148
149			Returns
150			-------
151			An array of the y coordinate data with NaN values interpolated.
152			"""
153			x_sort_args = np.argsort(x)
154			x = x[x_sort_args]
155			y = y[x_sort_args]
156			nans = np.isnan(y)
157			if kind is 'linear':
158			y[nans] = np.interp(x[nans], x[~nans], y[~nans])
159			elif kind is 'log':
160			y[nans] = np.interp(np.log(x[nans]), np.log(x[~nans]), y[~nans])
161			else:
162			raise ValueError('Unknown option for kind: {0}'.format(str(kind)))
163			return y[x_sort_args]
164
165
166			def _next_non_masked_element(a, idx):
167			"""Return the next non masked element of a masked array.
168
169			If an array is masked, return the next non-masked element (if the given index is masked).
170			If no other unmasked points are after the given masked point, returns none.
171
172			Parameters
173			----------
174			a : array-like
175			1-dimensional array of numeric values
176			idx : integer
177			index of requested element
178
179			Returns
180			-------
181			Index of next non-masked element and next non-masked element
182			"""
183			try:
184			next_idx = idx + a[idx:].mask.argmin()
185			if ma.is_masked(a[next_idx]):
186			return None, None
187			else:
188			return next_idx, a[next_idx]
189			except (AttributeError, TypeError, IndexError):
190			return idx, a[idx]
191
192
193			def delete_masked_points(*arrs):
194			"""Delete masked points from arrays.
195
196			Takes arrays and removes masked points to help with calculations and plotting.
197
198			Parameters
199			----------
200			arrs : one or more array-like
201			source arrays
202
203			Returns
204			-------
205			arrs : one or more array-like
206			arrays with masked elements removed
207
208			"""
209			if any(hasattr(a, 'mask') for a in arrs):
210			keep = ~functools.reduce(np.logical_or, (np.ma.getmaskarray(a) for a in arrs))
211			return tuple(ma.asarray(a[keep]) for a in arrs)
212			else:
213			return arrs
214
215
216			@exporter.export
217			def reduce_point_density(points, radius, priority=None):
218			r"""Return a mask to reduce the density of points in irregularly-spaced data.
219
220			This function is used to down-sample a collection of scattered points (e.g. surface
221			data), returning a mask that can be used to select the points from one or more arrays
222			(e.g. arrays of temperature and dew point). The points selected can be controlled by
223			providing an array of ``priority`` values (e.g. rainfall totals to ensure that
224			stations with higher precipitation remain in the mask).
225
226			Parameters
227			----------
228			points : (N, K) array-like
229			N locations of the points in K dimensional space
230			radius : float
231			minimum radius allowed between points
232			priority : (N, K) array-like, optional
233			If given, this should have the same shape as ``points``; these values will
234			be used to control selection priority for points.
235
236			Returns
237			-------
238			(N,) array-like of boolean values indicating whether points should be kept. This
239			can be used directly to index numpy arrays to return only the desired points.
240
241			Examples
242			--------
243			>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.)
244			array([ True, False, True], dtype=bool)
245			>>> metpy.calc.reduce_point_density(np.array([1, 2, 3]), 1.,
246			priority=np.array([0.1, 0.9, 0.3]))
247			array([False, True, False], dtype=bool)
248			"""
249			# Handle 1D input
250			if points.ndim < 2:
251			points = points.reshape(-1, 1)
252
253			# Make a kd-tree to speed searching of data.
254			tree = cKDTree(points)
255
256			# Need to use sorted indices rather than sorting the position
257			# so that the keep mask matches original order.
258			if priority is not None:
259			# Need to sort the locations in decreasing priority.
260			sorted_indices = np.argsort(priority)[::-1]
261			else:
262			# Take advantage of iterator nature of range here to avoid making big lists
263			sorted_indices = range(len(points))
264
265			# Keep all points initially
266			keep = np.ones(len(points), dtype=np.bool)
267
268			# Loop over all the potential points
269			for ind in sorted_indices:
270			# Only proceed if we haven't already excluded this point
271			if keep[ind]:
272			# Find the neighbors and eliminate them
273			neighbors = tree.query_ball_point(points[ind], radius)
274			keep[neighbors] = False
275
276			# We just removed ourselves, so undo that
277			keep[ind] = True
278
279			return keep
280

Unidata / MetPy

Pull Request — master (#378)

delete_masked_points() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like