Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add gradient and derivatives for irregular grids #692

Merged
merged 13 commits into from
Jan 4, 2018
Merged

Add gradient and derivatives for irregular grids #692

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
b4e09a2
MNT: Remove hyphen from reference tag
dopplershift Dec 29, 2017
6990115
MNT: Update matplotlib import in conftest
dopplershift Dec 29, 2017
fdac87e
ENH: Make unit versions of at_least1d/atleast_2d handle arrays
dopplershift Jan 2, 2018
43806e7
ENH: Add functions for derivatives on irregular grids
dopplershift Dec 30, 2017
b36f782
BUG: Fix incorrect deprecation
dopplershift Jan 2, 2018
09ca99d
MNT: Elminate divergence_vorticity
dopplershift Jan 2, 2018
8bf6020
MNT: Deprecate shearing_stretching_deformation
dopplershift Jan 2, 2018
fade954
MNT: Silence deprecation warnings in tests
dopplershift Jan 3, 2018
be27c60
MNT: Fix typos
dopplershift Jan 3, 2018
35d1824
ENH: Replace use of numpy.gradient wrapper
dopplershift Jan 3, 2018
2e32e2c
MNT: Disable printing test log output by default
dopplershift Jan 3, 2018
379a8bf
ENH: Promote diff wrapper to part of the API
dopplershift Jan 3, 2018
f73d4b8
MNT: Handle normalize_axis_index for old numpy
dopplershift Jan 3, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions conftest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
"""Configure pytest for metpy."""

import matplotlib
import matplotlib.pyplot
import numpy
import pint
import pytest
Expand Down
6 changes: 5 additions & 1 deletion docs/references.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@ References
doi:`10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2
<https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2>`_.

.. [Bowen2005] Bowen, M. K. and R. Smith, 2005: Derivative formulae and errors for
non-uniformly spaced points. *Proc. R. Soc. A*, **461**, 1975-1997,
doi:`10.1098/rspa.2004.1430 <https://doi.org/10.1098/rspa.2004.1430>`_.

.. [Bunkers2000] Bunkers, M. J., B. A. Klimowski, J. W. Zeitler, R. L. Thompson,
and M. L. Weisman, 2000: Predicting supercell motion using a new hodograph
technique. *Wea. Forecasting.*, **15**, 61–79,
Expand All @@ -31,7 +35,7 @@ References
doi:`10.1175/1520-0493(1959)087%3C0367:AOOAS%3E2.0.CO;2
<https://doi.org/10.1175/1520-0493(1959)087%3C0367:AOOAS%3E2.0.CO;2>`_.

.. [Davies-Jones2009] Davies-Jones, R., 2009: On Formulas for Equivalent Potential Temperature.
.. [DaviesJones2009] Davies-Jones, R., 2009: On Formulas for Equivalent Potential Temperature.
*Mon. Wea. Rev.*, **137**, 3137-3148,
doi: `10.1175/2009mwr2774.1 <https://doi.org/10.1175/2009MWR2774.1>`_.

Expand Down
154 changes: 60 additions & 94 deletions metpy/calc/kinematics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
import numpy as np
from pyproj import Geod

from ..calc.tools import get_layer_heights
from ..calc.tools import first_derivative, get_layer_heights, gradient
from ..cbook import is_string_like, iterable
from ..constants import Cp_d, g
from ..deprecation import deprecated
Expand All @@ -20,28 +20,10 @@
exporter = Exporter(globals())


def _gradient(f, *args, **kwargs):
"""Wrap :func:`numpy.gradient` to handle 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):
"""Return derivatives for components to simplify calculating divergence and vorticity."""
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:
Expand Down Expand Up @@ -112,8 +94,6 @@ def wrapper(*args, **kwargs):
def 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
Expand All @@ -132,10 +112,11 @@ def vorticity(u, v, dx, dy):

See Also
--------
divergence, divergence_vorticity
divergence

"""
_, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
dudy = first_derivative(u, delta=dy, axis=0)
dvdx = first_derivative(v, delta=dx, axis=1)
return dvdx - dudy


Expand All @@ -147,16 +128,16 @@ def v_vorticity(u, v, dx, dy, dim_order='xy'):
return vorticity(u, v, dx, dy, dim_order=dim_order)


