Merge pull request #1018 from CLIMADA-project/hotfix_bounding_box_from_countries

hotfix: bug in bounding_box_from_countries

Author: ValentinGebhart

climada/util/coordinates.py

@@ -46,7 +46,7 @@
 import shapely.vectorized
 import shapely.wkt
 from cartopy.io import shapereader
-from shapely.geometry import MultiPolygon, Point, Polygon, box
+from shapely.geometry import MultiPolygon, Point, box
 from sklearn.neighbors import BallTree
 
 import climada.util.hdf5_handler as u_hdf5
@@ -1723,19 +1723,20 @@ def bounding_box_from_countries(country_names, buffer=1.0):
     """
 
     country_geometry = get_country_geometries(country_names).geometry
-    longitudes, latitudes = [], []
-    for multipolygon in country_geometry:
-        if isinstance(multipolygon, Polygon):  # if entry is polygon
-            for coord in polygon.exterior.coords:  # Extract exterior coordinates
-                longitudes.append(coord[0])
-                latitudes.append(coord[1])
-        else:  # if entry is multipolygon
-            for polygon in multipolygon.geoms:
-                for coord in polygon.exterior.coords:  # Extract exterior coordinates
-                    longitudes.append(coord[0])
-                    latitudes.append(coord[1])
+    generator_country_geometry = (
+        poly if isinstance(poly, MultiPolygon) else MultiPolygon([poly])
+        for poly in country_geometry
+    )
+    lon, lat = np.concatenate(
+        [
+            np.array(pg.exterior.coords).T
+            for poly in generator_country_geometry
+            for pg in poly.geoms
+        ],
+        axis=1,
+    )
 
-    return latlon_bounds(np.array(latitudes), np.array(longitudes), buffer=buffer)
+    return latlon_bounds(lat, lon, buffer=buffer)
 
 
 def bounding_box_from_cardinal_bounds(*, northern, eastern, western, southern):
@@ -1756,12 +1757,13 @@ def bounding_box_from_cardinal_bounds(*, northern, eastern, western, southern):
     Returns
     -------
     tuple
-        The resulting normalized bounding box (min_lon, min_lat, max_lon, max_lat) with -180 <= min_lon < max_lon < 540
+        The resulting normalized bounding box (min_lon, min_lat, max_lon, max_lat)
+        with -180 <= min_lon < max_lon < 540
 
     """
 
     # latitude bounds check
-    if not ((90 >= northern > southern >= -90)):
+    if not 90 >= northern > southern >= -90:
         raise ValueError(
             "Given northern bound is below given southern bound or out of bounds"
         )

climada/util/test/test_coordinates.py

@@ -2305,6 +2305,7 @@ def test_bounding_box_global(self):
 
     def test_bounding_box_from_countries(self):
         """Test for a list of ISO country codes."""
+        # Italy is a multipolygon geometry
         result = u_coord.bounding_box_from_countries(
             ["ITA"], buffer=1.0
         )  # Testing with Italy (ITA)
@@ -2314,7 +2315,21 @@ def test_bounding_box_from_countries(self):
             34.48924388200004,
             19.517425977000073,
             48.08521494500006,
-        ]  # Italy's bounding box
+        ]  # Italy's bounding box with 1 degree buffer
+        np.testing.assert_array_almost_equal(result, expected)
+
+        # Switzerland is a polygon geometry
+        result = u_coord.bounding_box_from_countries(
+            ["CHE"], buffer=0.0
+        )  # Testing with Switzerland (CHE)
+        # Real expected bounds for Switzerland (calculated or manually known)
+        expected = [
+            5.954809204000128,
+            45.82071848599999,
+            10.466626831000013,
+            47.801166077000076,
+        ]  # CHE's bounding box with 0 degree buffer
+        np.testing.assert_array_almost_equal(result, expected)
 
         # invalid input
         with self.assertRaises(ValueError):