def lat_lon_grid_deltas(longitude, latitude, **kwargs):
    r"""Calculate the delta between grid points that are in a latitude/longitude format.

    Calculate the signed delta distance between grid points when the grid spacing is defined by
    delta lat/lon rather than delta x/y

    Parameters
    ----------
    longitude : array_like
        array of longitudes defining the grid
    latitude : array_like
        array of latitudes defining the grid
    kwargs
        Other keyword arguments to pass to :class:`~pyproj.Geod`

    Returns
    -------
    dx, dy:
        at least two dimensional arrays of signed deltas between grid points in the x and y
        direction

    Notes
    -----
    Accepts 1D, 2D, or higher arrays for latitude and longitude
    Assumes [..., Y, X] for >=2 dimensional arrays

    """
    from pyproj import Geod

    # Inputs must be the same number of dimensions
    if latitude.ndim != longitude.ndim:
        raise ValueError('Latitude and longitude must have the same number of dimensions.')

    # If we were given 1D arrays, make a mesh grid
    if latitude.ndim < 2:
        longitude, latitude = np.meshgrid(longitude, latitude)

    geod_args = {'ellps': 'sphere'}
    if kwargs:
        geod_args = kwargs

    g = Geod(**geod_args)

    forward_az, _, dy = g.inv(longitude[..., :-1, :], latitude[..., :-1, :],
                              longitude[..., 1:, :], latitude[..., 1:, :])
    dy[(forward_az < -90.) | (forward_az > 90.)] *= -1

    forward_az, _, dx = g.inv(longitude[..., :, :-1], latitude[..., :, :-1],
                              longitude[..., :, 1:], latitude[..., :, 1:])
    dx[(forward_az < 0.) | (forward_az > 180.)] *= -1

    return dx * units.meter, dy * units.meter