shearing_deformation() - Code Metrics - Inspection of "Adding shearing and stretching deformation calcs" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#402)

unknown

created 2017-04-24 14:47 UTC

shearing_deformation() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
dl	0
loc	29
rs	8.8571
c	0
b	0
f	0

# Copyright (c) 2008-2015 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains calculation of kinematic parameters (e.g. divergence or vorticity)."""
from __future__ import division

import functools
import warnings

import numpy as np

from ..cbook import is_string_like, iterable
from ..constants import g
from ..package_tools import Exporter
from ..units import atleast_2d, concatenate, units

exporter = Exporter(globals())


def _gradient(f, *args, **kwargs):
    """Wrapper for :func:`numpy.gradient` that handles units."""
    if len(args) < f.ndim:
        args = list(args)
        args.extend([units.Quantity(1., 'dimensionless')] * (f.ndim - len(args)))
    grad = np.gradient(f, *(a.magnitude for a in args), **kwargs)
    if f.ndim == 1:
        return units.Quantity(grad, f.units / args[0].units)
    return [units.Quantity(g, f.units / dx.units) for dx, g in zip(args, grad)]


def _stack(arrs):
    return concatenate([a[np.newaxis] for a in arrs], axis=0)


def _get_gradients(u, v, dx, dy):
    """Helper function for getting convergence and vorticity from 2D arrays."""
    dudy, dudx = _gradient(u, dy, dx)
    dvdy, dvdx = _gradient(v, dy, dx)
    return dudx, dudy, dvdx, dvdy


def _is_x_first_dim(dim_order):
    """Determine whether x is the first dimension based on the value of dim_order."""
    if dim_order is None:
        warnings.warn('dim_order is using the default setting (currently "xy"). This will '
                      'change to "yx" in the next version. It is recommended that you '
                      'specify the appropriate ordering ("xy", "yx") for your data by '
                      'passing the `dim_order` argument to the calculation.', FutureWarning)
        dim_order = 'xy'
    return dim_order == 'xy'


def _check_and_flip(arr):
    """Transpose array or list of arrays if they are 2D."""
    if hasattr(arr, 'ndim'):
        if arr.ndim >= 2:
            return arr.T
        else:
            return arr
    elif not is_string_like(arr) and iterable(arr):
        return tuple(_check_and_flip(a) for a in arr)
    else:
        return arr


def ensure_yx_order(func):
    """Wrap a function to ensure all array arguments are y, x ordered, based on kwarg."""
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        # Check what order we're given
        dim_order = kwargs.pop('dim_order', None)
        x_first = _is_x_first_dim(dim_order)

        # If x is the first dimension, flip (transpose) every array within the function args.
        if x_first:
            args = tuple(_check_and_flip(arr) for arr in args)
            for k, v in kwargs:
                kwargs[k] = _check_and_flip(v)

        ret = func(*args, **kwargs)

        # If we flipped on the way in, need to flip on the way out so that output array(s)
        # match the dimension order of the original input.
        if x_first:
            return _check_and_flip(ret)
        else:
            return ret

    # Inject a docstring for the dim_order argument into the function's docstring.
    dim_order_doc = '''
    dim_order : str or ``None``, optional
        The ordering of dimensions in passed in arrays. Can be one of ``None``, ``'xy'``,
        or ``'yx'``. ``'xy'`` indicates that the dimension corresponding to x is the leading
        dimension, followed by y. ``'yx'`` indicates that x is the last dimension, preceded
        by y. ``None`` indicates that the default ordering should be assumed,
        which will change in version 0.6 from 'xy' to 'yx'. Can only be passed as a keyword
        argument, i.e. func(..., dim_order='xy').'''

    # Find the first blank line after the start of the parameters section
    params = wrapper.__doc__.find('Parameters')
    blank = wrapper.__doc__.find('\n\n', params)
    wrapper.__doc__ = wrapper.__doc__[:blank] + dim_order_doc + wrapper.__doc__[blank:]

    return wrapper


@exporter.export
@ensure_yx_order
def v_vorticity(u, v, dx, dy):
    r"""Calculate the vertical vorticity of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    (M, N) ndarray
        vertical vorticity

    See Also
    --------
    h_convergence, convergence_vorticity
    """
    _, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
    return dvdx - dudy


