Fix geostationary bbox having inf values

mraspaud · Jan 26, 2023 · 9a9a16f · 9a9a16f
1 parent 8393134
commit 9a9a16f
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 3 deletions.
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
@@ -1556,7 +1556,7 @@ def _get_geo_boundary_sides(self, frequency=None):
         # Ensure an even number of vertices for side creation
         if (frequency % 2) != 0:
             frequency = frequency + 1
-        lons, lats = get_geostationary_bounding_box(self, nb_points=frequency)
+        lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=frequency)
         # Retrieve dummy sides for GEO (side1 and side3 always of length 2)
         side02_step = int(frequency / 2) - 1
         lon_sides = [lons[slice(0, side02_step + 1)],
@@ -2818,6 +2818,8 @@ def get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
     """
     x, y = get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
     lons, lats = Proj(geos_area.crs)(x, y, inverse=True)
+    lons[np.isinf(lons)] = np.nan
+    lats[np.isinf(lats)] = np.nan
     return lons, lats
 
 

diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
@@ -2606,7 +2606,7 @@ def test_get_geostationary_bbox(self):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
         # This musk be equal to lon.
         elon = np.array([-79.23372832, -78.19662326, -75.42516215, -70.22636028,
                          -56.89851775, 0., 56.89851775, 70.22636028,
@@ -2634,9 +2634,22 @@ def test_get_geostationary_bbox(self):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
         np.testing.assert_allclose(lon, elon + lon_0)
 
+    def test_get_geostationary_bbox_does_not_contain_inf(self):
+        """Ensure the geostationary bbox does not contain np.inf."""
+        projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos',
+                      'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
+        area_extent = (5567248.0742, 5570248.4773, -5570248.4773, 1393687.2705)
+        width = 3712
+        height = 1392
+        geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent)
+
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
+        assert not any(np.isinf(lon))
+        assert not any(np.isinf(lat))
+
     def test_get_geostationary_angle_extent(self):
         """Get max geostationary angles."""
         geos_area = MagicMock()