Repository: Unidata/MetPy

Issue:
`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
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)

```python
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](https://github.com/Unidata/MetPy/blob/7ea29cad33f881010b5c63b45b2dc3f30054ea88/src/metpy/calc/tools.py#L831-L837), this path is not triggered for xarray arguments, as the decorator casts them to `Quantity`s (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 `nan`s.

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).




Discussion:

Patch:
diff --git a/src/metpy/xarray.py b/src/metpy/xarray.py
index 545122c4c5..5559d97608 100644
--- a/src/metpy/xarray.py
+++ b/src/metpy/xarray.py
@@ -134,7 +134,12 @@ class MetPyDataArrayAccessor:
         if is_quantity(self._data_array.variable._data):
             return self._data_array.variable._data.units
         else:
-            return units.parse_units(self._data_array.attrs.get('units', 'dimensionless'))
+            axis = self._data_array.attrs.get('_metpy_axis', '')
+            if 'latitude' in axis or 'longitude' in axis:
+                default_unit = 'degrees'
+            else:
+                default_unit = 'dimensionless'
+            return units.parse_units(self._data_array.attrs.get('units', default_unit))
 
     @property
     def magnitude(self):


Test patch:
diff --git a/tests/test_xarray.py b/tests/test_xarray.py
index 8f336eab3e..00d9e16ea0 100644
--- a/tests/test_xarray.py
+++ b/tests/test_xarray.py
@@ -173,6 +173,20 @@ def test_convert_coordinate_units(test_ds_generic):
     assert result['b'].metpy.units == units.percent
 
 
+def test_latlon_default_units(test_var_multidim_full):
+    """Test that lat/lon are given degree units by default."""
+    del test_var_multidim_full.lat.attrs['units']
+    del test_var_multidim_full.lon.attrs['units']
+
+    lat = test_var_multidim_full.metpy.latitude.metpy.unit_array
+    assert lat.units == units.degrees
+    assert lat.max() > 50 * units.degrees
+
+    lon = test_var_multidim_full.metpy.longitude.metpy.unit_array
+    assert lon.units == units.degrees
+    assert lon.min() < -100 * units.degrees
+
+
 def test_quantify(test_ds_generic):
     """Test quantify method for converting data to Quantity."""
     original = test_ds_generic['test'].values