@exporter.export
@ensure_yx_order
def h_convergence(u, v, dx, dy):
    r"""Calculate the horizontal convergence of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    (M, N) ndarray
        The horizontal convergence

    See Also
    --------
    v_vorticity, convergence_vorticity
    """
    dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
    return dudx + dvdy


@exporter.export
@ensure_yx_order
def convergence_vorticity(u, v, dx, dy):
    r"""Calculate the horizontal convergence and vertical vorticity of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    convergence, vorticity : tuple of (M, N) ndarrays
        The horizontal convergence and vertical vorticity, respectively

    See Also
    --------
    v_vorticity, h_convergence

    Notes
    -----
    This is a convenience function that will do less work than calculating
    the horizontal convergence and vertical vorticity separately.
    """
    dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
    return dudx + dvdy, dvdx - dudy


@exporter.export
@ensure_yx_order
def shearing_deformation(u, v, dx, dy):
    r"""Calculate the shearing deformation of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    (M, N) ndarray
        Shearing Deformation

    See Also
    --------
    stretching_convergence, shearing_stretching_deformation
    """
    _, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
    return dvdx + dudy


@exporter.export
@ensure_yx_order
def stretching_deformation(u, v, dx, dy):
    r"""Calculate the stretching deformation of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    (M, N) ndarray
        Stretching Deformation

    See Also
    --------
    shearing_deformation, shearing_stretching_deformation
    """
    dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
    return dudx - dvdy


@exporter.export
@ensure_yx_order
def shearing_stretching_deformation(u, v, dx, dy):
    r"""Calculate the horizontal shearing and stretching deformation of the horizontal wind.

    The grid must have a constant spacing in each direction.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind
    v : (M, N) ndarray
        y component of the wind
    dx : float
        The grid spacing in the x-direction
    dy : float
        The grid spacing in the y-direction

    Returns
    -------
    shearing, strectching : tuple of (M, N) ndarrays
        The horizontal shearing and stretching deformation, respectively

    See Also
    --------
    shearing_deformation, stretching_deformation

    Notes
    -----
    This is a convenience function that will do less work than calculating
    the shearing and streching deformation terms separately.
    """
    dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
    return dvdx + dudy, dudx - dvdy


@exporter.export
@ensure_yx_order
def advection(scalar, wind, deltas):
    r"""Calculate the advection of a scalar field by the wind.

    The order of the dimensions of the arrays must match the order in which
    the wind components are given.  For example, if the winds are given [u, v],
    then the scalar and wind arrays must be indexed as x,y (which puts x as the
    rows, not columns).

    Parameters
    ----------
    scalar : N-dimensional array
        Array (with N-dimensions) with the quantity to be advected.
    wind : sequence of arrays
        Length N sequence of N-dimensional arrays.  Represents the flow,
        with a component of the wind in each dimension.  For example, for
        horizontal advection, this could be a list: [u, v], where u and v
        are each a 2-dimensional array.
    deltas : sequence
        A (length N) sequence containing the grid spacing in each dimension.

    Returns
    -------
    N-dimensional array
        An N-dimensional array containing the advection at all grid points.
    """
    # This allows passing in a list of wind components or an array.
    wind = _stack(wind)

    # If we have more than one component, we need to reverse the order along the first
    # dimension so that the wind components line up with the
    # order of the gradients from the ..., y, x ordered array.
    if wind.ndim > scalar.ndim:
        wind = wind[::-1]

    # Gradient returns a list of derivatives along each dimension. We convert
    # this to an array with dimension as the first index. Reverse the deltas to line up
    # with the order of the dimensions.
    grad = _stack(_gradient(scalar, *deltas[::-1]))

    # Make them be at least 2D (handling the 1D case) so that we can do the
    # multiply and sum below
    grad, wind = atleast_2d(grad, wind)

    return (-grad * wind).sum(axis=0)


