Merge pull request #475 from ghiggi/bugfix-duplicate-geo-bbox-lonlats

closes #474
Closes undefined

Author: djhoese

pyresample/geometry.py

@@ -2691,24 +2691,30 @@ def get_geostationary_angle_extent(geos_area):
 def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
     """Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`.
 
-    Args:
-      nb_points: Number of points on the polygon
+    Notes:
+    - The first and last element of the output vectors are equal.
+    - If nb_points is even, it will return x and y vectors of length nb_points + 1.
+
+    Parameters
+    ----------
+    nb_points : Number of points on the polygon.
+
     """
     xmax, ymax = get_geostationary_angle_extent(geos_area)
     h = get_geostationary_height(geos_area.crs)
 
     # generate points around the north hemisphere in satellite projection
     # make it a bit smaller so that we stay inside the valid area
-    x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0))) * (xmax - 0.0001)
-    y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0))) * (ymax - 0.0001)
+    x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (xmax - 0.0001)
+    y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (ymax - 0.0001)
 
     ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
 
     x *= h
     y *= h
-
-    x = np.clip(np.concatenate([x, x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
-    y = np.clip(np.concatenate([y, -y]), min(ll_y, ur_y), max(ll_y, ur_y))
+    # We remove one element with [:-1] to avoid duplicate values at the equator
+    x = np.clip(np.concatenate([x[:-1], x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
+    y = np.clip(np.concatenate([y[:-1], -y]), min(ll_y, ur_y), max(ll_y, ur_y))
 
     return x, y
 

pyresample/test/test_geometry.py

@@ -2603,19 +2603,21 @@ def test_get_geostationary_bbox(self):
 
         lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
         # This musk be equal to lon.
-        elon = np.array([-79.23372832, -77.9694809, -74.55229623, -67.32816598,
-                         -41.45591465, 41.45591465, 67.32816598, 74.55229623,
-                         77.9694809, 79.23372832, 79.23372832, 77.9694809,
-                         74.55229623, 67.32816598, 41.45591465, -41.45591465,
-                         -67.32816598, -74.55229623, -77.9694809, -79.23372832])
-        elat = np.array([6.94302533e-15, 1.97333299e+01, 3.92114217e+01, 5.82244715e+01,
-                         7.52409201e+01, 7.52409201e+01, 5.82244715e+01, 3.92114217e+01,
-                         1.97333299e+01, -0.00000000e+00, -6.94302533e-15, -1.97333299e+01,
-                         -3.92114217e+01, -5.82244715e+01, -7.52409201e+01, -7.52409201e+01,
-                         -5.82244715e+01, -3.92114217e+01, -1.97333299e+01, 0.0])
-
-        np.testing.assert_allclose(lon, elon)
-        np.testing.assert_allclose(lat, elat)
+        elon = np.array([-79.23372832, -78.19662326, -75.42516215, -70.22636028,
+                         -56.89851775, 0., 56.89851775, 70.22636028,
+                         75.42516215, 78.19662326, 79.23372832, 78.19662326,
+                         75.42516215, 70.22636028, 56.89851775, 0.,
+                        -56.89851775, -70.22636028, -75.42516215, -78.19662326,
+                        -79.23372832])
+        elat = np.array([0., 17.76879577, 35.34328897, 52.5978607,
+                         69.00533142, 79.14811219, 69.00533142, 52.5978607,
+                         35.34328897, 17.76879577, -0., -17.76879577,
+                         -35.34328897, -52.5978607, -69.00533142, -79.14811219,
+                         -69.00533142, -52.5978607, -35.34328897, -17.76879577,
+                         0.])
+
+        np.testing.assert_allclose(lon, elon, atol=1e-07)
+        np.testing.assert_allclose(lat, elat, atol=1e-07)
 
         geos_area = MagicMock()
         lon_0 = 10
@@ -2825,7 +2827,7 @@ def test_get_area_slices(self):
         assert isinstance(slice_x.start, int)
         assert isinstance(slice_y.start, int)
         self.assertEqual(slice(46, 3667, None), slice_x)
-        self.assertEqual(slice(52, 3663, None), slice_y)
+        self.assertEqual(slice(56, 3659, None), slice_y)
 
         area_to_cover = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
                                                 {'a': 6378144.0,
@@ -2883,7 +2885,7 @@ def test_get_area_slices(self):
             assert isinstance(slice_x.start, int)
             assert isinstance(slice_y.start, int)
             self.assertEqual(slice_x, slice(46, 3667, None))
-            self.assertEqual(slice_y, slice(52, 3663, None))
+            self.assertEqual(slice_y, slice(56, 3659, None))
 
     def test_get_area_slices_nongeos(self):
         """Check area slicing for non-geos projections."""