ENH: Add cumulative integration function.

src/metpy/calc/tools.py

@@ -1942,3 +1942,44 @@
     for v in variables:
         ret.append(v[~mask])
     return ret
+
+
+@exporter.export
+@xarray_derivative_wrap
+def cumulative_integrate(field, axis=None, x=None, delta=None):
+    """Return cumulative integral of field along axis.
+
+    Uses trapezoidal rule for integration.
+
+    Parameters
+    ----------
+    field : array-like
+        Array of values for which to calculate the integral
+    axis : int or str
+        The axis along which to integrate.  If `field` is an
+        `np.ndarray` or `pint.Quantity`, must be an integer.  Defaults
+        to zero.
+    x : array-like, optional
+        The coordinate values along which to integrate
+    delta : array-like, optional
+        Spacing between grid points in `field`.
+
+    Examples
+    --------
+    >>> cumulative_integrate(np.arange(5))
+    array([0. , 0.5, 2. , 4.5, 8. ])
+    >>> cumulative_integrate(xr.DataArray(np.arange(5), {"x": (("x",), np.ones(5), {"units": "m"})}, ("x",), "specific_humidity", {"units": "kg/kg"}), "x")
+    <DataArray dims: {"x": 5}
+    array([...])
+    Attrs:
+    units: kg m/kg
+    """
+    n, axis, delta = _process_deriv_args(f, axis, x, delta)
+    take = make_take(n, axis)
+
+    right = np.cumsum(
+        0.5 * (field[take(slice(1, None))] + field[take(slice(None, -1))]) * delta
+    )
+    left = np.zeros_like(field[take(0)])
+
+    return concatenate((left, right), axis=axis)

tests/calc/test_calc_tools.py

@@ -16,7 +16,8 @@
                         first_derivative, geospatial_gradient, get_layer, get_layer_heights,
                         gradient, laplacian, lat_lon_grid_deltas, nearest_intersection_idx,
                         parse_angle, pressure_to_height_std, reduce_point_density,
-                        resample_nn_1d, second_derivative, vector_derivative)
+                        resample_nn_1d, second_derivative, vector_derivative,
+                        cumulative_integrate)
 from metpy.calc.tools import (_delete_masked_points, _get_bound_pressure_height,
                               _greater_or_close, _less_or_close, _next_non_masked_element,
                               _remove_nans, azimuth_range_to_lat_lon, BASE_DEGREE_MULTIPLIER,
@@ -1557,3 +1558,30 @@
         u, v, longitude=lons, latitude=lats, crs=crs, return_only=return_only)
 
     assert len(ddx) == length
+
+
+def test_cumulative_integrate_numpy():
+    field = np.arange(5)
+    integral = cumulative_integrate(field)
+    assert integral == pytest.approx(np.array([0, 0.5, 2, 4.5, 8]))
+
+
+def test_cumulative_integrate_pint():
+    field = np.arange(6) * units("kg/m^3")
+    delta = np.array([1, 2, 3, 2, 1]) * units("cm")
+    integral = cumulative_integrate(field, delta=delta)
+    assert integral == pytest.approx(
+        np.array([0, 0.5, 3.5, 11, 18, 22.5]) / 100 * units("kg/m^2")
+    )
+
+
+def test_cumulative_integrate_xarray():
+    field = xr.DataArray(
+        np.arange(10) / 100,
+        {"x": (("x",), np.arange(100, 1001, 100), {"units": "hPa"})},
+        attrs={"units": "g/kg"}
+    )
+    integral = cumulative_integrate(field, "x")
+    assert integral.values == pytest.approx(
+        np.array([0, 0.5, 2, 4.5, 8, 12.5, 18, 24.5, 32, 40.5, 50])
+    )