@exporter.export
@ensure_yx_order
def geostrophic_wind(heights, f, dx, dy):
    r"""Calculate the geostrophic wind given from the heights or geopotential.

    Parameters
    ----------
    heights : (M, N) ndarray
        The height field, with either leading dimensions of (x, y) or trailing dimensions
        of (y, x), depending on the value of ``dim_order``.
    f : array_like
        The coriolis parameter.  This can be a scalar to be applied
        everywhere or an array of values.
    dx : scalar
        The grid spacing in the x-direction
    dy : scalar
        The grid spacing in the y-direction

    Returns
    -------
    A 2-item tuple of arrays
        A tuple of the u-component and v-component of the geostrophic wind.
    """
    if heights.dimensionality['[length]'] == 2.0:
        norm_factor = 1. / f
    else:
        norm_factor = g / f

    # If heights has more than 2 dimensions, we need to pass in some dummy
    # grid deltas so that we can still use np.gradient. It may be better to
    # to loop in this case, but that remains to be done.
    deltas = [dy, dx]
    if heights.ndim > 2:
        deltas = [units.Quantity(1., units.m)] * (heights.ndim - 2) + deltas

    grad = _gradient(heights, *deltas)
    dy, dx = grad[-2:]  # Throw away unused gradient components
    return -norm_factor * dy, norm_factor * dx


