CHANGELOG.md

+## Version 1.13.1 (2019/09/26)
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 218](https://github.com/pytroll/pyresample/pull/218) - Fix proj_str returning invalid PROJ strings when towgs84 was included
+* [PR 217](https://github.com/pytroll/pyresample/pull/217) - Fix get_geostationary_angle_extent assuming a/b definitions
+* [PR 216](https://github.com/pytroll/pyresample/pull/216) - Fix proj4 radius parameters for spherical cases
+
+In this release 3 pull requests were closed.
+
 ## Version 1.13.0 (2019/09/13)
 
 ### Issues Closed

debian/changelog

+pyresample (1.13.1-1) unstable; urgency=medium
+
+  * New upstream release.
+  * debian/patches:
+    - refresh all patches
+  * debian/control:
+    - explicitly specify Rules-Requires-Root: no
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Mon, 30 Sep 2019 08:01:23 +0200
+
 pyresample (1.13.0-1) unstable; urgency=medium
 
   * New upstream release.

debian/control

@@ -3,6 +3,7 @@ Maintainer: Debian GIS Project <pkg-grass-devel@lists.alioth.debian.org>
 Uploaders: Antonio Valentino <antonio.valentino@tiscali.it>
 Section: python
 Priority: optional
+Rules-Requires-Root: no
 Build-Depends: cython3,
                debhelper-compat (= 12),
                dh-python,

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

@@ -21,7 +21,7 @@ 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 ba3341e..7a9630b 100644
+index 387546f..2c6375b 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
 @@ -560,7 +560,7 @@ class Test(unittest.TestCase):

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

@@ -9,7 +9,7 @@ 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 7a9630b..0cd6ed2 100644
+index 2c6375b..c3f58dc 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
 @@ -9,7 +9,8 @@ import numpy as np
@@ -42,7 +42,7 @@ index 7a9630b..0cd6ed2 100644
      def test_get_proj_coords_dask(self):
          """Test get_proj_coords usage with dask arrays."""
          from pyresample import utils
-@@ -1424,6 +1431,8 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1456,6 +1463,8 @@ class TestSwathDefinition(unittest.TestCase):
          assert_np_dict_allclose(res.proj_dict, proj_dict)
          self.assertEqual(res.shape, (6, 3))
  
@@ -51,7 +51,7 @@ index 7a9630b..0cd6ed2 100644
      def test_aggregation(self):
          """Test aggregation on SwathDefinitions."""
          if (sys.version_info < (3, 0)):
-@@ -1441,6 +1450,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1473,6 +1482,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)
  

debian/patches/0003-Make-xarray-optional-for-testing.patch

@@ -9,7 +9,7 @@ Subject: Make xarray optional for testing
  3 files changed, 17 insertions(+), 1 deletion(-)
 
 diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
-index 0cd6ed2..de5fd7e 100644
+index c3f58dc..c63fa9c 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
 @@ -24,6 +24,11 @@ if sys.version_info < (2, 7):
@@ -24,7 +24,7 @@ index 0cd6ed2..de5fd7e 100644
  try:
      import dask
  except ImportError:
-@@ -1330,6 +1335,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1362,6 +1367,7 @@ class TestSwathDefinition(unittest.TestCase):
          self.assertFalse(
              swath_def == swath_def2, 'swath_defs are not expected to be equal')
  
@@ -32,7 +32,7 @@ index 0cd6ed2..de5fd7e 100644
      def test_compute_omerc_params(self):
          """Test omerc parameters computation."""
          lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
-@@ -1401,6 +1407,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1433,6 +1439,7 @@ class TestSwathDefinition(unittest.TestCase):
          np.testing.assert_allclose(lats, [80., 80., 80., 80., 80., 80., 80.,
                                            80., 80., 80., 80., 80.])
  
@@ -40,7 +40,7 @@ index 0cd6ed2..de5fd7e 100644
      def test_compute_optimal_bb(self):
          """Test computing the bb area."""
          from pyresample.utils import is_pyproj2
-@@ -1431,6 +1438,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1463,6 +1470,7 @@ class TestSwathDefinition(unittest.TestCase):
          assert_np_dict_allclose(res.proj_dict, proj_dict)
          self.assertEqual(res.shape, (6, 3))
  
@@ -48,7 +48,7 @@ index 0cd6ed2..de5fd7e 100644
      @unittest.skipIf(not hasattr(DataArray, 'coarsen'), 'DataArray.coarsen not available')
      @unittest.skipIf(not dask, 'dask not available')
      def test_aggregation(self):
-@@ -1450,6 +1458,7 @@ class TestSwathDefinition(unittest.TestCase):
+@@ -1482,6 +1490,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)
  
@@ -56,7 +56,7 @@ index 0cd6ed2..de5fd7e 100644
      @unittest.skipIf(not dask, 'dask not available')
      def test_striding(self):
          """Test striding."""
-@@ -1776,6 +1785,7 @@ class TestDynamicAreaDefinition(unittest.TestCase):
+@@ -1808,6 +1817,7 @@ class TestDynamicAreaDefinition(unittest.TestCase):
          self.assertEqual(result.width, 395)
          self.assertEqual(result.height, 539)
  

pyresample/geometry.py

@@ -36,7 +36,8 @@ from pyproj import Geod, transform
 from pyresample import CHUNK_SIZE
 from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP
 from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary
-from pyresample.utils import proj4_str_to_dict, proj4_dict_to_str, convert_proj_floats
+from pyresample.utils import (proj4_str_to_dict, proj4_dict_to_str,
+                              convert_proj_floats, proj4_radius_parameters)
 from pyresample.area_config import create_area_def
 
 try:
@@ -1346,7 +1347,17 @@ class AreaDefinition(BaseDefinition):
     @property
     def proj_str(self):
         """Return PROJ projection string."""
-        return proj4_dict_to_str(self.proj_dict, sort=True)
+        proj_dict = self.proj_dict.copy()
+        if 'towgs84' in proj_dict and isinstance(proj_dict['towgs84'], list):
+            # pyproj 2+ creates a list in the dictionary
+            # but the string should be comma-separated
+            if all(x == 0 for x in proj_dict['towgs84']):
+                # all 0s in towgs84 are technically equal to not having them
+                # specified, but PROJ considers them different
+                proj_dict.pop('towgs84')
+            else:
+                proj_dict['towgs84'] = ','.join(str(x) for x in proj_dict['towgs84'])
+        return proj4_dict_to_str(proj_dict, sort=True)
 
     def __str__(self):
         """Return string representation of the AreaDefinition."""
@@ -1904,8 +1915,9 @@ class AreaDefinition(BaseDefinition):
 def get_geostationary_angle_extent(geos_area):
     """Get the max earth (vs space) viewing angles in x and y."""
     # get some projection parameters
-    req = geos_area.proj_dict['a'] / 1000.0
-    rp = geos_area.proj_dict['b'] / 1000.0
+    a, b = proj4_radius_parameters(geos_area.proj_dict)
+    req = a / 1000.0
+    rp = b / 1000.0
     h = geos_area.proj_dict['h'] / 1000.0 + req
 
     # compute some constants

pyresample/test/test_geometry.py

@@ -1133,6 +1133,38 @@ class Test(unittest.TestCase):
                 area_extent=[-40000., -40000., 40000., 40000.])
             self.assertEqual(area.proj_str, expected_proj)
 