v_vorticity.__doc__ = vorticity.__doc__
v_vorticity.__doc__ = (vorticity.__doc__ +
'\n .. deprecated:: 0.7.0\n Function has been renamed to '
'`vorticity` and will be removed from MetPy in 0.9.0.')


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

The grid must have a constant spacing in each direction.

Parameters
----------
u : (M, N) ndarray
Expand All @@ -175,31 +156,34 @@ def divergence(u, v, dx, dy):

See Also
--------
vorticity, divergence_vorticity
vorticity

"""
dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
dudx = first_derivative(u, delta=dx, axis=1)
dvdy = first_derivative(v, delta=dy, axis=0)
return dudx + dvdy


@exporter.export
@deprecated('0.7', addendum=' This function has been replaced by divergence.',
pending=False)
def convergence(u, v, dx, dy, dim_order='xy'):
def h_convergence(u, v, dx, dy, dim_order='xy'):
"""Wrap divergence for deprecated convergence function."""
return divergence(u, v, dx, dy, dim_order=dim_order)


convergence.__doc__ = divergence.__doc__
h_convergence.__doc__ = (divergence.__doc__ +
'\n .. deprecated:: 0.7.0\n Function has been renamed to '
'`divergence` and will be removed from MetPy in 0.9.0.')


@exporter.export
@deprecated('0.7', addendum=' Use divergence and/or vorticity instead.',
pending=False)
@ensure_yx_order
def divergence_vorticity(u, v, dx, dy):
def convergence_vorticity(u, v, dx, dy, dim_order='xy'):
r"""Calculate the horizontal divergence and vertical vorticity of the horizontal wind.

The grid must have a constant spacing in each direction.

Parameters
----------
u : (M, N) ndarray
Expand All @@ -220,34 +204,24 @@ def divergence_vorticity(u, v, dx, dy):
--------
vorticity, divergence

Notes
-----
This is a convenience function that will do less work than calculating
the horizontal divergence and vertical vorticity separately.
.. deprecated:: 0.7.0
Function no longer has any performance benefit over individual calls to
`divergence` and `vorticity` and will be removed from MetPy in 0.9.0.


"""
dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
dudx = first_derivative(u, delta=dx, axis=1)
dudy = first_derivative(u, delta=dy, axis=0)
dvdx = first_derivative(v, delta=dx, axis=1)
dvdy = first_derivative(v, delta=dy, axis=0)
return dudx + dvdy, dvdx - dudy


@exporter.export
@deprecated('0.7', addendum=' This function has been replaced by divergence_vorticity.',
pending=False)
def convergence_vorticity(u, v, dx, dy, dim_order='xy'):
"""Wrap divergence_vorticity for deprecated convergence vorticity function."""
return divergence_vorticity(u, v, dx, dy, dim_order=dim_order)


convergence_vorticity.__doc__ = divergence_vorticity.__doc__


@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
Expand All @@ -266,10 +240,11 @@ def shearing_deformation(u, v, dx, dy):

See Also
--------
stretching_deformation, shearing_stretching_deformation
stretching_deformation, total_deformation

"""
_, dudy, dvdx, _ = _get_gradients(u, v, dx, dy)
dudy = first_derivative(u, delta=dy, axis=0)
dvdx = first_derivative(v, delta=dx, axis=1)
return dvdx + dudy


Expand All @@ -278,8 +253,6 @@ def shearing_deformation(u, v, dx, dy):
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
Expand All @@ -298,20 +271,21 @@ def stretching_deformation(u, v, dx, dy):

See Also
--------
shearing_deformation, shearing_stretching_deformation
shearing_deformation, total_deformation