1			# Copyright (c) 2008-2015 MetPy Developers.
2			# Distributed under the terms of the BSD 3-Clause License.
3			# SPDX-License-Identifier: BSD-3-Clause
4			"""Contains calculation of kinematic parameters (e.g. divergence or vorticity)."""
5			from __future__ import division
6
7			import functools
8			import warnings
9
10			import numpy as np
11
12			from ..cbook import is_string_like, iterable
13			from ..constants import g
14			from ..package_tools import Exporter
15			from ..units import atleast_2d, concatenate, units
16
17			exporter = Exporter(globals())
18
19
20			def _gradient(f, args, *kwargs):
21			"""Wrapper for :func:`numpy.gradient` that handles units."""
22			if len(args) < f.ndim:
23			args = list(args)
24			args.extend([units.Quantity(1., 'dimensionless')] * (f.ndim - len(args)))
25			grad = np.gradient(f, (a.magnitude for a in args), *kwargs)
26			if f.ndim == 1:
27			return units.Quantity(grad, f.units / args[0].units)
28			return [units.Quantity(g, f.units / dx.units) for dx, g in zip(args, grad)]
29
30
31			def _stack(arrs):
32			return concatenate([a[np.newaxis] for a in arrs], axis=0)
33
34
35			def _get_gradients(u, v, dx, dy):
36			"""Helper function for getting convergence and vorticity from 2D arrays."""
37			dudy, dudx = _gradient(u, dy, dx)
38			dvdy, dvdx = _gradient(v, dy, dx)
39			return dudx, dudy, dvdx, dvdy
40
41
42			def _is_x_first_dim(dim_order):
43			"""Determine whether x is the first dimension based on the value of dim_order."""
44			if dim_order is None:
45			warnings.warn('dim_order is using the default setting (currently "xy"). This will '
46			'change to "yx" in the next version. It is recommended that you '
47			'specify the appropriate ordering ("xy", "yx") for your data by '
48			'passing the `dim_order` argument to the calculation.', FutureWarning)
49			dim_order = 'xy'
50			return dim_order == 'xy'
51
52
53			def _check_and_flip(arr):
54			"""Transpose array or list of arrays if they are 2D."""
55			if hasattr(arr, 'ndim'):
56			if arr.ndim >= 2:
57			return arr.T
58			else:
59			return arr
60			elif not is_string_like(arr) and iterable(arr):
61			return tuple(_check_and_flip(a) for a in arr)
62			else:
63			return arr
64
65
66			def ensure_yx_order(func):
67			"""Wrap a function to ensure all array arguments are y, x ordered, based on kwarg."""
68			@functools.wraps(func)
69			def wrapper(args, *kwargs):
70			# Check what order we're given
71			dim_order = kwargs.pop('dim_order', None)
72			x_first = _is_x_first_dim(dim_order)
73
74			# If x is the first dimension, flip (transpose) every array within the function args.
75			if x_first:
76			args = tuple(_check_and_flip(arr) for arr in args)
77			for k, v in kwargs:
78			kwargs[k] = _check_and_flip(v)
79
80			ret = func(args, *kwargs)
81
82			# If we flipped on the way in, need to flip on the way out so that output array(s)
83			# match the dimension order of the original input.
84			if x_first:
85			return _check_and_flip(ret)
86			else:
87			return ret
88
89			# Inject a docstring for the dim_order argument into the function's docstring.
90			dim_order_doc = '''
91			dim_order : str or ``None``, optional
92			The ordering of dimensions in passed in arrays. Can be one of ``None``, ``'xy'``,
93			or ``'yx'``. ``'xy'`` indicates that the dimension corresponding to x is the leading
94			dimension, followed by y. ``'yx'`` indicates that x is the last dimension, preceded
95			by y. ``None`` indicates that the default ordering should be assumed,
96			which will change in version 0.6 from 'xy' to 'yx'. Can only be passed as a keyword
97			argument, i.e. func(..., dim_order='xy').'''
98
99			# Find the first blank line after the start of the parameters section
100			params = wrapper.__doc__.find('Parameters')
101			blank = wrapper.__doc__.find('\n\n', params)
102			wrapper.__doc__ = wrapper.__doc__[:blank] + dim_order_doc + wrapper.__doc__[blank:]
103
104			return wrapper
105
106
107			@exporter.export
108			@ensure_yx_order
109			def v_vorticity(u, v, dx, dy):
110			r"""Calculate the vertical vorticity of the horizontal wind.
111
112			The grid must have a constant spacing in each direction.
113
114			Parameters
115			----------
116			u : (M, N) ndarray
117			x component of the wind
118			v : (M, N) ndarray
119			y component of the wind
120			dx : float
121			The grid spacing in the x-direction
122			dy : float
123			The grid spacing in the y-direction
124
125			Returns
126			-------
127			(M, N) ndarray
128			vertical vorticity
129
130			See Also
131			--------
132			h_convergence, convergence_vorticity
133			"""
134			_, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
135			return dvdx - dudy
136
137
138			@exporter.export
139			@ensure_yx_order
140			def h_convergence(u, v, dx, dy):
141			r"""Calculate the horizontal convergence of the horizontal wind.
142
143			The grid must have a constant spacing in each direction.
144
145			Parameters
146			----------
147			u : (M, N) ndarray
148			x component of the wind
149			v : (M, N) ndarray
150			y component of the wind
151			dx : float
152			The grid spacing in the x-direction
153			dy : float
154			The grid spacing in the y-direction
155
156			Returns
157			-------
158			(M, N) ndarray
159			The horizontal convergence
160
161			See Also
162			--------
163			v_vorticity, convergence_vorticity
164			"""
165			dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
166			return dudx + dvdy
167
168
169			@exporter.export
170			@ensure_yx_order
171			def convergence_vorticity(u, v, dx, dy):
172			r"""Calculate the horizontal convergence and vertical vorticity of the horizontal wind.
173
174			The grid must have a constant spacing in each direction.
175
176			Parameters
177			----------
178			u : (M, N) ndarray
179			x component of the wind
180			v : (M, N) ndarray
181			y component of the wind
182			dx : float
183			The grid spacing in the x-direction
184			dy : float
185			The grid spacing in the y-direction
186
187			Returns
188			-------
189			convergence, vorticity : tuple of (M, N) ndarrays
190			The horizontal convergence and vertical vorticity, respectively
191
192			See Also
193			--------
194			v_vorticity, h_convergence
195
196			Notes
197			-----
198			This is a convenience function that will do less work than calculating
199			the horizontal convergence and vertical vorticity separately.
200			"""
201			dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
202			return dudx + dvdy, dvdx - dudy
203
204
205			@exporter.export
206			@ensure_yx_order
207			def shearing_deformation(u, v, dx, dy):
208			r"""Calculate the shearing deformation of the horizontal wind.
209
210			The grid must have a constant spacing in each direction.
211
212			Parameters
213			----------
214			u : (M, N) ndarray
215			x component of the wind
216			v : (M, N) ndarray
217			y component of the wind
218			dx : float
219			The grid spacing in the x-direction
220			dy : float
221			The grid spacing in the y-direction
222
223			Returns
224			-------
225			(M, N) ndarray
226			Shearing Deformation
227
228			See Also
229			--------
230			stretching_convergence, shearing_stretching_deformation
231			"""
232			_, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
233			return dvdx + dudy
234
235
236			@exporter.export
237			@ensure_yx_order
238			def stretching_deformation(u, v, dx, dy):
239			r"""Calculate the stretching deformation of the horizontal wind.
240
241			The grid must have a constant spacing in each direction.
242
243			Parameters
244			----------
245			u : (M, N) ndarray
246			x component of the wind
247			v : (M, N) ndarray
248			y component of the wind
249			dx : float
250			The grid spacing in the x-direction
251			dy : float
252			The grid spacing in the y-direction
253
254			Returns
255			-------
256			(M, N) ndarray
257			Stretching Deformation
258
259			See Also
260			--------
261			shearing_deformation, shearing_stretching_deformation
262			"""
263			dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
264			return dudx - dvdy
265
266
267			@exporter.export
268			@ensure_yx_order
269			def shearing_stretching_deformation(u, v, dx, dy):
270			r"""Calculate the horizontal shearing and stretching deformation of the horizontal wind.
271
272			The grid must have a constant spacing in each direction.
273
274			Parameters
275			----------
276			u : (M, N) ndarray
277			x component of the wind
278			v : (M, N) ndarray
279			y component of the wind
280			dx : float
281			The grid spacing in the x-direction
282			dy : float
283			The grid spacing in the y-direction
284
285			Returns
286			-------
287			shearing, strectching : tuple of (M, N) ndarrays
288			The horizontal shearing and stretching deformation, respectively
289
290			See Also
291			--------
292			shearing_deformation, stretching_deformation
293
294			Notes
295			-----
296			This is a convenience function that will do less work than calculating
297			the shearing and streching deformation terms separately.
298			"""
299			dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
300			return dvdx + dudy, dudx - dvdy
301
302
303			@exporter.export
304			@ensure_yx_order
305			def advection(scalar, wind, deltas):
306			r"""Calculate the advection of a scalar field by the wind.
307
308			The order of the dimensions of the arrays must match the order in which
309			the wind components are given. For example, if the winds are given [u, v],
310			then the scalar and wind arrays must be indexed as x,y (which puts x as the
311			rows, not columns).
312
313			Parameters
314			----------
315			scalar : N-dimensional array
316			Array (with N-dimensions) with the quantity to be advected.
317			wind : sequence of arrays
318			Length N sequence of N-dimensional arrays. Represents the flow,
319			with a component of the wind in each dimension. For example, for
320			horizontal advection, this could be a list: [u, v], where u and v
321			are each a 2-dimensional array.
322			deltas : sequence
323			A (length N) sequence containing the grid spacing in each dimension.
324
325			Returns
326			-------
327			N-dimensional array
328			An N-dimensional array containing the advection at all grid points.
329			"""
330			# This allows passing in a list of wind components or an array.
331			wind = _stack(wind)
332
333			# If we have more than one component, we need to reverse the order along the first
334			# dimension so that the wind components line up with the
335			# order of the gradients from the ..., y, x ordered array.
336			if wind.ndim > scalar.ndim:
337			wind = wind[::-1]
338
339			# Gradient returns a list of derivatives along each dimension. We convert
340			# this to an array with dimension as the first index. Reverse the deltas to line up
341			# with the order of the dimensions.
342			grad = _stack(_gradient(scalar, *deltas[::-1]))
343
344			# Make them be at least 2D (handling the 1D case) so that we can do the
345			# multiply and sum below
346			grad, wind = atleast_2d(grad, wind)
347
348			return (-grad * wind).sum(axis=0)
349
350
351			@exporter.export
352			@ensure_yx_order
353			def geostrophic_wind(heights, f, dx, dy):
354			r"""Calculate the geostrophic wind given from the heights or geopotential.
355
356			Parameters
357			----------
358			heights : (M, N) ndarray
359			The height field, with either leading dimensions of (x, y) or trailing dimensions
360			of (y, x), depending on the value of ``dim_order``.
361			f : array_like
362			The coriolis parameter. This can be a scalar to be applied
363			everywhere or an array of values.
364			dx : scalar
365			The grid spacing in the x-direction
366			dy : scalar
367			The grid spacing in the y-direction
368
369			Returns
370			-------
371			A 2-item tuple of arrays
372			A tuple of the u-component and v-component of the geostrophic wind.
373			"""
374			if heights.dimensionality['[length]'] == 2.0:
375			norm_factor = 1. / f
376			else:
377			norm_factor = g / f
378
379			# If heights has more than 2 dimensions, we need to pass in some dummy
380			# grid deltas so that we can still use np.gradient. It may be better to
381			# to loop in this case, but that remains to be done.
382			deltas = [dy, dx]
383			if heights.ndim > 2:
384			deltas = [units.Quantity(1., units.m)] * (heights.ndim - 2) + deltas
385
386			grad = _gradient(heights, *deltas)
387			dy, dx = grad[-2:] # Throw away unused gradient components
388			return -norm_factor * dy, norm_factor * dx
389

Unidata / MetPy

Pull Request — master (#402)

shearing_deformation() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like