Fix compatibility with cartopy 0.20.0+

pyresample/geometry.py

@@ -1606,21 +1606,28 @@ def __str__(self):
 
     def to_cartopy_crs(self):
         """Convert projection to cartopy CRS object."""
-        from pyresample.utils.cartopy import from_proj
+        import cartopy.crs as ccrs
+        if not issubclass(ccrs.Projection, CRS):
+            raise ImportError("Pyresample only supports converting to cartopy "
+                              "0.20.0+ CRS objects. Either update cartopy or "
+                              "downgrade to an older version of Pyresample "
+                              "(<1.22.0) that supports older versions of "
+                              "cartopy.")
+
+        # cartopy 0.20+ are subclasses of Pyproj CRS class
         bounds = (self.area_extent[0],
                   self.area_extent[2],
                   self.area_extent[1],
                   self.area_extent[3])
-        if self.crs.to_epsg() is not None:
-            proj_params = "EPSG:{}".format(self.crs.to_epsg())
-        else:
-            proj_params = self.crs.to_proj4()
-        if self.crs.is_geographic:
-            # Convert area extent from degrees to radians
-            bounds = np.deg2rad(bounds)
-        crs = from_proj(proj_params, bounds=bounds)
+        from pyresample.utils.cartopy import Projection
+        crs = Projection(self.crs, bounds=bounds)
         return crs
 
+    def _cartopy_proj_params(self):
+        if self.crs.to_epsg() is not None:
+            return "EPSG:{}".format(self.crs.to_epsg())
+        return self.crs.to_proj4()
+
     def create_areas_def(self):
         """Generate YAML formatted representation of this area."""
         warnings.warn("'create_areas_def' is deprecated. Please use `dump` instead, which also "

pyresample/test/test_geometry.py

@@ -137,12 +137,7 @@ def test_cartopy_crs(self):
                 projection=projection,
                 width=123, height=123,
                 area_extent=[-40000., -40000., 40000., 40000.])
-            with patch('pyresample.utils.cartopy.warnings.warn') as warn:
-                # Test that user warning has been issued (EPSG to proj4 string is potentially lossy)
-                area.to_cartopy_crs()
-                if projection.startswith('EPSG'):
-                    # we'll only get this for the new EPSG:XXXX syntax
-                    warn.assert_called()
+            area.to_cartopy_crs()
 
         # Bounds for latlong projection must be specified in radians
         latlong_crs = geometry.AreaDefinition(area_id='latlong',
@@ -152,7 +147,7 @@ def test_cartopy_crs(self):
                                               width=360,
                                               height=180,
                                               area_extent=(-180, -90, 180, 90)).to_cartopy_crs()
-        self.assertTrue(np.allclose(latlong_crs.bounds, [-np.pi, np.pi, -np.pi/2, np.pi/2]))
+        np.testing.assert_allclose(latlong_crs.bounds, [-180, 180, -90, 90])
 
     def test_dump(self):
         """Test exporting area defs."""

pyresample/utils/cartopy.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-# Copyright (C) 2018 PyTroll developers
+# Copyright (C) 2018-2021 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify it under
 # the terms of the GNU Lesser General Public License as published by the Free
@@ -14,101 +14,58 @@
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""Classes for geometry operations.
 
