How to use the metpy.calc.lat_lon_grid_deltas function in MetPy

[2.14238865, 0.88207403, -0.40652485, -1.67244834, -2.86837275]],
                       [[-1.99024801, -2.59146057, -3.04973279, -3.3296825, -3.42137476],
                        [-1.8711102, -2.71865804, -3.45952099, -4.05064148, -4.45309013],
                        [-0.99367383, -2.04299168, -3.02642031, -3.90252563, -4.62540783],
                        [0.547778, -0.63635567, -1.80391109, -2.90776869, -3.90375721],
                        [2.33967328, 1.12072805, -0.13066324, -1.3662872, -2.5404749]],
                       [[-1.75713411, -2.35840168, -2.82756083, -3.12859972, -3.24963361],
                        [-1.6178372, -2.44835688, -3.18270452, -3.77873886, -4.19768429],
                        [-0.7468435, -1.76693701, -2.72999246, -3.596093, -4.32045429],
                        [0.7730159, -0.37427326, -1.51160943, -2.59396958, -3.57843192],
                        [2.53695791, 1.35938207, 0.14519838, -1.06012605, -2.21257705]]])\
        * (units.meters / units.sec)
    lats = np.array([57.5, 57., 56.5, 56., 55.5]) * units.degrees
    lons = np.array([227.5, 228., 228.5, 229., 229.5]) * units.degrees

    dx, dy = lat_lon_grid_deltas(lons, lats)

    pvor = potential_vorticity_baroclinic(isentlevs[:, None, None], isentprs,
                                          isentu, isentv, dx[None, :, :], dy[None, :, :],
                                          lats[None, :, None])

    true_pv = np.array([[[2.97116898e-06, 3.38486331e-06, 3.81432403e-06, 4.24722471e-06,
                          4.64995688e-06],
                         [2.04235589e-06, 2.35739554e-06, 2.71138003e-06, 3.11803005e-06,
                          3.54655984e-06],
                         [1.41179481e-06, 1.60663306e-06, 1.85439220e-06, 2.17827401e-06,
                          2.55309150e-06],
                         [1.25933892e-06, 1.31915377e-06, 1.43444064e-06, 1.63067920e-06,
                          1.88631658e-06],
                         [1.37533104e-06, 1.31658998e-06, 1.30424716e-06, 1.36777872e-06,
                          1.49289942e-06]],
                        [[3.07674708e-06, 3.48172482e-06, 3.90371030e-06, 4.33207155e-06,

[17., 17., 17., 17.],
                      [16., 16., 16., 16.]],
                     [[12., 11., 10., 9.],
                      [11., 10.5, 10.5, 10.],
                      [10., 10.5, 10.5, 11.],
                      [9., 10., 11., 12.]],
                     [[-10., -10., -10., -10.],
                      [-10., -10., -10., -10.],
                      [-11., -11., -11., -11.],
                      [-11., -11., -11., -11.]]]) * units('degC')

    p = np.array([850., 700., 500.]) * units('hPa')

    lats = np.linspace(35., 40., 4) * units('degrees')
    lons = np.linspace(-100., -90., 4) * units('degrees')
    dx, dy = lat_lon_grid_deltas(lons, lats)

    return u, v, temp, p, dx, dy

lon = np.arange(-100, -90, 2.5)
    dx_truth = np.array([[212943.5585, 212943.5585, 212943.5585],
                         [204946.2305, 204946.2305, 204946.2305],
                         [196558.8269, 196558.8269, 196558.8269],
                         [187797.3216, 187797.3216, 187797.3216]]) * units.meter
    dy_truth = np.array([[277987.1857, 277987.1857, 277987.1857, 277987.1857],
                         [277987.1857, 277987.1857, 277987.1857, 277987.1857],
                         [277987.1857, 277987.1857, 277987.1857, 277987.1857]]) * units.meter
    if flip_order:
        lon = lon[::-1]
        lat = lat[::-1]
        dx_truth = -1 * dx_truth[::-1]
        dy_truth = -1 * dy_truth[::-1]

    lon, lat = np.meshgrid(lon, lat)
    dx, dy = lat_lon_grid_deltas(lon, lat)
    assert_almost_equal(dx, dx_truth, 4)
    assert_almost_equal(dy, dy_truth, 4)

def test_lat_lon_grid_deltas_mismatched_shape():
    """Test for lat_lon_grid_deltas for variable grid."""
    lat = np.arange(40, 50, 2.5)
    lon = np.array([[-100., -97.5, -95., -92.5],
                    [-100., -97.5, -95., -92.5],
                    [-100., -97.5, -95., -92.5],
                    [-100., -97.5, -95., -92.5]])
    with pytest.raises(ValueError):
        lat_lon_grid_deltas(lon, lat)

# ----------
#
# Use Xarray to access GFS data from THREDDS resource and uses
# metpy accessor to parse file to make it easy to pull data using
# common coordinate names (e.g., vertical) and attach units.
#

ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/'
                     'casestudies/python-gallery/NARR_19930313_1800.nc').metpy.parse_cf()

# Get lat/lon data from file
lats = ds.lat.data
lons = ds.lon.data

# Calculate variable dx, dy values for use in calculations
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)

# Get 700-hPa data and smooth
level = 700 * units.hPa
hght_700 = mpcalc.smooth_n_point(ds['Geopotential_height_isobaric'].metpy.sel(
    vertical=level).squeeze(), 9)
tmpk_700 = mpcalc.smooth_n_point(ds['Temperature_isobaric'].metpy.sel(
    vertical=level).squeeze(), 9)
uwnd_700 = mpcalc.smooth_n_point(
    ds['u-component_of_wind_isobaric'].metpy.sel(vertical=level).squeeze(), 9)