+        if utils.is_pyproj2():
+            # CRS with towgs84 in it
+            # we remove towgs84 if they are all 0s
+            projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000,
+                          'ellps': 'GRS80', 'towgs84': '0,0,0,0,0,0,0', 'units': 'm', 'no_defs': True}
+            area = geometry.AreaDefinition(
+                area_id='test_towgs84',
+                description='',
+                proj_id='',
+                projection=projection,
+                width=123, height=123,
+                area_extent=[-40000., -40000., 40000., 40000.])
+            self.assertEqual(area.proj_str,
+                             '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea '
+                             # '+towgs84=0.0,0.0,0.0,0.0,0.0,0.0,0.0 '
+                             '+type=crs +units=m '
+                             '+x_0=4321000 +y_0=3210000')
+            projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000,
+                          'ellps': 'GRS80', 'towgs84': '0,5,0,0,0,0,0', 'units': 'm', 'no_defs': True}
+            area = geometry.AreaDefinition(
+                area_id='test_towgs84',
+                description='',
+                proj_id='',
+                projection=projection,
+                width=123, height=123,
+                area_extent=[-40000., -40000., 40000., 40000.])
+            self.assertEqual(area.proj_str,
+                             '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea '
+                             '+towgs84=0.0,5.0,0.0,0.0,0.0,0.0,0.0 '
+                             '+type=crs +units=m '
+                             '+x_0=4321000 +y_0=3210000')
+
     def test_striding(self):
         """Test striding AreaDefinitions."""
         from pyresample import utils
