Use PyProj instead of Cartopy for internal coordinate transforms and …

…rename crs coordinate

src/metpy/calc/cross_sections.py

@@ -36,9 +36,7 @@ def distances_from_cross_section(cross):
     """
     if check_axis(cross.metpy.x, 'longitude') and check_axis(cross.metpy.y, 'latitude'):
         # Use pyproj to obtain x and y distances
-        from pyproj import Geod
-
-        g = Geod(cross.metpy.cartopy_crs.proj4_init)
+        g = cross.metpy.pyproj_crs.get_geod()
         lon = cross.metpy.x
         lat = cross.metpy.y
 
@@ -83,10 +81,13 @@ def latitude_from_cross_section(cross):
     if check_axis(y, 'latitude'):
         return y
     else:
-        import cartopy.crs as ccrs
-        latitude = ccrs.Geodetic().transform_points(cross.metpy.cartopy_crs,
-                                                    cross.metpy.x.values,
-                                                    y.values)[..., 1]
+        from pyproj import Proj
+        latitude = Proj(cross.metpy.pyproj_crs)(
+            cross.metpy.x.values,
+            y.values,
+            inverse=True,
+            radians=False
+        )[1]
         latitude = xr.DataArray(latitude * units.degrees_north, coords=y.coords, dims=y.dims)
         return latitude
 

src/metpy/calc/tools.py

@@ -751,7 +751,7 @@ def take(indexer):
 
 @exporter.export
 @preprocess_xarray
-def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs):
+def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, geod=None):
     r"""Calculate the actual delta between grid points that are in latitude/longitude format.
 
     Parameters
@@ -766,8 +766,9 @@ def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs):
         axis number for the y dimesion, defaults to -2.
     x_dim: int
         axis number for the x dimension, defaults to -1.
-    kwargs
-        Other keyword arguments to pass to :class:`~pyproj.Geod`
+    geod : pyproj.Geod or None
+        PyProj Geod to use for forward azimuth and distance calculations. If None, use a
+        default spherical ellipsoid.
 
     Returns
     -------
@@ -804,11 +805,10 @@ def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs):
     take_y = make_take(latitude.ndim, y_dim)
     take_x = make_take(latitude.ndim, x_dim)
 
-    geod_args = {'ellps': 'sphere'}
-    if kwargs:
-        geod_args = kwargs
-
-    g = Geod(**geod_args)
+    if geod is None:
+        g = Geod(ellps='sphere')
+    else:
+        g = geod
 
     forward_az, _, dy = g.inv(longitude[take_y(slice(None, -1))],
                               latitude[take_y(slice(None, -1))],
@@ -922,7 +922,7 @@ def grid_deltas_from_dataarray(f, kind='default'):
         (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, y_dim=y_dim, x_dim=x_dim,
-                                              initstring=f.metpy.cartopy_crs.proj4_init))
+                                              geod=f.metpy.pyproj_crs.get_geod()))
     else:
         # Obtain y/x coordinate differences
         y, x = f.metpy.coordinates('y', 'x')

src/metpy/interpolate/slices.py

@@ -74,8 +74,8 @@ def geodesic(crs, start, end, steps):
 
     Parameters
     ----------
-    crs: `cartopy.crs`
-        Cartopy Coordinate Reference System to use for the output
+    crs: `pyproj.CRS`
+        PyProj Coordinate Reference System to use for the output
     start: (2, ) array_like
         A latitude-longitude pair designating the start point of the geodesic (units are
         degrees north and degrees east).
