def cross_section(data, start, end, steps=100, interp_type='linear'):
    r"""Obtain an interpolated cross-sectional slice through gridded data.

    Utilizing the interpolation functionality in `xarray`, this function takes a vertical
    cross-sectional slice along a geodesic through the given data on a regular grid, which is
    given as an `xarray.DataArray` so that we can utilize its coordinate and projection
    metadata.

    Parameters
    ----------
    data: `xarray.DataArray` or `xarray.Dataset`
        Three- (or higher) dimensional field(s) to interpolate. The DataArray (or each
        DataArray in the Dataset) must have been parsed by MetPy and include both an x and
        y coordinate dimension and the added `crs` coordinate.
    start: (2, ) array_like
        A latitude-longitude pair designating the start point of the cross section (units are
        degrees north and degrees east).
    end: (2, ) array_like
        A latitude-longitude pair designating the end point of the cross section (units are
        degrees north and degrees east).
    steps: int, optional
        The number of points along the geodesic between the start and the end point
        (including the end points) to use in the cross section. Defaults to 100.
    interp_type: str, optional
        The interpolation method, either 'linear' or 'nearest' (see
        `xarray.DataArray.interp()` for details). Defaults to 'linear'.

    Returns
    -------
    `xarray.DataArray` or `xarray.Dataset`
        The interpolated cross section, with new index dimension along the cross-section.

    See Also
    --------
    interpolate_to_slice, geodesic

    """
    if isinstance(data, xr.Dataset):
        # Recursively apply to dataset
        return data.apply(cross_section, True, (start, end), steps=steps,
                          interp_type=interp_type)
    elif data.ndim == 0:
        # This has no dimensions, so it is likely a projection variable. In any case, there
        # are no data here to take the cross section with. Therefore, do nothing.
        return data
    else:

        # Get the projection and coordinates
        try:
            crs_data = data.metpy.cartopy_crs
            x = data.metpy.x
        except AttributeError:
            raise ValueError('Data missing required coordinate information. Verify that '
                             'your data have been parsed by MetPy with proper x and y '
                             'dimension coordinates and added crs coordinate of the '
                             'correct projection for each variable.')

        # Get the geodesic
        points_cross = geodesic(crs_data, start, end, steps)

        # Patch points_cross to match given longitude range, whether [0, 360) or (-180,  180]
        if CFConventionHandler.check_axis(x, 'lon') and (x > 180).any():
            points_cross[points_cross[:, 0] < 0, 0] += 360.

        # Return the interpolated data
        return interpolate_to_slice(data, points_cross, interp_type=interp_type)
r"""Obtain an interpolated cross-sectional slice through gridded data.

    Utilizing the interpolation functionality in `xarray`, this function takes a vertical
    cross-sectional slice along a geodesic through the given data on a regular grid, which is
    given as an `xarray.DataArray` so that we can utilize its coordinate and projection
    metadata.

    Parameters
    ----------
    data: `xarray.DataArray` or `xarray.Dataset`
        Three- (or higher) dimensional field(s) to interpolate. The DataArray (or each
        DataArray in the Dataset) must have been parsed by MetPy and include both an x and
        y coordinate dimension and the added `crs` coordinate.
    start: (2, ) array_like
        A latitude-longitude pair designating the start point of the cross section (units are
        degrees north and degrees east).
    end: (2, ) array_like
        A latitude-longitude pair designating the end point of the cross section (units are
        degrees north and degrees east).
    steps: int, optional
        The number of points along the geodesic between the start and the end point
        (including the end points) to use in the cross section. Defaults to 100.
    interp_type: str, optional
        The interpolation method, either 'linear' or 'nearest' (see
        `xarray.DataArray.interp()` for details). Defaults to 'linear'.

    Returns
    -------
    `xarray.DataArray` or `xarray.Dataset`
        The interpolated cross section, with new index dimension along the cross-section.

    See Also
    --------
    interpolate_to_slice, geodesic