vwnd_700 = mpcalc.smooth_n_point(
    ds['v-component_of_wind_isobaric'].metpy.sel(vertical=level).squeeze(), 9)

# Get 300-hPa data and
level = 300 * units.hPa
hght_300 = mpcalc.smooth_n_point(ds['Geopotential_height_isobaric'].metpy.sel(
    vertical=level).squeeze(), 9)

# Create a clean datetime object for plotting based on time of Geopotential heights
vtime = ds.time.data[0].astype('datetime64[ms]').astype('O')


######################################################################
# MetPy Absolute Vorticity Calculation
# ------------------------------------
#
# This code first uses MetPy to calcualte the grid deltas (sign aware) to
# use for derivative calculations with the funtcion
# ``lat_lon_grid_deltas()`` and then calculates ``absolute_vorticity()``
# using the wind components, grid deltas, and latitude values.
#

# Calculate grid spacing that is sign aware to use in absolute vorticity calculation
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)

# Calculate absolute vorticity from MetPy function
avor_500 = mpcalc.absolute_vorticity(uwnd_er, vwnd_er, dx, dy, lats * units.degrees,
                                     dim_order='yx')


######################################################################
# Map Creation
# ------------
#
# This next set of code creates the plot and draws contours on a Lambert
# Conformal map centered on -100 E longitude. The main view is over the
# CONUS with geopotential heights contoured every 60 m and absolute
# vorticity colorshaded (:math:`*10^5`).
#

######################################################################
# Calculate frontogenesis
# -----------------------
#
# Frontogenesis calculation in MetPy requires temperature, wind
# components, and grid spacings. First compute the grid deltas using MetPy
# functionality, then put it all together in the frontogenesis function.
#
# Note: MetPy will give the output with SI units, but typically
# frontogenesis (read: GEMPAK) output this variable with units of K per
# 100 km per 3 h; a conversion factor is included here to use at plot time
# to reflect those units.
#

dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)

fronto_850 = mpcalc.frontogenesis(thta_850, uwnd_850, vwnd_850, dx, dy, dim_order='yx')

# A conversion factor to get frontogensis units of K per 100 km per 3 h
convert_to_per_100km_3h = 1000*100*3600*3


######################################################################
# Plotting Frontogenesis
# ----------------------
#
# Using a Lambert Conformal projection from Cartopy to plot 850-hPa
# variables including frontogenesis.
#

# Set map projection

def _metpy(func, data, lat, lon, dim):
    """Wrapper for MetPy functions

    Parameters:
    func -- the MetPy function
    data -- the xarray or netcdf variable (already sliced)
    lat -- an array of latitudes, the shape must match that of data
    lon -- an array of longitudes, the shape must match that of data
    dim -- the dimension to return, a string, x or y
    """
    if hasattr(data, "dims"):
        dims = data.dims
    else:
        dims = data.dimensions

    dx, dy = metpy.calc.lat_lon_grid_deltas(np.array(lon), np.array(lat))
    dim_order = "".join([d for d in dims if d in 'yx'])

    if dim_order == "yx":
        deltas = [dy, dx]
    else:
        deltas = [dx, dy]

    if len(dims) > 2:
        axes = list(range(0, len(dims)))
        new_axes = list(axes)
        new_dims = list(dims)
        if dim_order == 'yx':
            new_axes += [new_axes.pop(new_dims.index('y'))]
            new_dims += [new_dims.pop(new_dims.index('y'))]
            new_axes += [new_axes.pop(new_dims.index('x'))]
            new_dims += [new_dims.pop(new_dims.index('x'))]

# Calculations
# ------------
#
# Most of the calculations in `metpy.calc` will accept DataArrays by converting them
# into their corresponding unit arrays. While this may often work without any issues, we must
# keep in mind that because the calculations are working with unit arrays and not DataArrays:
#
# - The calculations will return unit arrays rather than DataArrays
# - Broadcasting must be taken care of outside of the calculation, as it would only recognize
#   dimensions by order, not name
#
# As an example, we calculate geostropic wind at 500 hPa below:

lat, lon = xr.broadcast(y, x)
f = mpcalc.coriolis_parameter(lat)
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat, initstring=data_crs.proj4_init)
heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 500. * units.hPa}]
u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy)
print(u_geo)
print(v_geo)

#########################################################################
# Also, a limited number of calculations directly support xarray DataArrays or Datasets (they
# can accept *and* return xarray objects). Right now, this includes
#
# - Derivative functions
#     - ``first_derivative``
#     - ``second_derivative``
#     - ``gradient``
#     - ``laplacian``
# - Cross-section functions
#     - ``cross_section_components``

# Calculations
# ------------
#
# Most of the calculations in `metpy.calc` will accept DataArrays by converting them
# into their corresponding unit arrays. While this may often work without any issues, we must
# keep in mind that because the calculations are working with unit arrays and not DataArrays:
#
# - The calculations will return unit arrays rather than DataArrays
# - Broadcasting must be taken care of outside of the calculation, as it would only recognize
#   dimensions by order, not name
#
# As an example, we calculate geostropic wind at 500 hPa below:

lat, lon = xr.broadcast(y, x)
f = mpcalc.coriolis_parameter(lat)
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat, initstring=data_crs.proj4_init)
heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 500. * units.hPa}]
u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy)
print(u_geo)
print(v_geo)

#########################################################################
# Also, a limited number of calculations directly support xarray DataArrays or Datasets (they
# can accept *and* return xarray objects). Right now, this includes
#
# - Derivative functions
#     - ``first_derivative``
#     - ``second_derivative``
#     - ``gradient``
#     - ``laplacian``
# - Cross-section functions
#     - ``cross_section_components``