Update AreaDefinition equality to use pyproj CRS

pyresample/geometry.py

@@ -1817,7 +1817,7 @@ def create_areas_def_legacy(self):
     def __eq__(self, other):
         """Test for equality."""
         try:
-            return ((self.proj_str == other.proj_str) and
+            return ((self.crs == other.crs) and
                     (self.shape == other.shape) and
                     (np.allclose(self.area_extent, other.area_extent)))
         except AttributeError:

pyresample/test/test_files/areas.cfg

@@ -37,7 +37,7 @@ REGION: pc_world {
 REGION: ortho {
   NAME:    Ortho globe
   PCS_ID:  ortho_globe
-  PCS_DEF: proj=ortho, a=6370997.0, lon_0=40, lat_0=-40
+  PCS_DEF: proj=ortho, lon_0=40, lat_0=-40, ellps=sphere
   XSIZE: 640
   YSIZE: 480
   AREA_EXTENT:  (-10000000, -10000000, 10000000, 10000000)

pyresample/test/test_geometry.py

@@ -24,7 +24,6 @@
 
 import dask.array as da
 import numpy as np
-import pyproj
 import pytest
 import xarray as xr
 from pyproj import CRS
@@ -2128,18 +2127,17 @@ def test_append_area_defs(self, adef):
                                      area1.crs, area1.width, y_size, area_extent)
 
     @staticmethod
-    def _compare_area_defs(actual, expected):
-        actual_dict = actual.proj_dict
-        if 'EPSG' in actual_dict or 'init' in actual_dict:
-            # Use formal definition of EPSG projections to make them comparable to the base definition
-            proj_def = pyproj.Proj(actual_dict).definition_string().strip()
-            actual = actual.copy(projection=proj_def)
-
-            # Remove extra attributes from the formal definition
-            # for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']:
-            #     actual.proj_dict.pop(key, None)
-
-        assert actual == expected
+    def _compare_area_defs(actual, expected, use_proj4=False):
+        if use_proj4:
+            # some EPSG codes have a lot of extra metadata that makes the CRS
+            # unequal. Skip real area equality and use this as an approximation
+            actual_str = actual.crs.to_proj4()
+            expected_str = expected.crs.to_proj4()
+            assert actual_str == expected_str
+            assert actual.shape == expected.shape
+            np.allclose(actual.area_extent, expected.area_extent)
+        else:
+            assert actual == expected
 
     @pytest.mark.parametrize(
         'projection',
@@ -2206,7 +2204,7 @@ def test_create_area_def_base_combinations(self, projection, center, units):
             pytest.raises(ValueError, cad, *args, **kwargs)
         else:
             area_def = cad(*args, **kwargs)
-            self._compare_area_defs(area_def, base_def)
+            self._compare_area_defs(area_def, base_def, use_proj4="EPSG" in projection)
 
     def test_create_area_def_extra_combinations(self):
         """Test extra combinations of create_area_def parameters."""

pyresample/test/test_utils.py

@@ -492,8 +492,7 @@ def test_get_area_def_from_raster_non_georef_respects_proj_dict(self):
         transform = Affine(300.0379266750948, 0.0, 101985.0,
                            0.0, -300.041782729805, 2826915.0)
         source = tmptiff(transform=transform)
-        proj_dict = {'init': 'epsg:3857'}
-        area_def = utils.rasterio.get_area_def_from_raster(source, proj_dict=proj_dict)
+        area_def = utils.rasterio.get_area_def_from_raster(source, projection="EPSG:3857")
         self.assertEqual(area_def.crs, CRS(3857))
 
 

pyresample/utils/rasterio.py

@@ -19,7 +19,7 @@
 from . import proj4_str_to_dict
 
 
-def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, proj_dict=None):
+def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, projection=None):
     from pyresample.geometry import AreaDefinition
 
     # a: width of a pixel
@@ -33,46 +33,44 @@ def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, proj
         raise ValueError('Rotated rasters are not supported at this time.')
     area_extent = (c, f + e * dataset.RasterYSize, c + a * dataset.RasterXSize, f)
 
-    if proj_dict is None:
+    if projection is None:
         from osgeo import osr
         proj = dataset.GetProjection()
         if proj != '':
             sref = osr.SpatialReference(wkt=proj)
-            proj_dict = proj4_str_to_dict(sref.ExportToProj4())
+            projection = proj4_str_to_dict(sref.ExportToProj4())
         else:
             raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
 
         if proj_id is None:
             proj_id = proj.split('"')[1]
 
-    area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+    area_def = AreaDefinition(area_id, name, proj_id, projection,
                               dataset.RasterXSize, dataset.RasterYSize, area_extent)
     return area_def
 
 
-def _get_area_def_from_rasterio(dataset, area_id, name, proj_id=None, proj_dict=None):
+def _get_area_def_from_rasterio(dataset, area_id, name, proj_id=None, projection=None):
     from pyresample.geometry import AreaDefinition
 
     a, b, c, d, e, f, _, _, _ = dataset.transform
     if not (b == d == 0):
         raise ValueError('Rotated rasters are not supported at this time.')
 
-    if proj_dict is None:
-        crs = dataset.crs
-        if crs is not None:
-            proj_dict = dataset.crs.to_dict()
-        else:
+    if projection is None:
+        projection = dataset.crs
+        if projection is None:
             raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
 
         if proj_id is None:
-            proj_id = crs.wkt.split('"')[1]
+            proj_id = projection.wkt.split('"')[1]
 
-    area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+    area_def = AreaDefinition(area_id, name, proj_id, projection,
                               dataset.width, dataset.height, dataset.bounds)
     return area_def
 
 
-def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj_dict=None):
+def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, projection=None):
     """Construct AreaDefinition object from raster.
 
     Parameters
@@ -86,8 +84,9 @@ def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj
         Name of area
     proj_id : str, optional
         ID of projection
-    proj_dict : dict, optional
-        PROJ.4 parameters
+    projection : pyproj.CRS, dict, or string, optional
+        Projection parameters as a CRS object or a dictionary or string of
+        PROJ.4 parameters.
 
     Returns
     -------
@@ -114,8 +113,8 @@ def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj
 
     try:
         if rasterio is not None and isinstance(source, (rasterio.io.DatasetReader, rasterio.io.DatasetWriter)):
-            return _get_area_def_from_rasterio(source, area_id, name, proj_id, proj_dict)
-        return _get_area_def_from_gdal(source, area_id, name, proj_id, proj_dict)
+            return _get_area_def_from_rasterio(source, area_id, name, proj_id, projection)
+        return _get_area_def_from_gdal(source, area_id, name, proj_id, projection)
     finally:
         if cleanup_rasterio:
             source.close()