Merge pull request #1651 from jthielen/issue-1648

Remove misleading fallback to lat/lon CRS when only 2D lat/lon present

src/metpy/xarray.py

@@ -24,7 +24,7 @@
 import warnings
 
 import numpy as np
-from pyproj import Proj
+from pyproj import CRS, Proj
 import xarray as xr
 
 from ._vendor.xarray import either_dict_or_kwargs, expanded_indexer, is_dict_like
@@ -749,16 +749,26 @@ def parse_cf(self, varname=None, coordinates=None):
 
         if crs is None and not check_axis(var, 'latitude', 'longitude'):
             # This isn't a lat or lon coordinate itself, so determine if we need to fall back
-            # to creating a latitude_longitude CRS. We do so if there exists valid coordinates
-            # for latitude and longitude, even if they are not the dimension coordinates of
-            # the variable.
-            def _has_coord(coord_type):
-                return any(check_axis(coord_var, coord_type)
-                           for coord_var in var.coords.values())
-            if _has_coord('latitude') and _has_coord('longitude'):
-                crs = CFProjection({'grid_mapping_name': 'latitude_longitude'})
-                log.warning('Found valid latitude/longitude coordinates, assuming '
-                            'latitude_longitude for projection grid_mapping variable')
+            # to creating a latitude_longitude CRS. We do so if there exists valid *at most
+            # 1D* coordinates for latitude and longitude (usually dimension coordinates, but
+            # that is not strictly required, for example, for DSG's). What is required is that
+            # x == latitude and y == latitude (so that all assumptions about grid coordinates
+            # and CRS line up).
+            try:
+                latitude, y, longitude, x = var.metpy.coordinates(
+                    'latitude',
+                    'y',
+                    'longitude',
+                    'x'
+                )
+            except AttributeError:
+                # This means that we don't even have sufficient coordinates, so skip
+                pass
+            else:
+                if latitude.identical(y) and longitude.identical(x):
+                    crs = CFProjection({'grid_mapping_name': 'latitude_longitude'})
+                    log.warning('Found valid latitude/longitude coordinates, assuming '
+                                'latitude_longitude for projection grid_mapping variable')
 
         # Rebuild the coordinates of the dataarray, and return quantified DataArray
         var = self._rebuild_coords(var, crs)
@@ -1355,9 +1365,17 @@ def grid_deltas_from_dataarray(f, kind='default'):
     from metpy.calc import lat_lon_grid_deltas
 
     # Determine behavior
-    if kind == 'default' and f.metpy.crs['grid_mapping_name'] == 'latitude_longitude':
+    if (
+        kind == 'default'
+        and (
+            not hasattr(f.metpy, 'crs')
+            or f.metpy.crs['grid_mapping_name'] == 'latitude_longitude'
+        )
+    ):
+        # Use actual grid deltas by default with latitude_longitude or missing CRS
         kind = 'actual'
     elif kind == 'default':
+        # Otherwise, use grid deltas in projected grid space by default
         kind = 'nominal'
     elif kind not in ('actual', 'nominal'):
         raise ValueError('"kind" argument must be specified as "default", "actual", or '
@@ -1373,11 +1391,17 @@ def grid_deltas_from_dataarray(f, kind='default'):
             warnings.warn('y and x dimensions unable to be identified. Assuming [..., y, x] '
                           'dimension order.')
             y_dim, x_dim = -2, -1
+
+        # Get geod if it exists, otherwise fall back to PyProj default
+        try:
+            geod = f.metpy.pyproj_crs.get_geod()
+        except AttributeError:
+            geod = CRS.from_cf({'grid_mapping_name': 'latitude_longitude'}).get_geod()
         # Obtain grid deltas as xarray Variables
         (dx_var, dx_units), (dy_var, dy_units) = (
             (xr.Variable(dims=latitude.dims, data=deltas.magnitude), deltas.units)
             for deltas in lat_lon_grid_deltas(longitude, latitude, x_dim=x_dim, y_dim=y_dim,
-                                              geod=f.metpy.pyproj_crs.get_geod()))
+                                              geod=geod))
     else:
         # Obtain y/x coordinate differences
         y, x = f.metpy.coordinates('y', 'x')

tests/test_xarray.py

@@ -220,8 +220,8 @@ def test_radian_projection_coords():
     assert data_var.coords['y'].metpy.unit_array[1] == 3 * units.meter
 
 
-def test_missing_grid_mapping():
-    """Test falling back to implicit lat/lon projection."""
+def test_missing_grid_mapping_valid():
+    """Test falling back to implicit lat/lon projection when valid."""
     lon = xr.DataArray(-np.arange(3),
                        attrs={'standard_name': 'longitude', 'units': 'degrees_east'})
     lat = xr.DataArray(np.arange(2),
@@ -230,7 +230,16 @@ def test_missing_grid_mapping():
     ds = xr.Dataset({'data': data})
 
     data_var = ds.metpy.parse_cf('data')
-    assert 'metpy_crs' in data_var.coords
+    assert (
+        'metpy_crs' in data_var.coords
+        and data_var.metpy.crs['grid_mapping_name'] == 'latitude_longitude'
+    )
+
+
+def test_missing_grid_mapping_invalid(test_var_multidim_no_xy):
+    """Test not falling back to implicit lat/lon projection when invalid."""
+    data_var = test_var_multidim_no_xy.to_dataset(name='data').metpy.parse_cf('data')
+    assert 'metpy_crs' not in data_var.coords
 
 
 def test_missing_grid_mapping_var(caplog):