-"""Classes for geometry operations."""
+See the related Cartopy pull request that implements this functionality
+directly in Cartopy: https://github.com/SciTools/cartopy/pull/1888
+
+"""
+from __future__ import annotations
 
-from logging import getLogger
 import numpy as np
 import pyproj
-import warnings
 
 try:
     from xarray import DataArray
 except ImportError:
     DataArray = np.ndarray
 
-from pyresample.utils.proj4 import proj4_str_to_dict
 import cartopy.crs as ccrs
-import shapely.geometry as sgeom
-
-try:
-    from cartopy.crs import from_proj
-except ImportError:
-    from_proj = None
-
-logger = getLogger(__name__)
-
-_GLOBE_PARAMS = {'datum': 'datum',
-                 'ellps': 'ellipse',
-                 'a': 'semimajor_axis',
-                 'b': 'semiminor_axis',
-                 'f': 'flattening',
-                 'rf': 'inverse_flattening',
-                 'towgs84': 'towgs84',
-                 'nadgrids': 'nadgrids'}
-
-
-def _globe_from_proj4(proj4_terms):
-    """Create a `Globe` object from PROJ.4 parameters."""
-    globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS,
-                         proj4_terms.items())
-    globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in
-                          globe_terms})
-    return globe
-
-
-# copy of class in cartopy (before it was released)
-class _PROJ4Projection(ccrs.Projection):
 
-    def __init__(self, proj4_terms, globe=None, bounds=None):
-        if 'EPSG' in proj4_terms.upper():
-            warnings.warn('Converting EPSG projection to proj4 string, which is a potentially lossy transformation')
-            proj4_terms = pyproj.Proj(proj4_terms).definition_string().strip()
-        terms = proj4_str_to_dict(proj4_terms)
-        globe = _globe_from_proj4(terms) if globe is None else globe
 
-        other_terms = []
-        for term in terms.items():
-            if term[0] not in _GLOBE_PARAMS:
-                other_terms.append(term)
-        super(_PROJ4Projection, self).__init__(other_terms, globe)
+class Projection(ccrs.Projection):
+    """Flexible generic Projection with optionally specified bounds."""
+
+    def __init__(self,
+                 crs: pyproj.CRS,
+                 bounds: list[float, float, float, float] = None,
+                 transform_bounds: bool = False):
+        """Initialize CRS instance and compute bounds if possible."""
+        # NOTE: Skip the base cartopy Projection class __init__
+        super(ccrs.Projection, self).__init__(crs)
+        if bounds is None and self.area_of_use is not None:
+            bounds = (
+                self.area_of_use.west,
+                self.area_of_use.east,
+                self.area_of_use.south,
+                self.area_of_use.north,
+            )
+            transform_bounds = True
 
         self.bounds = bounds
-
-    def __repr__(self):
-        return '_PROJ4Projection({})'.format(self.proj4_init)
-
-    @property
-    def boundary(self):
-        x0, x1, y0, y1 = self.bounds
-        return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0),
-                                 (x0, y0)])
-
-    @property
-    def x_limits(self):
-        x0, x1, y0, y1 = self.bounds
-        return (x0, x1)
-
-    @property
-    def y_limits(self):
-        x0, x1, y0, y1 = self.bounds
-        return (y0, y1)
-
-    @property
-    def threshold(self):
-        x0, x1, y0, y1 = self.bounds
-        return min(abs(x1 - x0), abs(y1 - y0)) / 100.
-
-
-def _lesser_from_proj(proj4_terms, globe=None, bounds=None):
-    """Not-as-good version of cartopy's 'from_proj' function.
-
-    The user doesn't have a newer version of Cartopy so there is no
-    `from_proj` function to use which does a fancier job of creating CRS
-    objects from PROJ.4 strings than this does.
-    """
-    return _PROJ4Projection(proj4_terms, globe=globe, bounds=bounds)
-
-
-if from_proj is None:
-    from_proj = _lesser_from_proj
+        if bounds is not None and transform_bounds:
+            # Convert lat/lon bounds to projected bounds.
+            # Geographic area of the entire dataset referenced to WGS 84
+            # NB. We can't use a polygon transform at this stage because
+            # that relies on the existence of the map boundary... the very
+            # thing we're trying to work out! ;-)
+            x0, x1, y0, y1 = bounds
+            lons = np.array([x0, x0, x1, x1])
+            lats = np.array([y0, y1, y1, y0])
+            points = self.transform_points(self.as_geodetic(), lons, lats)
+            x = points[:, 0]
+            y = points[:, 1]
+            self.bounds = (x.min(), x.max(), y.min(), y.max())
+        if self.bounds is not None:
+            x0, x1, y0, y1 = self.bounds
+            self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.

setup.py

@@ -30,15 +30,15 @@
                 'pykdtree>=1.3.1', 'pyyaml', 'numpy>=1.10.0',
                 ]
 extras_require = {'numexpr': ['numexpr'],
-                  'quicklook': ['matplotlib', 'cartopy', 'pillow'],
+                  'quicklook': ['matplotlib', 'cartopy>=0.20.0', 'pillow'],
                   'rasterio': ['rasterio'],
                   'dask': ['dask>=0.16.1'],
                   'cf': ['xarray'],
                   'gradient_search': ['shapely'],
                   'xarray_bilinear': ['xarray', 'dask', 'zarr']}
 
 setup_requires = ['numpy>=1.10.0', 'cython']
-test_requires = ['rasterio', 'dask', 'xarray', 'cartopy', 'pillow', 'matplotlib', 'scipy', 'zarr']
+test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr']
 
 if sys.platform.startswith("win"):
     extra_compile_args = []