Merge pull request #221 from mraspaud/fix-#220

Make optimal bb computation work with numpy arrays
pytroll · Oct 8, 2019 · 245f75c · 245f75c
2 parents d9ff2c9 + 0687994
commit 245f75c
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 20 deletions.
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
@@ -698,18 +698,19 @@ def compute_bb_proj_params(self, proj_dict):
     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)
+
+        def notnull(arr):
+            try:
+                return arr.where(arr.notnull(), drop=True)
+            except AttributeError:
+                return arr[np.isfinite(arr)]
+
+        leftlons = notnull(self.lons[:, 0])
+        rightlons = notnull(self.lons[:, -1])
+        middlelons = notnull(self.lons[:, int(self.lons.shape[1] / 2)])
+        leftlats = notnull(self.lats[:, 0])
+        rightlats = notnull(self.lats[:, -1])
+        middlelats = notnull(self.lats[:, int(self.lats.shape[1] / 2)])
 
         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])

diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
@@ -1430,14 +1430,14 @@ def test_compute_optimal_bb(self):
         """Test computing the bb area."""
         from pyresample.utils import is_pyproj2
         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)
+        nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+                           [80.84000396728516, 60.74200439453125, 34.08500289916992],
+                           [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+        lats = xr.DataArray(nplats)
+        nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+                           [79.11000061035156, 7.284000396728516, -5.107000350952148],
+                           [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+        lons = xr.DataArray(nplons)
 
         area = geometry.SwathDefinition(lons, lats)
 
@@ -1456,6 +1456,23 @@ def test_compute_optimal_bb(self):
         assert_np_dict_allclose(res.proj_dict, proj_dict)
         self.assertEqual(res.shape, (6, 3))
 
+        area = geometry.SwathDefinition(nplons, nplats)
+
+        res = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'})
+
+        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}
+        if is_pyproj2():
+            # pyproj2 adds some extra defaults
+            proj_dict.update({'x_0': 0, 'y_0': 0, 'units': 'm',
+                              'k': 1, 'gamma': 0,
+                              'no_defs': None, 'type': 'crs'})
+        assert_np_dict_allclose(res.proj_dict, proj_dict)
+        self.assertEqual(res.shape, (6, 3))
+
     def test_aggregation(self):
         """Test aggregation on SwathDefinitions."""
         if (sys.version_info < (3, 0)):