Ensuring Unit Correctness with Metpy.calc.divergence()

[Tstolarz]

Hello! I am using metpy.calc.divergence() using HFR vector fields of ocean currents for my thesis work of analyzing divergence trends. I am calculating divergence with u and v's with lat and lon coordinates. Upon calculating the divergence of the data array I arrive at data with units of 1/s. This seems well and dandy but my advisor is wondering how this package exactly converts dx and dy with degrees lat and lon to a discrete length scale distance which leads to the units of 1/s if my vectors are a speed in m/s

My understanding is if I am plotting the data to a projection that already uses lat and lon coordinates (Mercator with transform=ccrs.PlateCarree on the plotting function) then it is not really a problem. But I could be wrong and I could be overlooking things. Is there anywhere that explicitly explains how metpy calculates the first_derivative dx on the longitude and dy on latitude to ensure the units are correct?

[DanielAdriaansen]

Here's the documentation of first_derivative():
https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.first_derivative.html#metpy.calc.first_derivative

And here's the actual source code that goes with it:

MetPy/src/metpy/calc/tools.py

Lines 952 to 1041 in df6f7c5

	def first_derivative(f, axis=None, x=None, delta=None):
	"""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, or `f` must be given as an `xarray.DataArray` with
	attached coordinate and projection information. If `f` is an `xarray.DataArray`, and `x` or
	`delta` are given, `f` will be converted to a `pint.Quantity` and the derivative returned
	as a `pint.Quantity`, otherwise, if neither `x` nor `delta` are given, the attached
	coordinate information belonging to `axis` will be used and the derivative will be returned
	as an `xarray.DataArray`.
	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 or str, optional
	The array axis along which to take the derivative. If `f` is ndarray-like, must be an
	integer. If `f` is a `DataArray`, can be a string (referring to either the coordinate
	dimension name or the axis type) or integer (referring to axis number), unless using
	implicit conversion to `pint.Quantity`, in which case it must be an integer. Defaults
	to 0. For reference, the current standard axis types are 'time', 'vertical', 'y', and
	'x'.
	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
	.. versionchanged:: 1.0
	Changed signature from ``(f, **kwargs)``
	See Also
	--------
	second_derivative
	"""
	n, axis, delta = _process_deriv_args(f, axis, x, delta)
	take = make_take(n, axis)
	# First handle centered case
	slice0 = take(slice(None, -2))
	slice1 = take(slice(1, -1))
	slice2 = take(slice(2, None))
	delta_slice0 = take(slice(None, -1))
	delta_slice1 = take(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 = take(slice(None, 1))
	slice1 = take(slice(1, 2))
	slice2 = take(slice(2, 3))
	delta_slice0 = take(slice(None, 1))
	delta_slice1 = take(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 = take(slice(-3, -2))
	slice1 = take(slice(-2, -1))
	slice2 = take(slice(-1, None))
	delta_slice0 = take(slice(-2, -1))
	delta_slice1 = take(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)

I think the key bit is here:

f must be given as an xarray.DataArray with attached coordinate and projection information.

So as long as your DataArray object has coordinate and projection information attached, then it will be using that.

Here's the actual reference for the mathematics used:
https://royalsocietypublishing.org/doi/10.1098/rspa.2004.1430

I did not write this code nor claim to be an expert in mathematics, but wanted to provide this information in case it is useful!

[dcamron]

Expanding on this, if you're working with lat/lon data, the actual determination of grid-spacing is coming from lat_lon_grid_deltas, which relies on Pyproj's inverse great circle calculations to calculate distances in meters between latitudes and longitudes. It will attempt to determine your relevant Geod from provided projection information, else it will assume Pyproj's default latitude_longitude Geod.