@@ -96,18 +96,19 @@ def geodesic(crs, start, end, steps):
     cross_section
 
     """
-    import cartopy.crs as ccrs
-    from pyproj import Geod
+    from pyproj import Proj
+
+    g = crs.get_geod()
+    p = Proj(crs)
 
     # Geod.npts only gives points *in between* the start and end, and we want to include
     # the endpoints.
-    g = Geod(crs.proj4_init)
     geodesic = np.concatenate([
         np.array(start[::-1])[None],
         np.array(g.npts(start[1], start[0], end[1], end[0], steps - 2)),
         np.array(end[::-1])[None]
     ]).transpose()
-    points = crs.transform_points(ccrs.Geodetic(), *geodesic)[:, :2]
+    points = np.stack(p(geodesic[0], geodesic[1], inverse=False, radians=False), axis=-1)
 
     return points
 
@@ -162,7 +163,7 @@ def cross_section(data, start, end, steps=100, interp_type='linear'):
 
         # Get the projection and coordinates
         try:
-            crs_data = data.metpy.cartopy_crs
+            crs_data = data.metpy.pyproj_crs
             x = data.metpy.x
         except AttributeError:
             raise ValueError('Data missing required coordinate information. Verify that '

src/metpy/plots/mapping.py

@@ -77,6 +77,12 @@ def to_cartopy(self):
 
         return proj_handler(self._attrs, globe)
 
+    def to_pyproj(self):
+        """Convert to a PyProj CRS."""
+        import pyproj
+
+        return pyproj.CRS.from_cf(self._attrs)
+
     def to_dict(self):
         """Get the dictionary of metadata attributes."""
         return self._attrs.copy()

src/metpy/xarray.py

@@ -212,8 +212,8 @@ def dequantify(self):
     @property
     def crs(self):
         """Return the coordinate reference system (CRS) as a CFProjection object."""
-        if 'crs' in self._data_array.coords:
-            return self._data_array.coords['crs'].item()
+        if 'metpy_crs' in self._data_array.coords:
+            return self._data_array.coords['metpy_crs'].item()
         raise AttributeError('crs attribute is not available.')
 
     @property
@@ -231,6 +231,11 @@ def cartopy_geodetic(self):
         """Return the Geodetic CRS associated with the native CRS globe."""
         return self.crs.cartopy_geodetic
 
+    @property
+    def pyproj_crs(self):
+        """Return the coordinate reference system (CRS) as a pyproj object."""
+        return self.crs.to_pyproj()
+
     def _fixup_coordinate_map(self, coord_map):
         """Ensure sure we have coordinate variables in map, not coordinate names."""
         new_coord_map = {}
@@ -569,7 +574,7 @@ def assign_latitude_longitude(self, force=False):
 
         Notes
         -----
-        A valid CRS coordinate must be present. Cartopy is used for the coordinate
+        A valid CRS coordinate must be present. PyProj is used for the coordinate
         transformations.
 
         """
@@ -605,7 +610,7 @@ def assign_y_x(self, force=False, tolerance=None):
 
         Notes
         -----
-        A valid CRS coordinate must be present. Cartopy is used for the coordinate
+        A valid CRS coordinate must be present. PyProj is used for the coordinate
         transformations.
 
         """
@@ -632,7 +637,7 @@ class MetPyDatasetAccessor:
         >>> import xarray as xr
         >>> from metpy.cbook import get_test_data
         >>> ds = xr.open_dataset(get_test_data('narr_example.nc', False)).metpy.parse_cf()
-        >>> print(ds['crs'].item())
+        >>> print(ds['metpy_crs'].item())
         Projection: lambert_conformal_conic
 
     """
@@ -675,6 +680,22 @@ def parse_cf(self, varname=None, coordinates=None):
 
         var = self._dataset[varname]
 
+        # Check for crs conflict
+        if varname == 'metpy_crs':
+            warnings.warn(
+                'Attempting to parse metpy_crs as a data variable. Unexpected merge conflicts '
+                'may occur.'
+            )
+        elif 'metpy_crs' in var.coords:
+            try:
+                assert isinstance(var.coords['metpy_crs'].item(), CFProjection)
+            except (ValueError, AssertionError):
+                # Catch non-scalar and non-CFProjection coordinates
+                warnings.warn(
+                    'metpy_crs already present as a non-CFProjection coordinate. Unexpected '
+                    'merge conflicts may occur.'
+                )
+
         # Assign coordinates if the coordinates argument is given
         if coordinates is not None:
             var = var.metpy.assign_coordinates(coordinates)
@@ -709,7 +730,7 @@ def _has_coord(coord_type):
         # Rebuild the coordinates of the dataarray, and return quantified DataArray
         var = self._rebuild_coords(var, crs)
         if crs is not None:
-            var = var.assign_coords(coords={'crs': crs})
+            var = var.assign_coords(coords={'metpy_crs': crs})
         return var
 
     def _rebuild_coords(self, var, crs):
@@ -795,7 +816,7 @@ def assign_latitude_longitude(self, force=False):
 
         Notes
         -----
