log_interpolation() - Code Metrics - Inspection of "Layer helpers" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#444)

unknown

created 2017-06-12 14:56 UTC

log_interpolation() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	1
c	1
b	0
f	0
dl	0
loc	23
rs	9.0856

# 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 ..calc import (dewpoint, potential_temperature, pressure_at_height_above_pressure,
                   saturation_mixing_ratio, vapor_pressure)
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 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


@exporter.export
@check_units('[pressure]')
def get_layer(p, datavar, bottom=None, depth=100 * units.hPa, interpolate=True):
    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
    datavar : array-like
        Atmospheric variable measured at 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 variable of the layer

    """
    # Make sure pressure and datavar are the same length
    if len(p) != len(datavar):
        raise ValueError('Pressure and data variable must have the same length.')

    # If no bottom is specified, use the first pressure value
    if not bottom:
        bottom_pressure = p[0]

    # If bottom is specified as a height (AGL), convert to pressure
    if bottom:
        if bottom.dimensionality == {'[length]': 1.0}:
            bottom_pressure = pressure_at_height_above_pressure(p[0], bottom)
        if bottom.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
            bottom_pressure = bottom

    # Calculate the pressure at the top of the layer
    if depth.dimensionality == {'[length]': 1.0}:
        top_pressure = pressure_at_height_above_pressure(bottom_pressure, depth)
    else:
        top_pressure = bottom_pressure - depth

    # Handle top or bottom values that are invalid
    if bottom_pressure > p[0]:
        raise ValueError(
            'Bottom of layer pressure is greater than the largest pressure in data.')
    if top_pressure < p[-1]:
        raise ValueError('Top of layer pressure is less than the lowest pressure in data.')
    if bottom_pressure < top_pressure:
        raise ValueError(
            'Pressure at the top of the layer is greater than that at the bottom.')

        # Mask the pressure values we have data for in the layer
    inds = (p <= bottom_pressure) & (p >= top_pressure)
    p_interp = p[inds]

    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)) * units.hPa
        if bottom_pressure not in p_interp:
            p_interp = np.sort(np.append(p_interp, bottom_pressure)) * units.hPa

        # Interpolate for the possibly missing bottom/top pressure values
        sort_args = np.argsort(p)
        datavar_interp = np.interp(p_interp, p[sort_args], datavar[sort_args]) * units.degC
        p = p_interp
        datavar = datavar_interp
    else:
        datavar = datavar[inds]
        p = p_interp

    return p, datavar


@exporter.export
@check_units('[pressure]')
def mixed_layer(p, datavar, bottom=None, depth=100 * units.hPa, interpolate=True):
    r"""Mix a variable over a layer using pressure as a mass-average.

    This function will integrate a data variable with respect to pressure and using the
    mean value theorem determine the average value.

    Parameters
    ----------
    p : array-like
        Atmospheric pressure profile
    datavar : array-like
        Atmospheric variable measured at 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`
        The mixed value of the data variable.

    """
    p_layer, datavar_layer = get_layer(p, datavar, bottom, depth, interpolate)
    actual_depth = abs(p_layer[0] - p_layer[-1])
    return (1. / actual_depth.m) * np.trapz(datavar_layer, p_layer) * datavar.units


@exporter.export
@check_units('[pressure]', '[temperature]', '[temperature]')
def mixed_parcel(p, T, Td, parcel_start_pressure=None, bottom=None, depth=100 * units.hPa,
                 interpolate=True):
    r"""Calculate the properties of a parcel mixed from a layer.

    Determines the properties of an air parcel that is the result of complete mixing of a
    given atmospheric layer.

    Parameters
    ----------
    p : array-like
        Atmospheric pressure profile
    T : array-like
        Atmospheric temperature profile
    Td : array-like
        Atmospheric dewpoint profile
    parcel_start_pressure : `pint.Quantity`
        Pressure at which the mixed parcel should begin
    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, pint.Quantity`
        The pressure, temperature, and dewpoint of the mixed parcel.

    """
    # If a parcel starting pressure is not provided, use the surface
    if not parcel_start_pressure:
        parcel_start_pressure = p[0]

    # Calculate the potential temperature and mixing ratio over the layer
    theta = potential_temperature(p, T)
    mixing_ratio = saturation_mixing_ratio(p, Td)

    # Mix the variables over the layer
    mean_theta = mixed_layer(p, theta, bottom, depth, interpolate)
    mean_mixing_ratio = mixed_layer(p, mixing_ratio, bottom, depth, interpolate)

    # Convert back to temperature
    mean_temperature = mean_theta / potential_temperature(parcel_start_pressure,
                                                                 1 * units.degK)

    # Convert back to dewpoint
    mean_vapor_pressure = vapor_pressure(parcel_start_pressure, mean_mixing_ratio)
    mean_dewpoint = dewpoint(mean_vapor_pressure)

    return parcel_start_pressure, mean_temperature.to(T.units), mean_dewpoint.to(Td.units)


def log_interpolation(x, xp, fp, **kwargs):
    """Interpolates data with logarithmic x-scale.

    Interpolation on a logarithmic x-scale for interpolating values in pressure coordinates.

    Parameters
    ----------
    x : array-like
        The x-coordinates of the interpolated values.
    xp : array-like
        The x-coordinates of the data points.
    fp : array-like
        The y-coordinates of the data points, same length as xp.

    Returns
    -------
    array-like
        The interpolated values, same shape as x.

    """
    sort_args = np.argsort(xp)
    interpolated_vals = np.interp(np.log(x.m), np.log(xp.m[sort_args]), fp[sort_args], **kwargs)
    return interpolated_vals * fp.units

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 ..calc import (dewpoint, potential_temperature, pressure_at_height_above_pressure,
13			saturation_mixing_ratio, vapor_pressure)
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 is 'linear':
165			y[nans] = np.interp(x[nans], x[~nans], y[~nans])
166			elif kind is '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			@check_units('[pressure]')
293			def get_layer(p, datavar, bottom=None, depth=100 * units.hPa, interpolate=True):
294			r"""Return an atmospheric layer from upper air data with the requested bottom and depth.
295
296			This function will subset an upper air dataset to contain only the specified layer. The
297			bottom of the layer can be specified with a pressure or height above the surface
298			pressure. The bottom defaults to the surface pressure. The depth of the layer can be
299			specified in terms of pressure or height above the bottom of the layer. If the top and
300			bottom of the layer are not in the data, they are interpolated by default.
301
302			Parameters
303			----------
304			p : array-like
305			Atmospheric pressure profile
306			datavar : array-like
307			Atmospheric variable measured at the given pressures
308			bottom : `pint.Quantity`
309			The bottom of the layer as a pressure or height above the surface pressure
310			depth : `pint.Quantity`
311			The thickness of the layer as a pressure or height above the bottom of the layer
312			interpolate : bool
313			Interpolate the top and bottom points if they are not in the given data
314
315			Returns
316			-------
317			`pint.Quantity, pint.Quantity`
318			The pressure and data variable of the layer
319
320			"""
321			# Make sure pressure and datavar are the same length
322			if len(p) != len(datavar):
323			raise ValueError('Pressure and data variable must have the same length.')
324
325			# If no bottom is specified, use the first pressure value
326			if not bottom:
327			bottom_pressure = p[0]
328
329			# If bottom is specified as a height (AGL), convert to pressure
330			if bottom:
331			if bottom.dimensionality == {'[length]': 1.0}:
332			bottom_pressure = pressure_at_height_above_pressure(p[0], bottom)
333			if bottom.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
334			bottom_pressure = bottom
335
336			# Calculate the pressure at the top of the layer
337			if depth.dimensionality == {'[length]': 1.0}:
338			top_pressure = pressure_at_height_above_pressure(bottom_pressure, depth)
339			else:
340			top_pressure = bottom_pressure - depth
341
342			# Handle top or bottom values that are invalid
343			if bottom_pressure > p[0]:
344			raise ValueError(
345			'Bottom of layer pressure is greater than the largest pressure in data.')
346			if top_pressure < p[-1]:
347			raise ValueError('Top of layer pressure is less than the lowest pressure in data.')
348			if bottom_pressure < top_pressure:
349			raise ValueError(
350			'Pressure at the top of the layer is greater than that at the bottom.')
351
352			# Mask the pressure values we have data for in the layer
353			inds = (p <= bottom_pressure) & (p >= top_pressure)
354			p_interp = p[inds]
355
356			if interpolate:
357			# If we don't have the bottom or top requested, append them
358			if top_pressure not in p_interp:
359			p_interp = np.sort(np.append(p_interp, top_pressure)) * units.hPa
360			if bottom_pressure not in p_interp:
361			p_interp = np.sort(np.append(p_interp, bottom_pressure)) * units.hPa
362
363			# Interpolate for the possibly missing bottom/top pressure values
364			sort_args = np.argsort(p)
365			datavar_interp = np.interp(p_interp, p[sort_args], datavar[sort_args]) * units.degC
366			p = p_interp
367			datavar = datavar_interp
368			else:
369			datavar = datavar[inds]
370			p = p_interp
371
372			return p, datavar
373
374
375			@exporter.export
376			@check_units('[pressure]')
377			def mixed_layer(p, datavar, bottom=None, depth=100 * units.hPa, interpolate=True):
378			r"""Mix a variable over a layer using pressure as a mass-average.
379
380			This function will integrate a data variable with respect to pressure and using the
381			mean value theorem determine the average value.
382
383			Parameters
384			----------
385			p : array-like
386			Atmospheric pressure profile
387			datavar : array-like
388			Atmospheric variable measured at the given pressures
389			bottom : `pint.Quantity`
390			The bottom of the layer as a pressure or height above the surface pressure
391			depth : `pint.Quantity`
392			The thickness of the layer as a pressure or height above the bottom of the layer
393			interpolate : bool
394			Interpolate the top and bottom points if they are not in the given data
395
396			Returns
397			-------
398			`pint.Quantity`
399			The mixed value of the data variable.
400
401			"""
402			p_layer, datavar_layer = get_layer(p, datavar, bottom, depth, interpolate)
403			actual_depth = abs(p_layer[0] - p_layer[-1])
404			return (1. / actual_depth.m) * np.trapz(datavar_layer, p_layer) * datavar.units
405
406
407			@exporter.export
408			@check_units('[pressure]', '[temperature]', '[temperature]')
409			def mixed_parcel(p, T, Td, parcel_start_pressure=None, bottom=None, depth=100 * units.hPa,
410			interpolate=True):
411			r"""Calculate the properties of a parcel mixed from a layer.
412
413			Determines the properties of an air parcel that is the result of complete mixing of a
414			given atmospheric layer.
415
416			Parameters
417			----------
418			p : array-like
419			Atmospheric pressure profile
420			T : array-like
421			Atmospheric temperature profile
422			Td : array-like
423			Atmospheric dewpoint profile
424			parcel_start_pressure : `pint.Quantity`
425			Pressure at which the mixed parcel should begin
426			bottom : `pint.Quantity`
427			The bottom of the layer as a pressure or height above the surface pressure
428			depth : `pint.Quantity`
429			The thickness of the layer as a pressure or height above the bottom of the layer
430			interpolate : bool
431			Interpolate the top and bottom points if they are not in the given data
432
433			Returns
434			-------
435			`pint.Quantity, pint.Quantity, pint.Quantity`
436			The pressure, temperature, and dewpoint of the mixed parcel.
437
438			"""
439			# If a parcel starting pressure is not provided, use the surface
440			if not parcel_start_pressure:
441			parcel_start_pressure = p[0]
442
443			# Calculate the potential temperature and mixing ratio over the layer
444			theta = potential_temperature(p, T)
445			mixing_ratio = saturation_mixing_ratio(p, Td)
446
447			# Mix the variables over the layer
448			mean_theta = mixed_layer(p, theta, bottom, depth, interpolate)
449			mean_mixing_ratio = mixed_layer(p, mixing_ratio, bottom, depth, interpolate)
450
451			# Convert back to temperature
452			mean_temperature = mean_theta / potential_temperature(parcel_start_pressure,
453			1 * units.degK)
454
455			# Convert back to dewpoint
456			mean_vapor_pressure = vapor_pressure(parcel_start_pressure, mean_mixing_ratio)
457			mean_dewpoint = dewpoint(mean_vapor_pressure)
458
459			return parcel_start_pressure, mean_temperature.to(T.units), mean_dewpoint.to(Td.units)
460
461
462			def log_interpolation(x, xp, fp, **kwargs):
463			"""Interpolates data with logarithmic x-scale.
464
465			Interpolation on a logarithmic x-scale for interpolating values in pressure coordinates.
466
467			Parameters
468			----------
469			x : array-like
470			The x-coordinates of the interpolated values.
471			xp : array-like
472			The x-coordinates of the data points.
473			fp : array-like
474			The y-coordinates of the data points, same length as xp.
475
476			Returns
477			-------
478			array-like
479			The interpolated values, same shape as x.
480
481			"""
482			sort_args = np.argsort(xp)
483			interpolated_vals = np.interp(np.log(x.m), np.log(xp.m[sort_args]), fp[sort_args], **kwargs)
484			return interpolated_vals * fp.units

Unidata / MetPy

Pull Request — master (#444)

log_interpolation() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like