CHANGELOG.md

+## Version 1.12.3 (2019/05/17)
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 193](https://github.com/pytroll/pyresample/pull/193) - Fix striding slicing in AreaDefinition
+
+In this release 1 pull request was closed.
+
+
+## Version 1.12.2 (2019/05/10)
+
+### Issues Closed
+
+* [Issue 187](https://github.com/pytroll/pyresample/issues/187) - Numerous `RuntimeWarning`s when resampling
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 190](https://github.com/pytroll/pyresample/pull/190) - Fix aggregate method using non-serializable internal function
+* [PR 189](https://github.com/pytroll/pyresample/pull/189) - Fix dask race condition in KDTree resampling
+
+#### Features added
+
+* [PR 183](https://github.com/pytroll/pyresample/pull/183) - Fix bb computation to generate areas with equal h and v resolutions
+
+In this release 3 pull requests were closed.
+
+
+## Version 1.12.1 (2019/04/24)
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 186](https://github.com/pytroll/pyresample/pull/186) - Fix support for pyproj-2 EPSG syntax
+
+#### Documentation changes
+
+* [PR 185](https://github.com/pytroll/pyresample/pull/185) - Fix argument order in get_area_def doc
+
+In this release 2 pull requests were closed.
+
+
 ## Version 1.12.0 (2019/04/06)
 
 ### Issues Closed
@@ -15,8 +63,9 @@ In this release 1 issue was closed.
 #### Features added
 
 * [PR 182](https://github.com/pytroll/pyresample/pull/182) - Implement striding and aggregation for Swath- and AreaDefinitions
+* [PR 180](https://github.com/pytroll/pyresample/pull/180) - Remove radians from create_area_def and allow compatibility with pyproj-2.0+
 
-In this release 2 pull requests were closed.
+In this release 3 pull requests were closed.
 
 
 ## Version 1.11.2 (2019/03/18)

debian/changelog

-pyresample (1.12.0-2) UNRELEASED; urgency=medium
+pyresample (1.12.3-1) experimental; urgency=medium
 
-  * Team upload.
+  [ Bas Couwenberg ]
   * Update gbp.conf to use --source-only-changes by default.
 
- -- Bas Couwenberg <sebastic@debian.org>  Sun, 07 Jul 2019 09:41:46 +0200
+  [ Antonio Valentino ]
+  * New upstream release.
+  * debian/patches:
+    - refrsh all patches
+  * Set distribution to unstable.
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Wed, 10 Jul 2019 06:06:22 +0000
 
 pyresample (1.12.0-1) experimental; urgency=medium
 

debian/patches/0001-fix-proj4-initialization.patch

@@ -21,10 +21,10 @@ index 063264d..620e665 100644
    YSIZE: 480
    AREA_EXTENT:  (-20037508.342789244, -10018754.171394622, 20037508.342789244, 10018754.171394622)
 diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
-index b916fa6..98ca468 100644
+index 267232c..a3f99ea 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
-@@ -450,7 +450,7 @@ class Test(unittest.TestCase):
+@@ -548,7 +548,7 @@ class Test(unittest.TestCase):
          swath_def = geometry.SwathDefinition(lons, lats)
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -33,7 +33,7 @@ index b916fa6..98ca468 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -475,7 +475,7 @@ class Test(unittest.TestCase):
+@@ -573,7 +573,7 @@ class Test(unittest.TestCase):
          data = np.array([1, 2, 3, 4])
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -42,7 +42,7 @@ index b916fa6..98ca468 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -510,7 +510,7 @@ class Test(unittest.TestCase):
+@@ -608,7 +608,7 @@ class Test(unittest.TestCase):
          data = np.dstack((data1, data2, data3))
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -51,7 +51,7 @@ index b916fa6..98ca468 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -579,7 +579,7 @@ class Test(unittest.TestCase):
+@@ -677,7 +677,7 @@ class Test(unittest.TestCase):
  
      def test_latlong_area(self):
          area_def = geometry.AreaDefinition('', '', '',

debian/patches/0002-Skip-dask-related-tests-if-dask-is-not-available.patch

@@ -9,11 +9,11 @@ Subject: Skip dask-related tests if dask is not available
  3 files changed, 18 insertions(+), 2 deletions(-)
 
 diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
-index 98ca468..8fbe1ed 100644
+index a3f99ea..931a13f 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
 @@ -9,7 +9,8 @@ import numpy as np
- from pyresample import geo_filter, geometry
+ from pyresample import geo_filter, geometry, parse_area_file
  from pyresample.geometry import (IncompatibleAreas,
                                   combine_area_extents_vertical,
 -                                 concatenate_area_defs)
@@ -34,7 +34,7 @@ index 98ca468..8fbe1ed 100644
  
  class Test(unittest.TestCase):
  
-@@ -725,6 +731,7 @@ class Test(unittest.TestCase):
+@@ -823,6 +829,7 @@ class Test(unittest.TestCase):
                                     np.array([-675286.976033, -682358.043845, -689429.111657, -696500.179469,
                                               -703571.247281]))
  
@@ -42,16 +42,16 @@ index 98ca468..8fbe1ed 100644
      def test_get_proj_coords_dask(self):
          """Test get_proj_coords usage with dask arrays."""
          from pyresample import utils
-@@ -1267,6 +1274,8 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1389,6 +1396,8 @@ class TestSwathDefinition(unittest.TestCase):
          assert_np_dict_allclose(res.proj_dict, proj_dict)
-         self.assertEqual(res.shape, (3, 3))
+         self.assertEqual(res.shape, (6, 3))
  
 +    @unittest.skipIf(not hasattr(DataArray, 'coarsen'), 'DataArray.coarsen not available')
 +    @unittest.skipIf(not dask, 'dask not available')
      def test_aggregation(self):
          """Test aggregation on SwathDefinitions."""
          if (sys.version_info < (3, 0)):
-@@ -1284,6 +1293,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1406,6 +1415,7 @@ class TestSwathDefinition(unittest.TestCase):
          np.testing.assert_allclose(res.lons, [[179, -179]])
          np.testing.assert_allclose(res.lats, [[0.5, 0.5]], atol=2e-5)
  

docs/source/geo_def.rst

@@ -51,11 +51,13 @@ where
 * **upper_right_x**: projection x coordinate of upper right corner of upper right pixel
 * **upper_right_y**: projection y coordinate of upper right corner of upper right pixel
 
-Below is an example of creating an ``AreaDefinition``:
+Below are three examples of creating an ``AreaDefinition``:
 
 .. doctest::
 
  >>> from pyresample.geometry import AreaDefinition
+
+ >>> # a) Using a projection dictionary
  >>> area_id = 'ease_sh'
  >>> description = 'Antarctic EASE grid'
  >>> proj_id = 'ease_sh'
@@ -74,6 +76,37 @@ Below is an example of creating an ``AreaDefinition``:
  Number of rows: 425
  Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
 
+ >>> # b) Using an explicit proj4 string
+ >>> proj_string = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
+ >>> area_def = AreaDefinition(area_id, description, proj_id, proj_string,
+ ...                           width, height, area_extent)
+ >>> print(area_def)
+ Area ID: ease_sh
+ Description: Antarctic EASE grid
+ Projection ID: ease_sh
+ Projection: {'a': '6371228.0', 'lat_0': '-90.0', 'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}
+ Number of columns: 425
+ Number of rows: 425
+ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+
+ >>> # c) Using an EPSG code in a proj4 string
+ >>> proj_string = '+init=EPSG:3409'  # Use 'EPSG:3409' with pyproj 2.0+
+ >>> area_def = AreaDefinition(area_id, description, proj_id, proj_string,
+ ...                           width, height, area_extent)
+ >>> print(area_def)
+ Area ID: ease_sh
+ Description: Antarctic EASE grid
+ Projection ID: ease_sh
+ Projection: {'init': 'EPSG:3409'}
+ Number of columns: 425
+ Number of rows: 425
+ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+
+.. note::
+
+  When using pyproj 2.0+, please use the new ``'EPSG:XXXX'`` syntax
+  as the old ``'+init=EPSG:XXXX'`` is no longer supported.
+
 Creating an ``AreaDefinition`` can be complex if you don't know everything
 about the region being described. Pyresample provides multiple utilities
 for creating areas as well as storing them on disk for repeated use. See

docs/source/geometry_utils.rst

@@ -212,10 +212,10 @@ from_ul_corner
 Loading from disk
 -----------------
 
-The :func:`~pyresample.utils.load_area` function can be used to
+The :func:`~pyresample.area_config.load_area` function can be used to
 parse area definitions from a configuration file by giving it the
-area file name and regions you wish to load. :func:`~pyresample.utils.load_area`
-takes advantage of :func:`~pyresample.utils.create_area_def`
+area file name and regions you wish to load. :func:`~pyresample.area_config.load_area`
+takes advantage of :func:`~pyresample.area_config.create_area_def`
 and hence allows for the same arguments in the on-disk file.
 Pyresample uses the YAML file format to store on-disk area definitions.
 Below is an example YAML configuration file showing the various ways
@@ -344,6 +344,14 @@ an area might be specified.
      resolution: 0.22542974631297721
      units: deg
 
+ epsg:
+   area_id: ease_sh
+   description: Example of making an area definition using EPSG codes
+   projection:
+     init: EPSG:3410
+   shape: [425, 425]
+   area_extent: [-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625]
+
 .. note::
 
   The `lower_left_xy` and `upper_right_xy` items give the coordinates of the
@@ -351,13 +359,18 @@ an area might be specified.
   projection coordinates are longitudes and latitudes, it is expected to
   provide the extent in `longitude, latitude` order.
 
+.. note::
+
+  When using pyproj 2.0+, please use the new ``'EPSG: XXXX'`` syntax
+  as the old ``'init: EPSG:XXXX'`` is no longer supported.
+
 If we assume the YAML content is stored in an ``areas.yaml`` file, we can
 read a single ``AreaDefinition`` named ``corner`` by doing:
 
 .. doctest::
 
- >>> from pyresample import utils
- >>> area_def = utils.load_area('areas.yaml', 'corner')
+ >>> from pyresample import load_area
+ >>> area_def = load_area('areas.yaml', 'corner')
  >>> print(area_def)
  Area ID: corner
  Description: Example of making an area definition using shape, upper_left_extent, and resolution
@@ -371,7 +384,7 @@ series of arguments:
 
 .. doctest::
 
- >>> corner, boundary = utils.load_area('areas.yaml', 'corner', 'boundary')
+ >>> corner, boundary = load_area('areas.yaml', 'corner', 'boundary')
  >>> print(boundary)
  Area ID: ease_sh
  Description: Example of making an area definition using shape and area_extent
@@ -410,8 +423,8 @@ An area definition dict can be read using
 
 .. doctest::
 
- >>> from pyresample import utils
- >>> area = utils.load_area('areas.cfg', 'ease_nh')
+ >>> from pyresample import load_area
+ >>> area = load_area('areas.cfg', 'ease_nh')
  >>> print(area)
  Area ID: ease_nh
  Description: Arctic EASE grid
@@ -427,7 +440,7 @@ Several area definitions can be read at once using the region names in an argume
 
 .. doctest::
 
- >>> nh_def, sh_def = utils.load_area('areas.cfg', 'ease_nh', 'ease_sh')
+ >>> nh_def, sh_def = load_area('areas.cfg', 'ease_nh', 'ease_sh')
  >>> print(sh_def)
  Area ID: ease_sh
  Description: Antarctic EASE grid

pyresample/_cartopy.py

@@ -19,6 +19,8 @@
 
 from logging import getLogger
 import numpy as np
+import pyproj
+import warnings
 
 try:
     from xarray import DataArray
@@ -59,6 +61,9 @@ def _globe_from_proj4(proj4_terms):
 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
 

pyresample/_spatial_mp.py

@@ -24,7 +24,7 @@ import pyproj
 import multiprocessing as mp
 import warnings
 
-from pyresample.utils import proj4_str_to_dict
+from pyresample.utils import proj4_str_to_dict, is_pyproj2
 
 
 try:
@@ -112,8 +112,12 @@ class cKDTree_MP(object):
 class BaseProj(pyproj.Proj):
 
     def __init__(self, projparams=None, preserve_units=True, **kwargs):
+        # Copy dict-type arguments as they will be modified in-place
+        if isinstance(projparams, dict):
+            projparams = projparams.copy()
+
         # Pyproj<2 uses __new__ to initiate data and does not define its own __init__ method.
-        if pyproj.__version__ >= '2':
+        if is_pyproj2():
             # If init is found in any of the data, override any other area parameters.
             if 'init' in kwargs:
                 warnings.warn('init="EPSG:XXXX" is no longer supported. Use "EPSG:XXXX" as a proj string instead')
@@ -127,10 +131,18 @@ class BaseProj(pyproj.Proj):
                     projparams = proj4_str_to_dict(projparams)
                 warnings.warn(warn_msg)
                 projparams = projparams.pop('init')
+            # New syntax 'EPSG:XXXX'
+            if 'EPSG' in kwargs or (isinstance(projparams, dict) and 'EPSG' in projparams):
+                if 'EPSG' in kwargs:
+                    epsg_code = kwargs.pop('EPSG')
+                else:
+                    epsg_code = projparams.pop('EPSG')
+                projparams = 'EPSG:{}'.format(epsg_code)
+
             super(BaseProj, self).__init__(projparams=projparams, preserve_units=preserve_units, **kwargs)
 
     def is_latlong(self):
-        if pyproj.__version__ >= '2':
+        if is_pyproj2():
             return self.crs.is_geographic
         return super(BaseProj, self).is_latlong()
 

pyresample/area_config.py

@@ -305,10 +305,10 @@ def get_area_def(area_id, area_name, proj_id, proj4_args, width, height, area_ex
     -----------
     area_id : str
         ID of area
-    proj_id : str
-        ID of projection
     area_name :str
         Description of area
+    proj_id : str
+        ID of projection
     proj4_args : list, dict, or str
         Proj4 arguments as list of arguments or string
     width : int
@@ -453,9 +453,12 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No
 
     # Fills in missing information to attempt to create an area definition.
     if area_extent is None or shape is None:
-        area_extent, shape = _extrapolate_information(area_extent, shape, center, radius, resolution,
-                                                      upper_left_extent, units, p, proj_dict)
-    return _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **kwargs)
+        area_extent, shape, resolution = \
+            _extrapolate_information(area_extent, shape, center, radius,
+                                     resolution, upper_left_extent, units,
+                                     p, proj_dict)
+    return _make_area(area_id, description, proj_id, proj_dict, shape,
+                      area_extent, resolution=resolution, **kwargs)
 
 
 def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **kwargs):
@@ -464,6 +467,7 @@ def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **k
     from pyresample.geometry import DynamicAreaDefinition
     # Remove arguments that are only for DynamicAreaDefinition.
     optimize_projection = kwargs.pop('optimize_projection', False)
+    resolution = kwargs.pop('resolution', None)
     # If enough data is provided, create an AreaDefinition. If only shape or area_extent are found, make a
     # DynamicAreaDefinition. If not enough information was provided, raise a ValueError.
     height, width = (None, None)
@@ -474,7 +478,7 @@ def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **k
 
     return DynamicAreaDefinition(area_id=area_id, description=description, projection=proj_dict, width=width,
                                  height=height, area_extent=area_extent, rotation=kwargs.get('rotation'),
-                                 optimize_projection=optimize_projection)
+                                 resolution=resolution, optimize_projection=optimize_projection)
 
 
 def _get_proj_data(projection):
@@ -523,7 +527,7 @@ def _distance_from_center_forward(var, center, p):
     # Since the distance between longitudes and latitudes is not constant in
     # most projections, there must be reference point to start from.
     if center is None:
-        raise ValueError('center must be given to convert radius or resolution from an angle to meters')
+        center = (0, 0)
 
     center_as_angle = p(*center, inverse=True, errcheck=True)
     pole = 90
@@ -700,7 +704,7 @@ def _extrapolate_information(area_extent, shape, center, radius, resolution, upp
             upper_left_extent[0], upper_left_extent[1] - 2 * radius[1], upper_left_extent[0] + 2 * radius[0],
             upper_left_extent[1])
         area_extent = _validate_variable(area_extent, new_area_extent, 'area_extent', ['upper_left_extent', 'radius'])
-    return area_extent, shape
+    return area_extent, shape, resolution
 
 
 def _format_list(var, name):

pyresample/geometry.py

@@ -31,7 +31,7 @@ from logging import getLogger
 
 import numpy as np
 import yaml
-from pyproj import Geod
+from pyproj import Geod, transform
 
 from pyresample import CHUNK_SIZE, utils
 from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP
@@ -551,6 +551,12 @@ class SwathDefinition(CoordinateDefinition):
         """Copy the current swath."""
         return SwathDefinition(self.lons, self.lats)
 
+    @staticmethod
+    def _do_transform(src, dst, lons, lats, alt):
+        """Helper for 'aggregate' method."""
+        x, y, z = transform(src, dst, lons, lats, alt)
+        return np.dstack((x, y, z))
+
     def aggregate(self, **dims):
         """Aggregate the current swath definition by averaging.
 
@@ -560,20 +566,15 @@ class SwathDefinition(CoordinateDefinition):
         import pyproj
         import dask.array as da
 
-        def do_transform(src, dst, lons, lats, alt):
-            import pyproj
-            x, y, z = pyproj.transform(src, dst, lons, lats, alt)
-            return np.dstack((x, y, z))
-
         geocent = pyproj.Proj(proj='geocent')
         latlong = pyproj.Proj(proj='latlong')
-        res = da.map_blocks(do_transform, latlong, geocent,
+        res = da.map_blocks(self._do_transform, latlong, geocent,
                             self.lons.data, self.lats.data,
                             da.zeros_like(self.lons.data), new_axis=[2],
                             chunks=(self.lons.chunks[0], self.lons.chunks[1], 3))
         res = DataArray(res, dims=['y', 'x', 'coord'], coords=self.lons.coords)
         res = res.coarsen(**dims).mean()
-        lonlatalt = da.map_blocks(do_transform, geocent, latlong,
+        lonlatalt = da.map_blocks(self._do_transform, geocent, latlong,
                                   res[:, :, 0].data, res[:, :, 1].data,
                                   res[:, :, 2].data, new_axis=[2],
                                   chunks=res.data.chunks)
@@ -676,6 +677,33 @@ class SwathDefinition(CoordinateDefinition):
             new_proj.update(proj_dict)
             return new_proj
 
+    def _compute_uniform_shape(self):
+        """Compute the height and width of a domain to have uniform resolution across dimensions."""
+        g = Geod(ellps='WGS84')
+        leftlons = self.lons[:, 0]
+        leftlons = leftlons.where(leftlons.notnull(), drop=True)
+        rightlons = self.lons[:, -1]
+        rightlons = rightlons.where(rightlons.notnull(), drop=True)
+        middlelons = self.lons[:, int(self.lons.shape[1] / 2)]
+        middlelons = middlelons.where(middlelons.notnull(), drop=True)
+        leftlats = self.lats[:, 0]
+        leftlats = leftlats.where(leftlats.notnull(), drop=True)
+        rightlats = self.lats[:, -1]
+        rightlats = rightlats.where(rightlats.notnull(), drop=True)
+        middlelats = self.lats[:, int(self.lats.shape[1] / 2)]
+        middlelats = middlelats.where(middlelats.notnull(), drop=True)
+
+        az1, az2, width1 = g.inv(leftlons[0], leftlats[0], rightlons[0], rightlats[0])
+        az1, az2, width2 = g.inv(leftlons[-1], leftlats[-1], rightlons[-1], rightlats[-1])
+        az1, az2, height = g.inv(middlelons[0], middlelats[0], middlelons[-1], middlelats[-1])
+        width = min(width1, width2)
+        vresolution = height * 1.0 / self.lons.shape[0]
+        hresolution = width * 1.0 / self.lons.shape[1]
+        resolution = min(vresolution, hresolution)
+        width = int(width * 1.1 / resolution)
+        height = int(height * 1.1 / resolution)
+        return height, width
+
     def compute_optimal_bb_area(self, proj_dict=None):
         """Compute the "best" bounding box area for this swath with `proj_dict`.
 
@@ -683,16 +711,16 @@ class SwathDefinition(CoordinateDefinition):
         which case the right projection angle `alpha` is computed from the
         swath centerline. For other projections, only the appropriate center of
         projection and area extents are computed.
+
+        The height and width are computed so that the resolution is
+        approximately the same across dimensions.
         """
         if proj_dict is None:
             proj_dict = {}
         projection = proj_dict.setdefault('proj', 'omerc')
         area_id = projection + '_otf'
         description = 'On-the-fly ' + projection + ' area'
-        lines, cols = self.lons.shape
-        width = int(cols * 1.1)
-        height = int(lines * 1.1)
-
+        height, width = self._compute_uniform_shape()
         proj_dict = self.compute_bb_proj_params(proj_dict)
 
         area = DynamicAreaDefinition(area_id, description, proj_dict)
@@ -768,10 +796,10 @@ class DynamicAreaDefinition(object):
                 x_resolution, y_resolution = resolution
             except TypeError:
                 x_resolution = y_resolution = resolution
-            width = int(np.rint((corners[2] - corners[0]) * 1.0 /
-                                 x_resolution + 1))
-            height = int(np.rint((corners[3] - corners[1]) * 1.0 /
-                                 y_resolution + 1))
+            width = int(np.rint((corners[2] - corners[0]) * 1.0
+                                / x_resolution + 1))
+            height = int(np.rint((corners[3] - corners[1]) * 1.0
+                                 / y_resolution + 1))
 
         area_extent = (corners[0] - x_resolution / 2,
                        corners[1] - y_resolution / 2,
@@ -1712,8 +1740,10 @@ class AreaDefinition(BaseDefinition):
         yslice, xslice = key
         # Get actual values, replace Nones
         yindices = yslice.indices(self.height)
+        total_rows = int((yindices[1] - yindices[0]) / yindices[2])
         ystopactual = yindices[1] - (yindices[1] - 1) % yindices[2]
         xindices = xslice.indices(self.width)
+        total_cols = int((xindices[1] - xindices[0]) / xindices[2])
         xstopactual = xindices[1] - (xindices[1] - 1) % xindices[2]
         yslice = slice(yindices[0], ystopactual, yindices[2])
         xslice = slice(xindices[0], xstopactual, xindices[2])
@@ -1725,8 +1755,8 @@ class AreaDefinition(BaseDefinition):
 
         new_area = AreaDefinition(self.area_id, self.description,
                                   self.proj_id, self.proj_dict,
-                                  int((xslice.stop - xslice.start) / xslice.step),
-                                  int((yslice.stop - yslice.start) / yslice.step),
+                                  total_cols,
+                                  total_rows,
                                   new_area_extent)
         new_area.crop_offset = (self.crop_offset[0] + yslice.start,
                                 self.crop_offset[1] + xslice.start)

pyresample/kd_tree.py

@@ -892,13 +892,13 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
 
 
 def lonlat2xyz(lons, lats):
-
     R = 6370997.0
-    x_coords = R * da.cos(da.deg2rad(lats)) * da.cos(da.deg2rad(lons))
-    y_coords = R * da.cos(da.deg2rad(lats)) * da.sin(da.deg2rad(lons))
-    z_coords = R * da.sin(da.deg2rad(lats))
+    x_coords = R * np.cos(np.deg2rad(lats)) * np.cos(np.deg2rad(lons))
+    y_coords = R * np.cos(np.deg2rad(lats)) * np.sin(np.deg2rad(lons))
+    z_coords = R * np.sin(np.deg2rad(lats))
 
-    return da.stack(
+    stack = np.stack if isinstance(lons, np.ndarray) else da.stack
+    return stack(
         (x_coords.ravel(), y_coords.ravel(), z_coords.ravel()), axis=-1)
 
 
@@ -922,7 +922,7 @@ def query_no_distance(target_lons, target_lats, valid_output_index,
 
     coords = lonlat2xyz(target_lons_valid, target_lats_valid)
     distance_array, index_array = kdtree.query(
-        coords.compute(),
+        coords,
         k=neighbours,
         eps=epsilon,
         distance_upper_bound=radius,
@@ -946,6 +946,20 @@ def query_no_distance(target_lons, target_lats, valid_output_index,
     return res_ia
 
 
+def _my_index(index_arr, vii, data_arr, vii_slices=None, ia_slices=None,
+              fill_value=np.nan):
+    """Helper function for 'get_sample_from_neighbour_info'."""
+    vii_slices = tuple(
+        x if x is not None else vii.ravel() for x in vii_slices)
+    mask_slices = tuple(
+        x if x is not None else (index_arr == -1) for x in ia_slices)
+    ia_slices = tuple(
+        x if x is not None else index_arr for x in ia_slices)
+    res = data_arr[vii_slices][ia_slices]
+    res[mask_slices] = fill_value
+    return res
+
+
 class XArrayResamplerNN(object):
     def __init__(self,
                  source_geo_def,
@@ -1176,18 +1190,6 @@ class XArrayResamplerNN(object):
         # FUTURE: when we allow more than one neighbor
         # neighbors_dim = i + 3
 
-        def _my_index(index_arr, vii, data_arr, vii_slices=None,
-                      ia_slices=None, fill_value=np.nan):
-            vii_slices = tuple(
-                x if x is not None else vii.ravel() for x in vii_slices)
-            mask_slices = tuple(
-                x if x is not None else (index_arr == -1) for x in ia_slices)
-            ia_slices = tuple(
-                x if x is not None else index_arr for x in ia_slices)
-            res = data_arr[vii_slices][ia_slices]
-            res[mask_slices] = fill_value
-            return res
-
         new_data = data.data.reshape(flat_src_shape)
         vii = vii.ravel()
         dst_adims = [dst_dim_to_ind[dim] for dim in dst_dims]

pyresample/test/__init__.py

@@ -37,6 +37,7 @@ from pyresample.test import (
     test_ewa_fornav,
     test_bilinear,
     test_data_reduce,
+    test_spatial_mp
 )
 
 import unittest
@@ -59,6 +60,7 @@ def suite():
     mysuite.addTests(test_ewa_fornav.suite())
     mysuite.addTests(test_bilinear.suite())
     mysuite.addTests(test_data_reduce.suite())
+    mysuite.addTests(test_spatial_mp.suite())
     return mysuite
 
 

pyresample/test/test_files/areas.yaml

@@ -219,3 +219,25 @@ test_latlong:
   area_extent:
     lower_left_xy: [-0.08115781021773638, 0.4038691889114878]
     upper_right_xy: [0.08115781021773638, 0.5427973973702365]
+
+test_dynamic_resolution:
+  description: Dynamic with resolution specified in meters
+  projection:
+    proj: lcc
+    lon_0: -95.0
+    lat_0: 25.0
+    lat_1: 25.0
+    ellps: WGS84
+  resolution: [1000, 1000]
+
+test_dynamic_resolution_ll:
+  description: Dynamic with resolution specified in degrees
+  projection:
+      proj: longlat
+      lat_0: 27.12
+      lon_0: -81.36
+      ellps: WGS84
+  resolution:
+    dy: 1
+    dx: 1
+    units: deg

pyresample/test/test_geometry.py

@@ -6,7 +6,7 @@ import sys
 
 import numpy as np
 
-from pyresample import geo_filter, geometry
+from pyresample import geo_filter, geometry, parse_area_file
 from pyresample.geometry import (IncompatibleAreas,
                                  combine_area_extents_vertical,
                                  concatenate_area_defs)
@@ -70,6 +70,8 @@ class Test(unittest.TestCase):
 
     def test_cartopy_crs(self):
         """Test conversion from area definition to cartopy crs"""
+        from pyresample import utils
+
         europe = geometry.AreaDefinition(area_id='areaD',
                                          description='Europe (3km, HRV, VTC)',
                                          proj_id='areaD',
@@ -111,7 +113,28 @@ class Test(unittest.TestCase):
                              np.fabs(area_def.area_extent[3] - area_def.area_extent[1])) / 100.
             self.assertEqual(crs.threshold, thresh_exp)
 
+        # EPSG projection
+        projections = ['+init=EPSG:6932']
+        if utils.is_pyproj2():
+            projections.append('EPSG:6932')
+
+        for projection in projections:
+            area = geometry.AreaDefinition(
+                area_id='ease-sh-2.0',
+                description='25km EASE Grid 2.0 (Southern Hemisphere)',
+                proj_id='ease-sh-2.0',
+                projection=projection,
+                width=123, height=123,
+                area_extent=[-40000., -40000., 40000., 40000.])
+            with patch('pyresample._cartopy.warnings.warn') as warn:
+                # Test that user warning has been issued (EPSG to proj4 string is potentially lossy)
+                area.to_cartopy_crs()
+                warn.assert_called()
+
     def test_create_areas_def(self):
+        from pyresample import utils
+        import yaml
+
         area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
                                            'areaD',
                                            {'a': '6378144.0',
@@ -126,9 +149,8 @@ class Test(unittest.TestCase):
                                             -909968.64000000001,
                                             1029087.28,
                                             1490031.3600000001])
-        import yaml
-        res = yaml.load(area_def.create_areas_def())
-        expected = yaml.load(('areaD:\n  description: Europe (3km, HRV, VTC)\n'
+        res = yaml.safe_load(area_def.create_areas_def())
+        expected = yaml.safe_load(('areaD:\n  description: Europe (3km, HRV, VTC)\n'
                                    '  projection:\n    a: 6378144.0\n    b: 6356759.0\n'
                                    '    lat_0: 50.0\n    lat_ts: 50.0\n    lon_0: 8.0\n'
                                    '    proj: stere\n  shape:\n    height: 800\n'
@@ -138,6 +160,82 @@ class Test(unittest.TestCase):
 
         self.assertDictEqual(res, expected)
 
+        # EPSG
+        projections = {'+init=epsg:3006': 'init: epsg:3006'}
+        if utils.is_pyproj2():
+            projections['EPSG:3006'] = 'EPSG: 3006'
+
+        for projection, epsg_yaml in projections.items():
+            area_def = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm',
+                                               'sweref99tm',
+                                               projection,
+                                               4667,
+                                               4667,
+                                               [-49739, 5954123, 1350361, 7354223])
+            res = yaml.safe_load(area_def.create_areas_def())
+            expected = yaml.safe_load(('baws300_sweref99tm:\n'
+                                       '  description: BAWS, 300m resolution, sweref99tm\n'
+                                       '  projection:\n'
+                                       '    {epsg}\n'
+                                       '  shape:\n'
+                                       '    height: 4667\n'
+                                       '    width: 4667\n'
+                                       '  area_extent:\n'
+                                       '    lower_left_xy: [-49739, 5954123]\n'
+                                       '    upper_right_xy: [1350361, 7354223]'.format(epsg=epsg_yaml)))
+            self.assertDictEqual(res, expected)
+
+    def test_parse_area_file(self):
+        from pyresample import utils
+
+        expected = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
+                                           'areaD',
+                                           {'a': '6378144.0',
+                                            'b': '6356759.0',
+                                            'lat_0': '50.00',
+                                            'lat_ts': '50.00',
+                                            'lon_0': '8.00',
+                                            'proj': 'stere'},
+                                           800,
+                                           800,
+                                           [-1370912.72,
+                                            -909968.64000000001,
+                                            1029087.28,
+                                            1490031.3600000001])
+        yaml_str = ('areaD:\n  description: Europe (3km, HRV, VTC)\n'
+                    '  projection:\n    a: 6378144.0\n    b: 6356759.0\n'
+                    '    lat_0: 50.0\n    lat_ts: 50.0\n    lon_0: 8.0\n'
+                    '    proj: stere\n  shape:\n    height: 800\n'
+                    '    width: 800\n  area_extent:\n'
+                    '    lower_left_xy: [-1370912.72, -909968.64]\n'
+                    '    upper_right_xy: [1029087.28, 1490031.36]\n')
+        area_def = parse_area_file(yaml_str, 'areaD')[0]
+        self.assertEqual(area_def, expected)
+
+        # EPSG
+        projections = {'+init=epsg:3006': 'init: epsg:3006'}
+        if utils.is_pyproj2():
+            projections['EPSG:3006'] = 'EPSG: 3006'
+        for projection, epsg_yaml in projections.items():
+            expected = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm',
+                                               'sweref99tm',
+                                               projection,
+                                               4667,
+                                               4667,
+                                               [-49739, 5954123, 1350361, 7354223])
+            yaml_str = ('baws300_sweref99tm:\n'
+                        '  description: BAWS, 300m resolution, sweref99tm\n'
+                        '  projection:\n'
+                        '    {epsg}\n'
+                        '  shape:\n'
+                        '    height: 4667\n'
+                        '    width: 4667\n'
+                        '  area_extent:\n'
+                        '    lower_left_xy: [-49739, 5954123]\n'
+                        '    upper_right_xy: [1350361, 7354223]'.format(epsg=epsg_yaml))
+            area_def = parse_area_file(yaml_str, 'baws300_sweref99tm')[0]
+            self.assertEqual(area_def, expected)
+
     def test_base_type(self):
         lons1 = np.arange(-135., +135, 50.)
         lats = np.ones_like(lons1) * 70.
@@ -916,9 +1014,14 @@ class Test(unittest.TestCase):
         self.assertEqual(slice_y, slice(158, 515, None))
 
         # totally different area
-        area_to_cover = geometry.AreaDefinition('epsg4326', 'Global equal latitude/longitude grid for global sphere',
+        projections = [{"init": 'EPSG:4326'}]
+        if utils.is_pyproj2():
+            projections.append('EPSG:4326')
+        for projection in projections:
+            area_to_cover = geometry.AreaDefinition(
+                'epsg4326', 'Global equal latitude/longitude grid for global sphere',
                 'epsg4326',
-                                                {"init": 'EPSG:4326'},
+                projection,
                 8192,
                 4096,
                 [-180.0, -90.0, 180.0, 90.0])
@@ -964,6 +1067,8 @@ class Test(unittest.TestCase):
 
     def test_proj_str(self):
         from collections import OrderedDict
+        from pyresample import utils
+
         proj_dict = OrderedDict()
         proj_dict['proj'] = 'stere'
         proj_dict['a'] = 6378144.0
@@ -985,6 +1090,20 @@ class Test(unittest.TestCase):
         self.assertEqual(area.proj_str,
                          '+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +no_rot +proj=stere')
 
+        # EPSG
+        projections = ['+init=EPSG:6932']
+        if utils.is_pyproj2():
+            projections.append('EPSG:6932')
+        for projection in projections:
+            area = geometry.AreaDefinition(
+                area_id='ease-sh-2.0',
+                description='25km EASE Grid 2.0 (Southern Hemisphere)',
+                proj_id='ease-sh-2.0',
+                projection=projection,
+                width=123, height=123,
+                area_extent=[-40000., -40000., 40000., 40000.])
+            self.assertEqual(area.proj_str, projection)
+
     def test_striding(self):
         """Test striding AreaDefinitions."""
         from pyresample import utils
@@ -1009,6 +1128,7 @@ class Test(unittest.TestCase):
                                                               area_extent[1] + 3 * area_def.pixel_size_y,
                                                               area_extent[2] - 3 * area_def.pixel_size_x,
                                                               area_extent[3]))
+        self.assertEqual(reduced_area.shape, (928, 928))
 
 
 def assert_np_dict_allclose(dict1, dict2):
@@ -1247,25 +1367,27 @@ class TestSwathDefinition(unittest.TestCase):
 
     def test_compute_optimal_bb(self):
         """Test computing the bb area."""
+        import xarray as xr
         lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
                          [80.84000396728516, 60.74200439453125, 34.08500289916992],
                          [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
-
+        lats = xr.DataArray(lats)
         lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
                          [79.11000061035156, 7.284000396728516, -5.107000350952148],
                          [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+        lons = xr.DataArray(lons)
 
         area = geometry.SwathDefinition(lons, lats)
 
         res = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'})
 
-        np.testing.assert_allclose(res.area_extent, [-2348379.728104, 2284625.526467,
-                                                     2432121.058435, 11719235.223912])
+        np.testing.assert_allclose(res.area_extent, [-2348379.728104, 3228086.496211,
+                                                     2432121.058435, 10775774.254169])
         proj_dict = {'gamma': 0.0, 'lonc': -11.391744043133668,
                      'ellps': 'WGS84', 'proj': 'omerc',
                      'alpha': 9.185764390923012, 'lat_0': -0.2821013754097188}
         assert_np_dict_allclose(res.proj_dict, proj_dict)
-        self.assertEqual(res.shape, (3, 3))
+        self.assertEqual(res.shape, (6, 3))
 
     def test_aggregation(self):
         """Test aggregation on SwathDefinitions."""
@@ -1472,11 +1594,16 @@ class TestStackedAreaDefinition(unittest.TestCase):
         from pyresample.geometry import DynamicAreaDefinition
         from pyresample.area_config import DataArray
         from pyresample.area_config import create_area_def as cad
+        from pyresample import utils
+        import pyproj
 
         area_id = 'ease_sh'
         description = 'Antarctic EASE grid'
         projection_list = [{'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'},
-                           '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m']
+                           '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m',
+                           '+init=EPSG:3409']
+        if utils.is_pyproj2():
+            projection_list.append('EPSG:3409')
         proj_id = 'ease_sh'
         shape = (425, 850)
         upper_left_extent = (-5326849.0625, 5326849.0625)
@@ -1511,7 +1638,7 @@ class TestStackedAreaDefinition(unittest.TestCase):
                                      radius=essentials[1], description=description, units=units, rotation=45))
             except ValueError:
                 pass
-        self.assertEqual(len(area_list), 4)
+        self.assertEqual(len(area_list), 8 if utils.is_pyproj2() else 6)
 
         # Tests that specifying units through xarrays works.
         area_list.append(cad(area_id, projection_list[1], shape=shape,
@@ -1544,6 +1671,18 @@ class TestStackedAreaDefinition(unittest.TestCase):
         self.assertEqual(area_def.shape, (101, 90))
         # Checks every area definition made
         for area_def in area_list:
+            if 'EPSG' in area_def.proj_dict or 'init' in area_def.proj_dict:
+                # Use formal definition of EPSG projections to make them comparable to the base definition
+                proj_def = pyproj.Proj(area_def.proj_str).definition_string().strip()
+                area_def = area_def.copy(projection=proj_def)
+
+                # Remove extra attributes from the formal definition
+                if 'R' in area_def.proj_dict:
+                    # pyproj < 2
+                    area_def.proj_dict['a'] = area_def.proj_dict.pop('R')
+                for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']:
+                    area_def.proj_dict.pop(key, None)
+
             self.assertEqual(area_def, base_def)
 
         # Makes sure if shape or area_extent is found/given, a DynamicAreaDefinition is made.
@@ -1588,14 +1727,15 @@ class TestDynamicAreaDefinition(unittest.TestCase):
         lats = [[66, 67, 68, 69.],
                 [58, 59, 60, 61],
                 [50, 51, 52, 53]]
-        sdef = geometry.SwathDefinition(lons, lats)
+        import xarray as xr
+        sdef = geometry.SwathDefinition(xr.DataArray(lons), xr.DataArray(lats))
         result = area.freeze(sdef,
                              resolution=1000)
         np.testing.assert_allclose(result.area_extent,
-                                   [-336277.698941, 5047207.008079,
-                                    192456.651909, 8215588.023806])
+                                   [-336277.698941, 5513145.392745,
+                                    192456.651909, 7749649.63914])
         self.assertEqual(result.x_size, 4)
-        self.assertEqual(result.y_size, 3)
+        self.assertEqual(result.y_size, 18)
 
     def test_compute_domain(self):
         """Test computing size and area extent."""

pyresample/test/test_spatial_mp.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# pyresample, Resampling of remote sensing image data in python
+#
+# Copyright (C) 2014-2019 PyTroll 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
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+# details.
+#
+# 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/>.
+"""Testing the _spatial_mp module"""
+
+
+try:
+    from unittest import mock
+except ImportError:
+    # separate mock package py<3.3
+    import mock
+import pyproj
+import unittest
+
+from pyresample._spatial_mp import BaseProj
+
+
+class SpatialMPTest(unittest.TestCase):
+    @mock.patch('pyresample._spatial_mp.pyproj.Proj.__init__', return_value=None)
+    def test_base_proj_epsg(self, proj_init):
+        """Test Proj creation with EPSG codes"""
+        if pyproj.__version__ < '2':
+            return self.skipTest(reason='pyproj 2+ only')
+
+        args = [
+            [None, {'init': 'EPSG:6932'}],
+            [{'init': 'EPSG:6932'}, {}],
+            [None, {'EPSG': '6932'}],
+            [{'EPSG': '6932'}, {}]
+        ]
+        for projparams, kwargs in args:
+            BaseProj(projparams, **kwargs)
+            proj_init.assert_called_with(projparams='EPSG:6932', preserve_units=mock.ANY)
+            proj_init.reset_mock()
+
+
+def suite():
+    """The test suite.
+    """
+    loader = unittest.TestLoader()
+    mysuite = unittest.TestSuite()
+    mysuite.addTest(loader.loadTestsFromTestCase(SpatialMPTest))
+    return mysuite
+
+
+if __name__ == '__main__':
+    unittest.main()

pyresample/test/test_utils.py

@@ -136,6 +136,27 @@ Number of rows: 4058
 Area extent: (-0.0812, 0.4039, 0.0812, 0.5428)"""
         self.assertEqual(test_latlong.__str__(), latlong_str)
 
+    def test_dynamic_area_parser_yaml(self):
+        """Test YAML area parser on dynamic areas."""
+        from pyresample import parse_area_file
+        from pyresample.geometry import DynamicAreaDefinition
+        test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+        test_area = parse_area_file(test_area_file, 'test_dynamic_resolution')[0]
+
+        self.assertIsInstance(test_area, DynamicAreaDefinition)
+        self.assertTrue(hasattr(test_area, 'resolution'))
+        self.assertEqual(test_area.resolution, (1000.0, 1000.0))
+
+        # lat/lon
+        from pyresample import parse_area_file
+        from pyresample.geometry import DynamicAreaDefinition
+        test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+        test_area = parse_area_file(test_area_file, 'test_dynamic_resolution_ll')[0]
+
+        self.assertIsInstance(test_area, DynamicAreaDefinition)
+        self.assertTrue(hasattr(test_area, 'resolution'))
+        self.assertEqual(test_area.resolution, (1.0, 1.0))
+
     def test_multiple_file_content(self):
         from pyresample import parse_area_file
         area_list = ["""ease_sh:
@@ -291,8 +312,23 @@ class TestMisc(unittest.TestCase):
         np.testing.assert_almost_equal(a, 6378137.)
         np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
 
+    def test_convert_proj_floats(self):
+        from collections import OrderedDict
+        import pyresample.utils as utils
+
+        pairs = [('proj', 'lcc'), ('ellps', 'WGS84'), ('lon_0', '-95'), ('no_defs', True)]
+        expected = OrderedDict([('proj', 'lcc'), ('ellps', 'WGS84'), ('lon_0', -95.0), ('no_defs', True)])
+        self.assertDictEqual(utils._proj4.convert_proj_floats(pairs), expected)
+
+        # EPSG
+        pairs = [('init', 'EPSG:4326'), ('EPSG', 4326)]
+        for pair in pairs:
+            expected = OrderedDict([pair])
+            self.assertDictEqual(utils._proj4.convert_proj_floats([pair]), expected)
+
     def test_proj4_str_dict_conversion(self):
         from pyresample import utils
+
         proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +no_defs"
         proj_dict = utils._proj4.proj4_str_to_dict(proj_str)
         proj_str2 = utils._proj4.proj4_dict_to_str(proj_dict)
@@ -301,6 +337,18 @@ class TestMisc(unittest.TestCase):
         self.assertIsInstance(proj_dict['lon_0'], float)
         self.assertIsInstance(proj_dict2['lon_0'], float)
 
+        # EPSG
+        expected = {'+init=EPSG:4326': {'init': 'EPSG:4326'},
+                    'EPSG:4326': {'EPSG': 4326}}
+
+        for proj_str, proj_dict_exp in expected.items():
+            proj_dict = utils._proj4.proj4_str_to_dict(proj_str)
+            self.assertEqual(proj_dict, proj_dict_exp)
+            self.assertEqual(utils._proj4.proj4_dict_to_str(proj_dict), proj_str)  # round-trip
+
+        # Invalid EPSG code (pyproj-2 syntax only)
+        self.assertRaises(ValueError, utils._proj4.proj4_str_to_dict, 'EPSG:XXXX')
+
     def test_def2yaml_converter(self):
         from pyresample import parse_area_file, convert_def_to_yaml
         import tempfile

pyresample/utils/__init__.py

@@ -18,6 +18,7 @@
 """Miscellaneous utility functions for pyresample."""
 from collections import Mapping
 import numpy as np
+import pyproj
 import warnings
 
 from ._proj4 import (proj4_dict_to_str, proj4_str_to_dict, convert_proj_floats, proj4_radius_parameters)  # noqa
@@ -234,3 +235,8 @@ def recursive_dict_update(d, u):
         else:
             d[k] = u[k]
     return d
+
+
+def is_pyproj2():
+    """Determine whether the current pyproj version is >= 2.0"""
+    return pyproj.__version__ >= '2'

pyresample/utils/_proj4.py

@@ -16,6 +16,7 @@
 # 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/>.
 from collections import OrderedDict
+import six
 
 
 def convert_proj_floats(proj_pairs):
@@ -25,6 +26,9 @@ def convert_proj_floats(proj_pairs):
         if len(x) == 1 or x[1] is True:
             proj_dict[x[0]] = True
             continue
+        if x[0] == 'EPSG':
+            proj_dict[x[0]] = x[1]
+            continue
 
         try:
             proj_dict[x[0]] = float(x[1])
@@ -39,6 +43,13 @@ def proj4_str_to_dict(proj4_str):
 
     Note: Key only parameters will be assigned a value of `True`.
     """
+    if proj4_str.startswith('EPSG:'):
+        try:
+            code = int(proj4_str.split(':', 1)[1])
+        except ValueError as err:
+            six.raise_from(ValueError("Invalid EPSG code '{}': {}".format(proj4_str, err)),
+                           None)  # Suppresses original exception context in python 3
+        return OrderedDict(EPSG=code)
     pairs = (x.split('=', 1) for x in proj4_str.replace('+', '').split(" "))
     return convert_proj_floats(pairs)
 
@@ -50,6 +61,11 @@ def proj4_dict_to_str(proj4_dict, sort=False):
         items = sorted(items)
     params = []
     for key, val in items:
+        if key == 'EPSG':
+            # If EPSG code is present, ignore other parameters
+            params = ['EPSG:{}'.format(val)]
+            break
+
         key = str(key) if key.startswith('+') else '+' + str(key)
         if key in ['+no_defs', '+no_off', '+no_rot']:
             param = key

pyresample/version.py

@@ -23,9 +23,9 @@ def get_keywords():
     # setup.py/versioneer.py will grep for the variable names, so they must
     # each be defined on a line of their own. _version.py will just call
     # get_keywords().
-    git_refnames = " (tag: v1.12.0)"
-    git_full = "5c4480e87898aaed2950e775324fec0bd408853c"
-    git_date = "2019-04-06 20:47:48 -0500"
+    git_refnames = " (tag: v1.12.3)"
+    git_full = "6b4367ce09ca6040c6107d2fafc31524eb45d1f8"
+    git_date = "2019-05-17 10:58:40 -0500"
     keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
     return keywords