diff --git a/conftest.py b/conftest.py
index d140f83a2e2..dd2ae12e831 100644
--- a/conftest.py
+++ b/conftest.py
@@ -4,6 +4,7 @@
 """Configure pytest for metpy."""
 
 import matplotlib
+import matplotlib.pyplot
 import numpy
 import pint
 import pytest
diff --git a/docs/references.rst b/docs/references.rst
index 3dcab24cf9a..5b32d9495a3 100644
--- a/docs/references.rst
+++ b/docs/references.rst
@@ -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,
@@ -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>`_.
 
diff --git a/metpy/calc/kinematics.py b/metpy/calc/kinematics.py
index 9fdc123a6f9..02c17ae709c 100644
--- a/metpy/calc/kinematics.py
+++ b/metpy/calc/kinematics.py
@@ -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
@@ -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:
@@ -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
@@ -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
 
 
@@ -147,7 +128,9 @@ 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
@@ -155,8 +138,6 @@ def v_vorticity(u, v, dx, dy, dim_order='xy'):
 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
@@ -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
@@ -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
@@ -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
 
 
@@ -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
@@ -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.',
+            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
@@ -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
 
 
@@ -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
@@ -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)
 
 
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/metpy/calc/tests/test_kinematics.py b/metpy/calc/tests/test_kinematics.py
index 617ca4f46b4..e198b79721f 100644
--- a/metpy/calc/tests/test_kinematics.py
+++ b/metpy/calc/tests/test_kinematics.py
@@ -6,13 +6,14 @@
 import numpy as np
 import pytest
 
-from metpy.calc import (advection, convergence, convergence_vorticity, divergence,
-                        divergence_vorticity, frontogenesis, geostrophic_wind,
-                        get_wind_components, lat_lon_grid_spacing, montgomery_streamfunction,
+from metpy.calc import (advection, convergence_vorticity, divergence,
+                        frontogenesis, geostrophic_wind, get_wind_components, h_convergence,
+                        lat_lon_grid_spacing, montgomery_streamfunction,
                         shearing_deformation, shearing_stretching_deformation,
                         storm_relative_helicity, stretching_deformation, total_deformation,
                         v_vorticity, vorticity)
 from metpy.constants import g, omega, Re
+from metpy.deprecation import MetpyDeprecationWarning
 from metpy.testing import assert_almost_equal, assert_array_equal
 from metpy.units import concatenate, units
 
@@ -21,13 +22,14 @@ def test_default_order_warns():
     """Test that using the default array ordering issues a warning."""
     u = np.ones((3, 3)) * units('m/s')
     with pytest.warns(FutureWarning):
-        divergence_vorticity(u, u, 1 * units.meter, 1 * units.meter)
+        vorticity(u, u, 1 * units.meter, 1 * units.meter)
 
 
 def test_zero_gradient():
     """Test divergence_vorticity when there is no gradient in the field."""
     u = np.ones((3, 3)) * units('m/s')
-    c, v = divergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c, v = convergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
     truth = np.zeros_like(u) / units.sec
     assert_array_equal(c, truth)
     assert_array_equal(v, truth)
@@ -37,7 +39,8 @@ def test_cv_zero_vorticity():
     """Test divergence_vorticity when there is only divergence."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    c, v = divergence_vorticity(u, u.T, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c, v = convergence_vorticity(u, u.T, 1 * units.meter, 1 * units.meter, dim_order='xy')
     true_c = 2. * np.ones_like(u) / units.sec
     true_v = np.zeros_like(u) / units.sec
     assert_array_equal(c, true_c)
@@ -48,7 +51,8 @@ def test_divergence_vorticity():
     """Test of vorticity and divergence calculation for basic case."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    c, v = divergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c, v = convergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
     true_c = np.ones_like(u) / units.sec
     true_v = np.ones_like(u) / units.sec
     assert_array_equal(c, true_c)
@@ -59,15 +63,18 @@ def test_vorticity_divergence_asym():
     """Test vorticity and divergence calculation with a complicated field."""
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
-    c, vort = divergence_vorticity(u, v, 1 * units.meters, 2 * units.meters, dim_order='yx')
-    true_c = np.array([[0., 4., 0.], [1., 0.5, -0.5], [2., 0., 5.]]) / units.sec
-    true_vort = np.array([[-1., 2., 7.], [3.5, -1.5, -6.], [-2., 0., 1.]]) / units.sec
+    with pytest.warns(MetpyDeprecationWarning):
+        c, vort = convergence_vorticity(u, v, 1 * units.meters, 2 * units.meters,
+                                        dim_order='yx')
+    true_c = np.array([[-2, 5.5, -2.5], [2., 0.5, -1.5], [3., -1.5, 8.5]]) / units.sec
+    true_vort = np.array([[-2.5, 3.5, 13.], [8.5, -1.5, -11.], [-5.5, -1.5, 0.]]) / units.sec
     assert_array_equal(c, true_c)
     assert_array_equal(vort, true_vort)
 
     # Now try for xy ordered
-    c, vort = divergence_vorticity(u.T, v.T, 1 * units.meters, 2 * units.meters,
-                                   dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c, vort = convergence_vorticity(u.T, v.T, 1 * units.meters, 2 * units.meters,
+                                        dim_order='xy')
     assert_array_equal(c, true_c.T)
     assert_array_equal(vort, true_vort.T)
 
@@ -95,7 +102,7 @@ def test_vorticity_asym():
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
     vort = vorticity(u, v, 1 * units.meters, 2 * units.meters, dim_order='yx')
-    true_vort = np.array([[-1., 2., 7.], [3.5, -1.5, -6.], [-2., 0., 1.]]) / units.sec
+    true_vort = np.array([[-2.5, 3.5, 13.], [8.5, -1.5, -11.], [-5.5, -1.5, 0.]]) / units.sec
     assert_array_equal(vort, true_vort)
 
     # Now try for xy ordered
@@ -126,7 +133,7 @@ def test_divergence_asym():
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
     c = divergence(u, v, 1 * units.meters, 2 * units.meters, dim_order='yx')
-    true_c = np.array([[0., 4., 0.], [1., 0.5, -0.5], [2., 0., 5.]]) / units.sec
+    true_c = np.array([[-2, 5.5, -2.5], [2., 0.5, -1.5], [3., -1.5, 8.5]]) / units.sec
     assert_array_equal(c, true_c)
 
     # Now try for xy ordered
@@ -137,8 +144,9 @@ def test_divergence_asym():
 def test_shst_zero_gradient():
     """Test shear_stretching_deformation when there is zero gradient."""
     u = np.ones((3, 3)) * units('m/s')
-    sh, st = shearing_stretching_deformation(u, u, 1 * units.meter, 1 * units.meter,
-                                             dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        sh, st = shearing_stretching_deformation(u, u, 1 * units.meter, 1 * units.meter,
+                                                 dim_order='xy')
     truth = np.zeros_like(u) / units.sec
     assert_array_equal(sh, truth)
     assert_array_equal(st, truth)
@@ -148,8 +156,9 @@ def test_shst_zero_stretching():
     """Test shear_stretching_deformation when there is only shearing."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    sh, st = shearing_stretching_deformation(u, u.T, 1 * units.meter, 1 * units.meter,
-                                             dim_order='yx')
+    with pytest.warns(MetpyDeprecationWarning):
+        sh, st = shearing_stretching_deformation(u, u.T, 1 * units.meter, 1 * units.meter,
+                                                 dim_order='yx')
     true_sh = 2. * np.ones_like(u) / units.sec
     true_st = np.zeros_like(u) / units.sec
     assert_array_equal(sh, true_sh)
@@ -160,8 +169,9 @@ def test_shst_deformation():
     """Test of shearing and stretching deformation calculation for basic case."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    sh, st = shearing_stretching_deformation(u, u, 1 * units.meter, 1 * units.meter,
-                                             dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        sh, st = shearing_stretching_deformation(u, u, 1 * units.meter, 1 * units.meter,
+                                                 dim_order='xy')
     true_sh = np.ones_like(u) / units.sec
     true_st = np.ones_like(u) / units.sec
     assert_array_equal(sh, true_st)
@@ -172,16 +182,18 @@ def test_shst_deformation_asym():
     """Test shearing and stretching deformation calculation with a complicated field."""
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
-    sh, st = shearing_stretching_deformation(u, v, 1 * units.meters, 2 * units.meters,
-                                             dim_order='yx')
-    true_sh = np.array([[-3., 0., 1.], [4.5, -0.5, -6.], [2., 4., 7.]]) / units.sec
-    true_st = np.array([[4., 2., 8.], [3., 1.5, 0.5], [2., 4., -1.]]) / units.sec
+    with pytest.warns(MetpyDeprecationWarning):
+        sh, st = shearing_stretching_deformation(u, v, 1 * units.meters, 2 * units.meters,
+                                                 dim_order='yx')
+    true_sh = np.array([[-7.5, -1.5, 1.], [9.5, -0.5, -11.], [1.5, 5.5, 12.]]) / units.sec
+    true_st = np.array([[4., 0.5, 12.5], [4., 1.5, -0.5], [1., 5.5, -4.5]]) / units.sec
     assert_array_equal(sh, true_sh)
     assert_array_equal(st, true_st)
 
     # Now try for yx ordered
-    sh, st = shearing_stretching_deformation(u.T, v.T, 1 * units.meters, 2 * units.meters,
-                                             dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        sh, st = shearing_stretching_deformation(u.T, v.T, 1 * units.meters, 2 * units.meters,
+                                                 dim_order='xy')
     assert_array_equal(sh, true_sh.T)
     assert_array_equal(st, true_st.T)
 
@@ -191,7 +203,7 @@ def test_shearing_deformation_asym():
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
     sh = shearing_deformation(u, v, 1 * units.meters, 2 * units.meters, dim_order='yx')
-    true_sh = np.array([[-3., 0., 1.], [4.5, -0.5, -6.], [2., 4., 7.]]) / units.sec
+    true_sh = np.array([[-7.5, -1.5, 1.], [9.5, -0.5, -11.], [1.5, 5.5, 12.]]) / units.sec
     assert_array_equal(sh, true_sh)
 
     # Now try for yx ordered
@@ -205,7 +217,7 @@ def test_stretching_deformation_asym():
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
     st = stretching_deformation(u, v, 1 * units.meters, 2 * units.meters, dim_order='yx')
-    true_st = np.array([[4., 2., 8.], [3., 1.5, 0.5], [2., 4., -1.]]) / units.sec
+    true_st = np.array([[4., 0.5, 12.5], [4., 1.5, -0.5], [1., 5.5, -4.5]]) / units.sec
     assert_array_equal(st, true_st)
 
     # Now try for yx ordered
@@ -220,8 +232,8 @@ def test_total_deformation_asym():
     v = np.array([[6, 4, 8], [2, 6, 0], [2, 2, 6]]) * units('m/s')
     tdef = total_deformation(u, v, 1 * units.meters, 2 * units.meters,
                              dim_order='yx')
-    true_tdef = np.array([[5., 2., 8.06225775], [5.40832691, 1.58113883, 6.02079729],
-                          [2.82842712, 5.65685425, 7.07106781]]) / units.sec
+    true_tdef = np.array([[8.5, 1.58113883, 12.5399362], [10.30776406, 1.58113883, 11.0113578],
+                          [1.80277562, 7.7781746, 12.8160056]]) / units.sec
     assert_almost_equal(tdef, true_tdef)
 
     # Now try for xy ordered
@@ -237,9 +249,9 @@ def test_frontogenesis_asym():
     theta = np.array([[303, 295, 305], [308, 310, 312], [299, 293, 289]]) * units('K')
     fronto = frontogenesis(theta, u, v, 1 * units.meters, 2 * units.meters,
                            dim_order='yx')
-    true_fronto = np.array([[-20.93890452, -7.83070042, -36.43293256],
-                            [0.89442719, -2.12218672, -8.94427191],
-                            [-16.8, -7.65600391, -61.65921479]]
+    true_fronto = np.array([[-52.4746386, -37.3658646, -50.3996939],
+                            [3.5777088, -2.1221867, -16.9941166],
+                            [-23.1417334, 26.0499143, -158.4839684]]
                            ) * units.K / units.meter / units.sec
     assert_almost_equal(fronto, true_fronto)
 
@@ -291,7 +303,7 @@ def test_advection_2d():
     v = 2 * np.ones((3, 3)) * units('m/s')
     s = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) * units.kelvin
     a = advection(s, [u, v], (1 * units.meter, 1 * units.meter), dim_order='xy')
-    truth = np.array([[-3, -2, 1], [-4, 0, 4], [-1, 2, 3]]) * units('K/sec')
+    truth = np.array([[-6, -4, 2], [-8, 0, 8], [-2, 4, 6]]) * units('K/sec')
     assert_array_equal(a, truth)
 
 
@@ -301,7 +313,7 @@ def test_advection_2d_asym():
     v = 2 * u
     s = np.array([[1, 2, 4], [4, 8, 4], [8, 6, 4]]) * units.kelvin
     a = advection(s, [u, v], (2 * units.meter, 1 * units.meter), dim_order='yx')
-    truth = np.array([[0, -12.75, -2], [-27., -16., 10.], [-42, 35, 8]]) * units('K/sec')
+    truth = np.array([[0, -20.75, -2.5], [-33., -16., 20.], [-48, 91., 8]]) * units('K/sec')
     assert_array_equal(a, truth)
 
     # Now try xy ordered
@@ -315,7 +327,7 @@ def test_geostrophic_wind():
     # Using g as the value for f allows it to cancel out
     ug, vg = geostrophic_wind(z, g.magnitude / units.sec,
                               100. * units.meter, 100. * units.meter, dim_order='xy')
-    true_u = np.array([[-1, 0, 1]] * 3) * units('m/s')
+    true_u = np.array([[-2, 0, 2]] * 3) * units('m/s')
     true_v = -true_u.T
     assert_array_equal(ug, true_u)
     assert_array_equal(vg, true_v)
@@ -327,8 +339,8 @@ def test_geostrophic_wind_asym():
     # Using g as the value for f allows it to cancel out
     ug, vg = geostrophic_wind(z, g.magnitude / units.sec,
                               200. * units.meter, 100. * units.meter, dim_order='yx')
-    true_u = -np.array([[6, 12, 0], [7, 4, 0], [8, -4, 0]]) * units('m/s')
-    true_v = np.array([[1, 1.5, 2], [4, 0, -4], [-2, -2, -2]]) * units('m/s')
+    true_u = -np.array([[5, 20, 0], [7, 4, 0], [9, -12, 0]]) * units('m/s')
+    true_v = np.array([[0.5, 1.5, 2.5], [8, 0, -8], [-2, -2, -2]]) * units('m/s')
     assert_array_equal(ug, true_u)
     assert_array_equal(vg, true_v)
 
@@ -344,7 +356,7 @@ def test_geostrophic_geopotential():
     z = np.array([[48, 49, 48], [49, 50, 49], [48, 49, 48]]) * 100. * units('m^2/s^2')
     ug, vg = geostrophic_wind(z, 1 / units.sec, 100. * units.meter, 100. * units.meter,
                               dim_order='xy')
-    true_u = np.array([[-1, 0, 1]] * 3) * units('m/s')
+    true_u = np.array([[-2, 0, 2]] * 3) * units('m/s')
     true_v = -true_u.T
     assert_array_equal(ug, true_u)
     assert_array_equal(vg, true_v)
@@ -357,7 +369,7 @@ def test_geostrophic_3d():
     z3d = np.dstack((z, z)) * units.meter
     ug, vg = geostrophic_wind(z3d, g.magnitude / units.sec,
                               100. * units.meter, 100. * units.meter, dim_order='xy')
-    true_u = np.array([[-1, 0, 1]] * 3) * units('m/s')
+    true_u = np.array([[-2, 0, 2]] * 3) * units('m/s')
     true_v = -true_u.T
 
     true_u = concatenate((true_u[..., None], true_u[..., None]), axis=2)
@@ -521,7 +533,8 @@ def test_v_vorticity():
     """Test that v_vorticity wrapper works (deprecated in 0.7)."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    v = v_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        v = v_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
     true_v = np.ones_like(u) / units.sec
     assert_array_equal(v, true_v)
 
@@ -530,7 +543,8 @@ def test_convergence():
     """Test that convergence wrapper works (deprecated in 0.7)."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    c = convergence(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c = h_convergence(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
     true_c = np.ones_like(u) / units.sec
     assert_array_equal(c, true_c)
 
@@ -539,7 +553,8 @@ def test_convergence_vorticity():
     """Test that convergence_vorticity wrapper works (deprecated in 0.7)."""
     a = np.arange(3)
     u = np.c_[a, a, a] * units('m/s')
-    c, v = convergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
+    with pytest.warns(MetpyDeprecationWarning):
+        c, v = convergence_vorticity(u, u, 1 * units.meter, 1 * units.meter, dim_order='xy')
     true_c = np.ones_like(u) / units.sec
     true_v = np.ones_like(u) / units.sec
     assert_array_equal(c, true_c)
diff --git a/metpy/calc/tests/test_tools.py b/metpy/calc/tests/test_tools.py
index b3cc22a4f0e..be227c6577a 100644
--- a/metpy/calc/tests/test_tools.py
+++ b/metpy/calc/tests/test_tools.py
@@ -3,14 +3,16 @@
 # SPDX-License-Identifier: BSD-3-Clause
 """Test the `tools` module."""
 
+from collections import namedtuple
+
 import numpy as np
 import numpy.ma as ma
 import pytest
 
-from metpy.calc import (find_intersections, get_layer, get_layer_heights,
-                        interp, interpolate_nans, log_interp,
+from metpy.calc import (find_intersections, first_derivative, get_layer, get_layer_heights,
+                        interp, interpolate_nans, laplacian, log_interp,
                         nearest_intersection_idx, pressure_to_height_std,
-                        reduce_point_density, resample_nn_1d)
+                        reduce_point_density, resample_nn_1d, second_derivative)
 from metpy.calc.tools import (_get_bound_pressure_height, _greater_or_close, _less_or_close,
                               _next_non_masked_element, delete_masked_points)
 from metpy.testing import assert_array_almost_equal, assert_array_equal
@@ -526,3 +528,103 @@ def test_get_layer_heights_agl_bottom_no_interp():
     data_true = np.array([50, 60, 70, 80, 90, 100]) * units.degC
     assert_array_almost_equal(heights_true, heights, 6)
     assert_array_almost_equal(data_true, data, 6)
+
+
+@pytest.fixture()
+def deriv_1d_data():
+    """Return 1-dimensional data for testing derivative functions."""
+    return namedtuple('D_1D_Test_Data', 'x values')(np.array([0, 1.25, 3.75]) * units.cm,
+                                                    np.array([13.5, 12, 10]) * units.degC)
+
+
+@pytest.fixture()
+def deriv_2d_data():
+    """Return 2-dimensional data for analytic function for testing derivative functions."""
+    ret = namedtuple('D_2D_Test_Data', 'x y x0 y0 a b f')(
+        np.array([0., 2., 7.]), np.array([1., 5., 11., 13.]), 3, 1.5, 0.5, 0.25, 0)
+
+    # Makes a value array with y changing along rows (axis 0) and x along columns (axis 1)
+    return ret._replace(f=ret.a * (ret.x - ret.x0)**2 + ret.b * (ret.y[:, None] - ret.y0)**2)
+
+
+def test_first_derivative(deriv_1d_data):
+    """Test first_derivative with a simple 1D array."""
+    dv_dx = first_derivative(deriv_1d_data.values, x=deriv_1d_data.x)
+
+    # Worked by hand and taken from Chapra and Canale 23.2
+    truth = np.array([-1.333333, -1.06666667, -0.5333333]) * units('delta_degC / cm')
+    assert_array_almost_equal(dv_dx, truth, 5)
+
+
+def test_first_derivative_2d(deriv_2d_data):
+    """Test first_derivative with a full 2D array."""
+    df_dx = first_derivative(deriv_2d_data.f, x=deriv_2d_data.x, axis=1)
+    df_dx_analytic = np.tile(2 * deriv_2d_data.a * (deriv_2d_data.x - deriv_2d_data.x0),
+                             (deriv_2d_data.f.shape[0], 1))
+    assert_array_almost_equal(df_dx, df_dx_analytic, 5)
+
+    df_dy = first_derivative(deriv_2d_data.f, x=deriv_2d_data.y, axis=0)
+    # Repeat each row, then flip to get variation along rows
+    df_dy_analytic = np.tile(2 * deriv_2d_data.b * (deriv_2d_data.y - deriv_2d_data.y0),
+                             (deriv_2d_data.f.shape[1], 1)).T
+    assert_array_almost_equal(df_dy, df_dy_analytic, 5)
+
+
+def test_first_derivative_too_small(deriv_1d_data):
+    """Test first_derivative with too small an array."""
+    with pytest.raises(ValueError):
+        first_derivative(deriv_1d_data.values[None, :].T, x=deriv_1d_data.x, axis=1)
+
+
+def test_first_derivative_scalar_delta():
+    """Test first_derivative with a scalar passed for a delta."""
+    df_dx = first_derivative(np.arange(3), delta=1)
+    assert_array_almost_equal(df_dx, np.array([1., 1., 1.]), 6)
+
+
+def test_second_derivative(deriv_1d_data):
+    """Test second_derivative with a simple 1D array."""
+    d2v_dx2 = second_derivative(deriv_1d_data.values, x=deriv_1d_data.x)
+
+    # Worked by hand
+    truth = np.ones_like(deriv_1d_data.values) * 0.2133333 * units('delta_degC/cm**2')
+    assert_array_almost_equal(d2v_dx2, truth, 5)
+
+
+def test_second_derivative_2d(deriv_2d_data):
+    """Test second_derivative with a full 2D array."""
+    df2_dx2 = second_derivative(deriv_2d_data.f, x=deriv_2d_data.x, axis=1)
+    assert_array_almost_equal(df2_dx2,
+                              np.ones_like(deriv_2d_data.f) * (2 * deriv_2d_data.a), 5)
+
+    df2_dy2 = second_derivative(deriv_2d_data.f, x=deriv_2d_data.y, axis=0)
+    assert_array_almost_equal(df2_dy2,
+                              np.ones_like(deriv_2d_data.f) * (2 * deriv_2d_data.b), 5)
+
+
+def test_second_derivative_too_small(deriv_1d_data):
+    """Test second_derivative with too small an array."""
+    with pytest.raises(ValueError):
+        second_derivative(deriv_1d_data.values[None, :].T, x=deriv_1d_data.x, axis=1)
+
+
+def test_second_derivative_scalar_delta():
+    """Test second_derivative with a scalar passed for a delta."""
+    df_dx = second_derivative(np.arange(3), delta=1)
+    assert_array_almost_equal(df_dx, np.array([0., 0., 0.]), 6)
+
+
+def test_laplacian(deriv_1d_data):
+    """Test laplacian with simple 1D data."""
+    laplac = laplacian(deriv_1d_data.values, x=(deriv_1d_data.x,))
+
+    # Worked by hand
+    truth = np.ones_like(deriv_1d_data.values) * 0.2133333 * units('delta_degC/cm**2')
+    assert_array_almost_equal(laplac, truth, 5)
+
+
+def test_laplacian_2d(deriv_2d_data):
+    """Test lapacian with full 2D arrays."""
+    laplac_true = 2 * (np.ones_like(deriv_2d_data.f) * (deriv_2d_data.a + deriv_2d_data.b))
+    laplac = laplacian(deriv_2d_data.f, x=(deriv_2d_data.y, deriv_2d_data.x))
+    assert_array_almost_equal(laplac, laplac_true, 5)
diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py
index 1d8da7e3396..4a3fb1b7599 100644
--- a/metpy/calc/thermo.py
+++ b/metpy/calc/thermo.py
@@ -596,7 +596,7 @@ def equivalent_potential_temperature(pressure, temperature, dewpoint):
     Notes
     -----
     [Bolton1980]_ formula for Theta-e is used, since according to
-    [Davies-Jones2009]_ it is the most accurate non-iterative formulation
+    [DaviesJones2009]_ it is the most accurate non-iterative formulation
     available.
 
     """
diff --git a/metpy/calc/tools.py b/metpy/calc/tools.py
index 7c90123d1d1..fbc2f0dc338 100644
--- a/metpy/calc/tools.py
+++ b/metpy/calc/tools.py
@@ -2,17 +2,24 @@
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Contains a collection of generally useful calculation tools."""
+from __future__ import division
 
 import functools
 import warnings
 
 import numpy as np
+try:
+    from numpy.core.numeric import normalize_axis_index
+except ImportError:  # Only available in numpy >=1.13.0
+    def normalize_axis_index(a, n):
+        """No op version of :func:`numpy.core.numeric.normalize_axis_index`."""
+        return a
 import numpy.ma as ma
 from scipy.spatial import cKDTree
 
 from . import height_to_pressure_std, pressure_to_height_std
 from ..package_tools import Exporter
-from ..units import check_units, units
+from ..units import atleast_1d, check_units, concatenate, diff, units
 
 exporter = Exporter(globals())
 
@@ -831,3 +838,290 @@ def _less_or_close(a, value, **kwargs):
 
     """
     return np.less(a, value) | np.isclose(a, value, **kwargs)
+
+
+@exporter.export
+def first_derivative(f, **kwargs):
+    """Calculate the first derivative of a grid of values.
+
+    Works for both regularly-spaced data and grids with varying spacing.
+
+    Either `x` or `delta` must be specified. This uses 3 points to calculate the
+    derivative, using forward or backward at the edges of the grid as appropriate, and
+    centered elsewhere. The irregular spacing is handled explicitly, using the formulation
+    as specified by [Bowen2005]_.
+
+    Parameters
+    ----------
+    f : array-like
+        Array of values of which to calculate the derivative
+    axis : int, optional
+        The array axis along which to take the derivative. Defaults to 0.
+    x : array-like, optional
+        The coordinate values corresponding to the grid points in `f`.
+    delta : array-like, optional
+        Spacing between the grid points in `f`. Should be one item less than the size
+        of `f` along `axis`.
+
+    Returns
+    -------
+    array-like
+        The first derivative calculated along the selected axis.
+
+    See Also
+    --------
+    second_derivative
+
+    """
+    n, axis, delta = _process_deriv_args(f, kwargs)
+
+    # create slice objects --- initially all are [:, :, ..., :]
+    slice0 = [slice(None)] * n
+    slice1 = [slice(None)] * n
+    slice2 = [slice(None)] * n
+    delta_slice0 = [slice(None)] * n
+    delta_slice1 = [slice(None)] * n
+
+    # First handle centered case
+    slice0[axis] = slice(None, -2)
+    slice1[axis] = slice(1, -1)
+    slice2[axis] = slice(2, None)
+    delta_slice0[axis] = slice(None, -1)
+    delta_slice1[axis] = slice(1, None)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    delta_diff = delta[delta_slice1] - delta[delta_slice0]
+    center = (- delta[delta_slice1] / (combined_delta * delta[delta_slice0]) * f[slice0] +
+              delta_diff / (delta[delta_slice0] * delta[delta_slice1]) * f[slice1] +
+              delta[delta_slice0] / (combined_delta * delta[delta_slice1]) * f[slice2])
+
+    # Fill in "left" edge with forward difference
+    slice0[axis] = slice(None, 1)
+    slice1[axis] = slice(1, 2)
+    slice2[axis] = slice(2, 3)
+    delta_slice0[axis] = slice(None, 1)
+    delta_slice1[axis] = slice(1, 2)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    big_delta = combined_delta + delta[delta_slice0]
+    left = (- big_delta / (combined_delta * delta[delta_slice0]) * f[slice0] +
+            combined_delta / (delta[delta_slice0] * delta[delta_slice1]) * f[slice1] -
+            delta[delta_slice0] / (combined_delta * delta[delta_slice1]) * f[slice2])
+
+    # Now the "right" edge with backward difference
+    slice0[axis] = slice(-3, -2)
+    slice1[axis] = slice(-2, -1)
+    slice2[axis] = slice(-1, None)
+    delta_slice0[axis] = slice(-2, -1)
+    delta_slice1[axis] = slice(-1, None)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    big_delta = combined_delta + delta[delta_slice1]
+    right = (delta[delta_slice1] / (combined_delta * delta[delta_slice0]) * f[slice0] -
+             combined_delta / (delta[delta_slice0] * delta[delta_slice1]) * f[slice1] +
+             big_delta / (combined_delta * delta[delta_slice1]) * f[slice2])
+
+    return concatenate((left, center, right), axis=axis)
+
+
+@exporter.export
+def second_derivative(f, **kwargs):
+    """Calculate the second derivative of a grid of values.
+
+    Works for both regularly-spaced data and grids with varying spacing.
+
+    Either `x` or `delta` must be specified.
+
+    Either `x` or `delta` must be specified. This uses 3 points to calculate the
+    derivative, using forward or backward at the edges of the grid as appropriate, and
+    centered elsewhere. The irregular spacing is handled explicitly, using the formulation
+    as specified by [Bowen2005]_.
+
+    Parameters
+    ----------
+    f : array-like
+        Array of values of which to calculate the derivative
+    axis : int, optional
+        The array axis along which to take the derivative. Defaults to 0.
+    x : array-like, optional
+        The coordinate values corresponding to the grid points in `f`.
+    delta : array-like, optional
+        Spacing between the grid points in `f`. There should be one item less than the size
+        of `f` along `axis`.
+
+    Returns
+    -------
+    array-like
+        The second derivative calculated along the selected axis.
+
+    See Also
+    --------
+    first_derivative
+
+    """
+    n, axis, delta = _process_deriv_args(f, kwargs)
+
+    # create slice objects --- initially all are [:, :, ..., :]
+    slice0 = [slice(None)] * n
+    slice1 = [slice(None)] * n
+    slice2 = [slice(None)] * n
+    delta_slice0 = [slice(None)] * n
+    delta_slice1 = [slice(None)] * n
+
+    # First handle centered case
+    slice0[axis] = slice(None, -2)
+    slice1[axis] = slice(1, -1)
+    slice2[axis] = slice(2, None)
+    delta_slice0[axis] = slice(None, -1)
+    delta_slice1[axis] = slice(1, None)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    center = 2 * (f[slice0] / (combined_delta * delta[delta_slice0]) -
+                  f[slice1] / (delta[delta_slice0] * delta[delta_slice1]) +
+                  f[slice2] / (combined_delta * delta[delta_slice1]))
+
+    # Fill in "left" edge
+    slice0[axis] = slice(None, 1)
+    slice1[axis] = slice(1, 2)
+    slice2[axis] = slice(2, 3)
+    delta_slice0[axis] = slice(None, 1)
+    delta_slice1[axis] = slice(1, 2)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    left = 2 * (f[slice0] / (combined_delta * delta[delta_slice0]) -
+                f[slice1] / (delta[delta_slice0] * delta[delta_slice1]) +
+                f[slice2] / (combined_delta * delta[delta_slice1]))
+
+    # Now the "right" edge
+    slice0[axis] = slice(-3, -2)
+    slice1[axis] = slice(-2, -1)
+    slice2[axis] = slice(-1, None)
+    delta_slice0[axis] = slice(-2, -1)
+    delta_slice1[axis] = slice(-1, None)
+
+    combined_delta = delta[delta_slice0] + delta[delta_slice1]
+    right = 2 * (f[slice0] / (combined_delta * delta[delta_slice0]) -
+                 f[slice1] / (delta[delta_slice0] * delta[delta_slice1]) +
+                 f[slice2] / (combined_delta * delta[delta_slice1]))
+
+    return concatenate((left, center, right), axis=axis)
+
+
+def gradient(f, **kwargs):
+    """Calculate the gradient of a grid of values.
+
+    Works for both regularly-spaced data, and grids with varying spacing.
+
+    Either `x` or `deltas` must be specified.
+
+    Parameters
+    ----------
+    f : array-like
+        Array of values of which to calculate the derivative
+    x : array-like, optional
+        The coordinate values corresponding to the grid points in `f`
+    deltas : array-like, optional
+        Spacing between the grid points in `f`. There should be one item less than the size
+        of `f` along `axis`.
+
+    Returns
+    -------
+    array-like
+        The first derivative calculated along each axis in the original array
+
+    See Also
+    --------
+    laplacian
+
+    """
+    pos_kwarg, positions = _process_gradient_args(kwargs)
+    return tuple(first_derivative(f, axis=ind, **{pos_kwarg: positions[ind]})
+                 for ind, pos in enumerate(positions))
+
+
+@exporter.export
+def laplacian(f, **kwargs):
+    """Calculate the laplacian of a grid of values.
+
+    Works for both regularly-spaced data, and grids with varying spacing.
+
+    Either `x` or `deltas` must be specified.
+
+    Parameters
+    ----------
+    f : array-like
+        Array of values of which to calculate the derivative
+    x : array-like, optional
+        The coordinate values corresponding to the grid points in `f`
+    deltas : array-like, optional
+        Spacing between the grid points in `f`. There should be one item less than the size
+        of `f` along `axis`.
+
+    Returns
+    -------
+    array-like
+        The laplacian
+
+    See Also
+    --------
+    gradient
+
+    """
+    pos_kwarg, positions = _process_gradient_args(kwargs)
+    return sum(second_derivative(f, axis=ind, **{pos_kwarg: positions[ind]})
+               for ind, pos in enumerate(positions))
+
+
+def _broadcast_to_axis(arr, axis, ndim):
+    """Handle reshaping coordinate array to have proper dimensionality.
+
+    This puts the values along the specified axis.
+    """
+    if arr.ndim == 1 and arr.ndim < ndim:
+        new_shape = [1] * ndim
+        new_shape[axis] = arr.size
+        arr = arr.reshape(*new_shape)
+    return arr
+
+
+def _process_gradient_args(kwargs):
+    """Handle common processing of arguments for gradient and gradient-like functions."""
+    if 'deltas' in kwargs:
+        if 'x' in kwargs:
+            raise ValueError('Cannot specify both "x" and "deltas".')
+        return 'delta', kwargs['deltas']
+    elif 'x' in kwargs:
+        return 'x', kwargs['x']
+    else:
+        raise ValueError('Must specify either "x" or "delta" for value positions.')
+
+
+def _process_deriv_args(f, kwargs):
+    """Handle common processing of arguments for derivative functions."""
+    n = f.ndim
+    axis = normalize_axis_index(kwargs.get('axis', 0), n)
+
+    if f.shape[axis] < 3:
+        raise ValueError('f must have at least 3 point along the desired axis.')
+
+    if 'delta' in kwargs:
+        if 'x' in kwargs:
+            raise ValueError('Cannot specify both "x" and "delta".')
+
+        delta = atleast_1d(kwargs['delta'])
+        if delta.size == 1:
+            diff_size = list(f.shape)
+            diff_size[axis] -= 1
+            delta_units = getattr(delta, 'units', None)
+            delta = np.broadcast_to(delta, diff_size, subok=True)
+            if delta_units is not None:
+                delta = delta * delta_units
+        else:
+            delta = _broadcast_to_axis(delta, axis, n)
+    elif 'x' in kwargs:
+        x = _broadcast_to_axis(kwargs['x'], axis, n)
+        delta = diff(x, axis=axis)
+    else:
+        raise ValueError('Must specify either "x" or "delta" for value positions.')
+
+    return n, axis, delta
diff --git a/metpy/tests/test_units.py b/metpy/tests/test_units.py
index 2738031865c..c0b22bc63cf 100644
--- a/metpy/tests/test_units.py
+++ b/metpy/tests/test_units.py
@@ -6,10 +6,11 @@
 import sys
 
 import matplotlib.pyplot as plt
+import numpy as np
 import pytest
 
-from metpy.testing import set_agg_backend  # noqa: F401
-from metpy.units import check_units, units
+from metpy.testing import assert_array_equal, set_agg_backend  # noqa: F401
+from metpy.units import atleast_1d, atleast_2d, check_units, diff, units
 
 
 @pytest.mark.mpl_image_compare(tolerance=0, remove_text=True)
@@ -31,6 +32,23 @@ def test_axvline():
     ax.set_xlabel('')
     return fig
 
+
+def test_atleast1d_without_units():
+    """Test that atleast_1d wrapper can handle plain arrays."""
+    assert_array_equal(atleast_1d(1), np.array([1]))
+
+
+def test_atleast2d_without_units():
+    """Test that atleast_2d wrapper can handle plain arrays."""
+    assert_array_equal(atleast_2d(1), np.array([[1]]))
+
+
+def test_units_diff():
+    """Test our diff handles units properly."""
+    assert_array_equal(diff(np.arange(20, 22) * units.degC),
+                       np.array([1]) * units.delta_degC)
+
+
 #
 # Tests for unit-checking decorator
 #
diff --git a/metpy/units.py b/metpy/units.py
index 7ec695129af..100d48c2379 100644
--- a/metpy/units.py
+++ b/metpy/units.py
@@ -69,6 +69,41 @@ def concatenate(arrs, axis=0):
     return units.Quantity(np.concatenate(data, axis=axis), dest)
 
 
+def diff(x, **kwargs):
+    """Calculate the n-th discrete difference along given axis.
+
+    Wraps :func:`numpy.diff` to handle units.
+
+    Parameters
+    ----------
+    x : array-like
+        Input data
+    n : int, optional
+        The number of times values are differenced.
+    axis : int, optional
+        The axis along which the difference is taken, default is the last axis.
+
+    Returns
+    -------
+    diff : ndarray
+        The n-th differences. The shape of the output is the same as `a`
+        except along `axis` where the dimension is smaller by `n`. The
+        type of the output is the same as that of the input.
+
+    See Also
+    --------
+    numpy.diff
+
+    """
+    ret = np.diff(x, **kwargs)
+    if hasattr(x, 'units'):
+        # Can't just use units because of how things like temperature work
+        it = x.flat
+        true_units = (next(it) - next(it)).units
+        ret = ret * true_units
+    return ret
+
+
 def atleast_1d(*arrs):
     r"""Convert inputs to arrays with at least one dimension.
 
@@ -87,12 +122,15 @@ def atleast_1d(*arrs):
         A single quantity or a list of quantities, matching the number of inputs.
 
     """
-    mags = [a.magnitude for a in arrs]
-    orig_units = [a.units for a in arrs]
+    mags = [a.magnitude if hasattr(a, 'magnitude') else a for a in arrs]
+    orig_units = [a.units if hasattr(a, 'units') else None for a in arrs]
     ret = np.atleast_1d(*mags)
     if len(mags) == 1:
-        return units.Quantity(ret, orig_units[0])
-    return [units.Quantity(m, u) for m, u in zip(ret, orig_units)]
+        if orig_units[0] is not None:
+            return units.Quantity(ret, orig_units[0])
+        else:
+            return ret
+    return [units.Quantity(m, u) if u is not None else m for m, u in zip(ret, orig_units)]
 
 
 def atleast_2d(*arrs):
@@ -113,12 +151,15 @@ def atleast_2d(*arrs):
         A single quantity or a list of quantities, matching the number of inputs.
 
     """
-    mags = [a.magnitude for a in arrs]
-    orig_units = [a.units for a in arrs]
+    mags = [a.magnitude if hasattr(a, 'magnitude') else a for a in arrs]
+    orig_units = [a.units if hasattr(a, 'units') else None for a in arrs]
     ret = np.atleast_2d(*mags)
     if len(mags) == 1:
-        return units.Quantity(ret, orig_units[0])
-    return [units.Quantity(m, u) for m, u in zip(ret, orig_units)]
+        if orig_units[0] is not None:
+            return units.Quantity(ret, orig_units[0])
+        else:
+            return ret
+    return [units.Quantity(m, u) if u is not None else m for m, u in zip(ret, orig_units)]
 
 
 def masked_array(data, data_units=None, **kwargs):
diff --git a/setup.cfg b/setup.cfg
index 622029c6344..059c55918cf 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -19,6 +19,7 @@ flake8-ignore = *.py F405
 flake8-max-line-length = 95
 doctest_optionflags = NORMALIZE_WHITESPACE
 mpl-results-path = test_output
+log_print = False
 
 [doc8]
 ignore-path = docs/build,docs/api,docs/_templates