"""
dudx, _, _, dvdy = _get_gradients(u, v, dx, dy)
dudx = first_derivative(u, delta=dx, axis=1)
dvdy = first_derivative(v, delta=dy, axis=0)
return dudx - dvdy


@exporter.export
@deprecated('0.7', addendum=' Use stretching_deformation and/or shearing_deformation instead.',
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should also probably note in the docs like we do elsewhere:

    .. deprecated:: 0.6.0
      Function has been moved to the Siphon library and will be removed from MetPy in 0.8.0.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I missed that we did that. Will fix.

pending=False)
@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
Expand All @@ -332,13 +306,17 @@ def shearing_stretching_deformation(u, v, dx, dy):
--------
shearing_deformation, stretching_deformation

Notes
-----
This is a convenience function that will do less work than calculating
the shearing and streching deformation terms separately.

.. deprecated:: 0.7.0
Function no longer has any performance benefit over individual calls to
`shearing_deformation` and `stretching_deformation` and will be removed from
MetPy in 0.9.0.

"""
dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
dudx = first_derivative(u, delta=dx, axis=1)
dudy = first_derivative(u, delta=dy, axis=0)
dvdx = first_derivative(v, delta=dx, axis=1)
dvdy = first_derivative(v, delta=dy, axis=0)
return dvdx + dudy, dudx - dvdy


Expand All @@ -347,8 +325,6 @@ def shearing_stretching_deformation(u, v, dx, dy):
def total_deformation(u, v, dx, dy):
r"""Calculate the horizontal total deformation of the horizontal wind.

The grid must have a constant spacing in each direction.

Parameters
----------
u : (M, N) ndarray
Expand All @@ -367,16 +343,13 @@ def total_deformation(u, v, dx, dy):

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

Notes
-----
This is a convenience function that will do less work than calculating
the shearing and streching deformation terms separately and calculating the
total deformation "by hand".
shearing_deformation, stretching_deformation

"""
dudx, dudy, dvdx, dvdy = _get_gradients(u, v, dx, dy)
dudx = first_derivative(u, delta=dx, axis=1)
dudy = first_derivative(u, delta=dy, axis=0)
dvdx = first_derivative(v, delta=dx, axis=1)
dvdy = first_derivative(v, delta=dy, axis=0)
return np.sqrt((dvdx + dudy)**2 + (dudx - dvdy)**2)


Expand Down Expand Up @@ -420,7 +393,7 @@ def advection(scalar, wind, deltas):
# 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]))
grad = _stack(gradient(scalar, deltas=deltas[::-1]))

# Make them be at least 2D (handling the 1D case) so that we can do the
# multiply and sum below
Expand All @@ -445,10 +418,6 @@ def frontogenesis(thta, u, v, dx, dy, dim_order='yx'):
* :math:`\beta` is the angle between the axis of dilitation and the isentropes
* :math:`\delta` is the divergence

Notes
-----
Assumes dim_order='yx', unless otherwise specified.

Parameters
----------
thta : (M, N) ndarray
Expand All @@ -465,22 +434,26 @@ def frontogenesis(thta, u, v, dx, dy, dim_order='yx'):
Returns
-------
(M, N) ndarray
2D Frotogenesis in [temperature units]/m/s
2D Frontogenesis in [temperature units]/m/s

Notes
-----
Assumes dim_order='yx', unless otherwise specified.

Conversion factor to go from [temperature units]/m/s to [tempature units/100km/3h]
Conversion factor to go from [temperature units]/m/s to [temperature units/100km/3h]
:math:`1.08e4*1.e5`

"""
# Get gradients of potential temperature in both x and y
grad = _gradient(thta, dy, dx)
ddy_thta, ddx_thta = grad[-2:] # Throw away unused gradient components
ddy_thta = first_derivative(thta, delta=dy, axis=-2)
ddx_thta = first_derivative(thta, delta=dx, axis=-1)

# Compute the magnitude of the potential temperature gradient
mag_thta = np.sqrt(ddx_thta**2 + ddy_thta**2)

# Get the shearing, stretching, and total deformation of the wind field
shrd, strd = shearing_stretching_deformation(u, v, dx, dy, dim_order=dim_order)
shrd = shearing_deformation(u, v, dx, dy, dim_order=dim_order)
strd = stretching_deformation(u, v, dx, dy, dim_order=dim_order)
tdef = total_deformation(u, v, dx, dy, dim_order=dim_order)

# Get the divergence of the wind field
Expand Down Expand Up @@ -522,16 +495,9 @@ def geostrophic_wind(heights, f, dx, dy):
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
dhdy = first_derivative(heights, delta=dy, axis=-2)
dhdx = first_derivative(heights, delta=dx, axis=-1)
return -norm_factor * dhdy, norm_factor * dhdx


@exporter.export
Expand Down
Loading