-        A valid CRS coordinate must be present. Cartopy is used for the coordinate
+        A valid CRS coordinate must be present. PyProj is used for the coordinate
         transformations.
 
         """
@@ -840,7 +861,7 @@ def assign_y_x(self, force=False, tolerance=None):
 
         Notes
         -----
-        A valid CRS coordinate must be present. Cartopy is used for the coordinate
+        A valid CRS coordinate must be present. PyProj is used for the coordinate
         transformations.
 
         """
@@ -994,14 +1015,16 @@ def _assign_crs(xarray_object, cf_attributes, cf_kwargs):
     attrs = cf_attributes if cf_attributes is not None else cf_kwargs
 
     # Assign crs coordinate to xarray object
-    return xarray_object.assign_coords(crs=CFProjection(attrs))
+    return xarray_object.assign_coords(metpy_crs=CFProjection(attrs))
 
 
 def _build_latitude_longitude(da):
     """Build latitude/longitude coordinates from DataArray's y/x coordinates."""
+    from pyproj import Proj
+
     y, x = da.metpy.coordinates('y', 'x')
     xx, yy = np.meshgrid(x.values, y.values)
-    lonlats = da.metpy.cartopy_geodetic.transform_points(da.metpy.cartopy_crs, xx, yy)
+    lonlats = np.stack(Proj(da.metpy.pyproj_crs)(xx, yy, inverse=True, radians=False), axis=-1)
     longitude = xr.DataArray(lonlats[..., 0], dims=(y.name, x.name),
                              coords={y.name: y, x.name: x},
                              attrs={'units': 'degrees_east', 'standard_name': 'longitude'})
@@ -1013,6 +1036,8 @@ def _build_latitude_longitude(da):
 
 def _build_y_x(da, tolerance):
     """Build y/x coordinates from DataArray's latitude/longitude coordinates."""
+    from pyproj import Proj
+
     # Initial sanity checks
     latitude, longitude = da.metpy.coordinates('latitude', 'longitude')
     if latitude.dims != longitude.dims:
@@ -1022,9 +1047,12 @@ def _build_y_x(da, tolerance):
                          'must be 2D')
 
     # Convert to projected y/x
-    xxyy = da.metpy.cartopy_crs.transform_points(da.metpy.cartopy_geodetic,
-                                                 longitude.values,
-                                                 latitude.values)
+    xxyy = np.stack(Proj(da.metpy.pyproj_crs)(
+        longitude.values,
+        latitude.values,
+        inverse=False,
+        radians=False
+    ), axis=-1)
 
     # Handle tolerance
     tolerance = 1 if tolerance is None else tolerance.m_as('m')
@@ -1051,7 +1079,7 @@ def _build_y_x(da, tolerance):
     else:
         raise ValueError('Projected y and x coordinates cannot be collapsed to 1D within '
                          'tolerance. Verify that your latitude and longitude coordinates '
-                         'correpsond to your CRS coordinate.')
+                         'correspond to your CRS coordinate.')
 
 
 def preprocess_xarray(func):

tests/calc/test_calc_tools.py

@@ -497,9 +497,11 @@ def test_lat_lon_grid_deltas_mismatched_shape():
 @needs_pyproj
 def test_lat_lon_grid_deltas_geod_kwargs():
     """Test that geod kwargs are overridden by users #774."""
+    from pyproj import Geod
+
     lat = np.arange(40, 50, 2.5)
     lon = np.arange(-100, -90, 2.5)