@@ -1867,6 +1899,13 @@ class TestCrop(unittest.TestCase):
 
         expected = (0.15185342867090912, 0.15133555510297725)
 
+        np.testing.assert_allclose(expected,
+                                   geometry.get_geostationary_angle_extent(geos_area))
+
+        geos_area.proj_dict = {'ellps': 'GRS80',
+                               'h': 35785831.00}
+        expected = (0.15185277703584374, 0.15133971368991794)
+
         np.testing.assert_allclose(expected,
                                    geometry.get_geostationary_angle_extent(geos_area))
 

pyresample/test/test_utils.py

@@ -330,6 +330,7 @@ class TestMisc(unittest.TestCase):
                            1000, 1000, (-1000, -1000, 1000, 1000))
 
     def test_proj4_radius_parameters_provided(self):
+        """Test proj4_radius_parameters with a/b."""
         from pyresample import utils
         a, b = utils._proj4.proj4_radius_parameters(
             '+proj=stere +a=6378273 +b=6356889.44891',
@@ -338,6 +339,7 @@ class TestMisc(unittest.TestCase):
         np.testing.assert_almost_equal(b, 6356889.44891)
 
     def test_proj4_radius_parameters_ellps(self):
+        """Test proj4_radius_parameters with ellps."""
         from pyresample import utils
         a, b = utils._proj4.proj4_radius_parameters(
             '+proj=stere +ellps=WGS84',
@@ -346,6 +348,7 @@ class TestMisc(unittest.TestCase):
         np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
 
     def test_proj4_radius_parameters_default(self):
+        """Test proj4_radius_parameters with default parameters."""
         from pyresample import utils
         a, b = utils._proj4.proj4_radius_parameters(
             '+proj=lcc',
@@ -354,6 +357,15 @@ class TestMisc(unittest.TestCase):
         np.testing.assert_almost_equal(a, 6378137.)
         np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
 
+    def test_proj4_radius_parameters_spherical(self):
+        """Test proj4_radius_parameters in case of a spherical earth."""
+        from pyresample import utils
+        a, b = utils._proj4.proj4_radius_parameters(
+            '+proj=stere +R=6378273',
+        )
+        np.testing.assert_almost_equal(a, 6378273.)
+        np.testing.assert_almost_equal(b, 6378273.)
+
     def test_convert_proj_floats(self):
         from collections import OrderedDict
         import pyresample.utils as utils

pyresample/utils/_proj4.py

@@ -113,6 +113,9 @@ def proj4_radius_parameters(proj4_dict):
             new_info['b'] = float(new_info['a']) * (1 - float(new_info['f']))
         elif 'b' in new_info and 'f' in new_info:
             new_info['a'] = float(new_info['b']) / (1 - float(new_info['f']))
+        elif 'R' in new_info:
+            new_info['a'] = new_info['R']
+            new_info['b'] = new_info['R']
         else:
             geod = Geod(**{'ellps': 'WGS84'})
             new_info['a'] = geod.a

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 = " (HEAD -> master, tag: v1.13.0)"
-    git_full = "790a0bae85cb243c17f0150011e30c834f244e04"
-    git_date = "2019-09-13 08:32:20 +0200"
+    git_refnames = " (HEAD -> master, tag: v1.13.1)"
+    git_full = "d9ff2c9012a4dbfd0f39172fde2e21bb34c04e4a"
+    git_date = "2019-09-26 20:22:21 +0200"
     keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
     return keywords