lat_lon_grid_deltas (and in turn xarray grid argument magic) fails silently on xarray input with missing lat/lon units due to incorrect unit conversion

[jthielen]

In the process of helping someone troubleshoot, I came across an example for which MetPy's automatic xarray grid argument handling failed and gave an all-nan result:

(simplified and cleaned up from original wrf-python output)

coords = {
    'longitude': (('y', 'x'), [
        [-110.52898 , -110.46803 , -110.407074],
        [-110.53708 , -110.476074, -110.415085],
        [-110.54518 , -110.484146, -110.42311 ]
    ]),
    'latitude': (('y', 'x'), [
        [30.301296, 30.30827 , 30.315193],
        [30.353909, 30.360882, 30.36782 ],
        [30.40654 , 30.41352 , 30.42046 ]
    ])
}
u = xr.DataArray(
    [
        [-1.1592122 , -0.9725316 , -0.79281145],
        [-1.1712375 , -0.9994456 , -0.8228045 ],
        [-1.2131113 , -1.0486624 , -0.8700033 ]
    ],
    dims=('y', 'x'),
    coords=coords,
    attrs={'units': 'm s-1'}
)
v = xr.DataArray(
    [
        [-1.7016519, -1.7435555, -1.7542726],
        [-1.6374593, -1.6676924, -1.6710446],
        [-1.6060846, -1.6344578, -1.6379466]
    ],
    dims=('y', 'x'),
    coords=coords,
    attrs={'units': 'm s-1'}
)
print(mpcalc.divergence(u, v))

gives

<xarray.DataArray (y: 3, x: 3)>
<Quantity([[nan nan nan]
 [nan nan nan]
 [nan nan nan]], '1 / second')>
Coordinates:
    longitude  (y, x) float64 -110.5 -110.5 -110.4 ... -110.5 -110.5 -110.4
    latitude   (y, x) float64 30.3 30.31 30.32 30.35 ... 30.37 30.41 30.41 30.42
Dimensions without coordinates: y, x

whereas the expected is

<xarray.DataArray (y: 3, x: 3)>
<Quantity([[4.58534201e-05 4.74986356e-05 4.82096719e-05]
 [3.67433861e-05 3.87269479e-05 4.01718706e-05]
 [2.91369994e-05 3.10237683e-05 3.27716528e-05]], '1 / second')>
Coordinates:
    longitude  (y, x) float64 -110.5 -110.5 -110.4 ... -110.5 -110.5 -110.4
    latitude   (y, x) float64 30.3 30.31 30.32 30.35 ... 30.37 30.41 30.41 30.42
Dimensions without coordinates: y, x

or an error/warning either about missing units or passing too large of radian values.

---

While mpcalc.lat_lon_grid_deltas assumes units of degrees when no units are given, this path is not triggered for xarray arguments, as the decorator casts them to Quantitys (which are 'dimensionless' if no units string is present on the DataArray). And so, these latitudes and longitudes are treated as if they are in radians, which makes pyproj and the rest of the grid delta calculation give nans.

Given that there is already a path for assuming degrees when units are missing (that is, non-quantity arguments), I'd lean towards moving xarray handling inside lat_lon_grid_deltas and not using the default decorator. This would allow situations like the above to work seamlessly via the same fallback. However, I also see the point that this can be viewed as incorrect metadata on the input, and we shouldn't silently guess that it is degrees. In that case, we should probably include something like the _check_radians helper inside lat_lon_grid_deltas. Alternatively, if we want to enforce a more general rule of "xarray coordinates need to have correct units", we could include unit-checking in the add_grid_arguments_to_dataarray decorator (or other suitable helper).

[jthielen]

These grid argument utils are almost certainly going to affected by the changes towards #893, so would it make sense to fix this along the way in (or as a quick patch after) the eventual PR for #893?