-    dx, dy = lat_lon_grid_deltas(lon, lat, a=4370997)
+    dx, dy = lat_lon_grid_deltas(lon, lat, geod=Geod(a=4370997))
     dx_truth = np.array([[146095.76101984, 146095.76101984, 146095.76101984],
                          [140608.9751528, 140608.9751528, 140608.9751528],
                          [134854.56713287, 134854.56713287, 134854.56713287],
@@ -1095,29 +1097,33 @@ def test_grid_deltas_from_dataarray_xy(test_da_xy):
     assert_array_almost_equal(dy, true_dy, 5)
 
 
-def test_grid_deltas_from_dataarray_actual_xy(test_da_xy, ccrs):
+@needs_pyproj
+def test_grid_deltas_from_dataarray_actual_xy(test_da_xy):
     """Test grid_deltas_from_dataarray with a xy grid and kind='actual'."""
     # Construct lon/lat coordinates
+    from pyproj import Proj
     y, x = xr.broadcast(*test_da_xy.metpy.coordinates('y', 'x'))
-    lonlat = (ccrs.Geodetic(test_da_xy.metpy.cartopy_globe)
-              .transform_points(test_da_xy.metpy.cartopy_crs, x.values, y.values))
-    lon = lonlat[..., 0]
-    lat = lonlat[..., 1]
+    lon, lat = Proj(test_da_xy.metpy.pyproj_crs)(
+        x.values,
+        y.values,
+        inverse=True,
+        radians=False
+    )
     test_da_xy = test_da_xy.assign_coords(
         longitude=xr.DataArray(lon, dims=('y', 'x'), attrs={'units': 'degrees_east'}),
         latitude=xr.DataArray(lat, dims=('y', 'x'), attrs={'units': 'degrees_north'}))
 
     # Actually test calculation
     dx, dy = grid_deltas_from_dataarray(test_da_xy, kind='actual')
-    true_dx = [[[[494426.3249766, 493977.6028005, 493044.0656467],
-                 [498740.2046073, 498474.9771064, 497891.6588559],
-                 [500276.2649627, 500256.3440237, 500139.9484845],
-                 [498740.6956936, 499045.0391707, 499542.7244501]]]] * units.m
-    true_dy = [[[[496862.4106337, 496685.4729999, 496132.0732114, 495137.8882404],
-                 [499774.9676486, 499706.3354977, 499467.5546773, 498965.2587818],
-                 [499750.8962991, 499826.2263137, 500004.4977747, 500150.9897759]]]] * units.m
-    assert_array_almost_equal(dx, true_dx, 3)
-    assert_array_almost_equal(dy, true_dy, 3)
+    true_dx = [[[[494152.626, 493704.152, 492771.132],
+                 [498464.118, 498199.037, 497616.042],
+                 [499999.328, 499979.418, 499863.087],
+                 [498464.608, 498768.783, 499266.193]]]] * units.m
+    true_dy = [[[[496587.363, 496410.523, 495857.430, 494863.795],
+                 [499498.308, 499429.714, 499191.065, 498689.047],
+                 [499474.250, 499549.538, 499727.711, 499874.122]]]] * units.m
+    assert_array_almost_equal(dx, true_dx, 2)
+    assert_array_almost_equal(dy, true_dy, 2)
 
 
 def test_grid_deltas_from_dataarray_nominal_lonlat(test_da_lonlat):
@@ -1177,7 +1183,7 @@ def test_first_derivative_xarray_lonlat(test_da_lonlat):
         coords=(('lat', test_da_lonlat['lat']),)
     )
     _, truth = xr.broadcast(test_da_lonlat, partial)
-    truth.coords['crs'] = test_da_lonlat['crs']
+    truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs']
     truth.attrs['units'] = 'kelvin / meter'
     truth = truth.metpy.quantify()
 
@@ -1231,7 +1237,7 @@ def test_second_derivative_xarray_lonlat(test_da_lonlat):
         coords=(('lat', test_da_lonlat['lat']),)
     )
     _, truth = xr.broadcast(test_da_lonlat, partial)
-    truth.coords['crs'] = test_da_lonlat['crs']
+    truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs']
     truth.attrs['units'] = 'kelvin / meter^2'
     truth = truth.metpy.quantify()
 
@@ -1259,7 +1265,7 @@ def test_gradient_xarray(test_da_xy):
         coords=(('isobaric', test_da_xy['isobaric']),)
     )
     _, truth_p = xr.broadcast(test_da_xy, partial)
-    truth_p.coords['crs'] = test_da_xy['crs']
+    truth_p.coords['metpy_crs'] = test_da_xy['metpy_crs']
     truth_p.attrs['units'] = 'kelvin / hectopascal'
     truth_p = truth_p.metpy.quantify()
 
@@ -1343,7 +1349,7 @@ def test_laplacian_xarray_lonlat(test_da_lonlat):
         coords=(('lat', test_da_lonlat['lat']),)
     )
     _, truth = xr.broadcast(test_da_lonlat, partial)
-    truth.coords['crs'] = test_da_lonlat['crs']
+    truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs']
     truth.attrs['units'] = 'kelvin / meter^2'
     truth = truth.metpy.quantify()