Merge pull request #980 from CLIMADA-project/feature/bounds_NESW_func…

…tions

add bounding box functions to util.coordinates

CHANGELOG.md

@@ -12,6 +12,9 @@ Code freeze date: YYYY-MM-DD
 
 ### Added
 
+- `climada.util.coordinates.bounding_box_global` function [#980](https://github.com/CLIMADA-project/climada_python/pull/980)
+- `climada.util.coordinates.bounding_box_from_countries` function [#980](https://github.com/CLIMADA-project/climada_python/pull/980)
+- `climada.util.coordinates.bounding_box_from_cardinal_bounds` function [#980](https://github.com/CLIMADA-project/climada_python/pull/980)
 - `climada.engine.impact.Impact.local_return_period` method [#971](https://github.com/CLIMADA-project/climada_python/pull/971)
 - `doc.tutorial.climada_util_local_exceedance_values.ipynb` tutorial explaining `Hazard.local_exceedance_intensity`, `Hazard.local_return_period`, `Impact.local_exceedance_impact`, and `Impact.local_return_period` methods [#971](https://github.com/CLIMADA-project/climada_python/pull/971)
 - `Hazard.local_exceedance_intensity`, `Hazard.local_return_period` and `Impact.local_exceedance_impact`, that all use the `climada.util.interpolation` module [#918](https://github.com/CLIMADA-project/climada_python/pull/918)
@@ -28,6 +31,7 @@ Code freeze date: YYYY-MM-DD
 
 ### Changed
 
+- `climada.util.coordinates.get_country_geometries` function: Now throwing a ValueError if unregognized ISO country code is given (before, the invalid ISO code was ignored) [#980](https://github.com/CLIMADA-project/climada_python/pull/980)
 - Improved scaling factors implemented in `climada.hazard.trop_cyclone.apply_climate_scenario_knu` to model the impact of climate changes to tropical cyclones [#734](https://github.com/CLIMADA-project/climada_python/pull/734)
 - In `climada.util.plot.geo_im_from_array`, NaNs are plotted in gray while cells with no centroid are not plotted [#929](https://github.com/CLIMADA-project/climada_python/pull/929)
 - Renamed `climada.util.plot.subplots_from_gdf` to `climada.util.plot.plot_from_gdf` [#929](https://github.com/CLIMADA-project/climada_python/pull/929)

climada/util/coordinates.py

@@ -783,6 +783,12 @@ def get_country_geometries(
     if country_names:
         if isinstance(country_names, str):
             country_names = [country_names]
+
+        # raise error if a country name is not recognized
+        for country_name in country_names:
+            if not country_name in nat_earth[["ISO_A3", "WB_A3", "ADM0_A3"]].values:
+                raise ValueError(f"ISO code {country_name} not recognized.")
+
         country_mask = np.isin(
             nat_earth[["ISO_A3", "WB_A3", "ADM0_A3"]].values,
             country_names,
@@ -1687,6 +1693,89 @@ def _ensure_utf8(val):
     return admin1_info, admin1_shapes
 
 
+def bounding_box_global():
+    """
+    Return global bounds in EPSG 4326
+
+    Returns
+    -------
+    tuple:
+        The global bounding box as (min_lon, min_lat, max_lon, max_lat)
+    """
+    return (-180, -90, 180, 90)
+
+
+def bounding_box_from_countries(country_names, buffer=1.0):
+    """
+    Return bounding box in EPSG 4326 containing given countries.
+
+    Parameters
+    ----------
+    country_names : list or str
+        list with ISO 3166 alpha-3 codes of countries, e.g ['ZWE', 'GBR', 'VNM', 'UZB']
+    buffer : float, optional
+        Buffer to add to both sides of the bounding box. Default: 1.0.
+
+    Returns
+    -------
+    tuple
+        The bounding box containing all given coutries as (min_lon, min_lat, max_lon, max_lat)
+    """
+
+    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])
+
+    return latlon_bounds(np.array(latitudes), np.array(longitudes), buffer=buffer)
+
+
+def bounding_box_from_cardinal_bounds(*, northern, eastern, western, southern):
+    """
+    Return and normalize bounding box in EPSG 4326 from given cardinal bounds.
+
+    Parameters
+    ----------
+    northern : (int, float)
+        Northern boundary of bounding box
+    eastern : (int, float)
+        Eastern boundary of bounding box
+    western : (int, float)
+        Western boundary of bounding box
+    southern : (int, float)
+        Southern boundary of bounding box
+
+    Returns
+    -------
+    tuple
+        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)):
+        raise ValueError(
+            "Given northern bound is below given southern bound or out of bounds"
+        )
+
+    eastern = (eastern + 180) % 360 - 180
+    western = (western + 180) % 360 - 180
+
+    # Ensure eastern > western
+    if western > eastern:
+        eastern += 360
+
+    return (western, southern, eastern, northern)
+
+
 def get_admin1_geometries(countries):
     """
     return geometries, names and codes of admin 1 regions in given countries

climada/util/test/test_coordinates.py

@@ -2294,6 +2294,80 @@ def test_mask_raster_with_geometry(self):
         )
 
 
+class TestBoundsFromUserInput(unittest.TestCase):
+    """Unit tests for the bounds_from_user_input function."""
+
+    def test_bounding_box_global(self):
+        """Test for 'global' area selection."""
+        result = u_coord.bounding_box_global()
+        expected = (-180, -90, 180, 90)
+        np.testing.assert_almost_equal(result, expected)
+
+    def test_bounding_box_from_countries(self):
+        """Test for a list of ISO country codes."""
+        result = u_coord.bounding_box_from_countries(
+            ["ITA"], buffer=1.0
+        )  # Testing with Italy (ITA)
+        # Real expected bounds for Italy (calculated or manually known)
+        expected = [
+            5.6027283120000675,
+            34.48924388200004,
+            19.517425977000073,
+            48.08521494500006,
+        ]  # Italy's bounding box
+
+        # invalid input
+        with self.assertRaises(ValueError):
+            u_coord.bounding_box_from_countries(["invalid_ISO", "DEU"])
+
+    def test_bounding_box_from_cardinal_bounds(self):
+        """Test for conversion from cardinal bounds to bounds."""
+        np.testing.assert_array_almost_equal(
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern=90, southern=-20, eastern=30, western=20
+            ),
+            (20, -20, 30, 90),
+        )
+        np.testing.assert_array_almost_equal(
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern=90, southern=-20, eastern=20, western=30
+            ),
+            (30, -20, 380, 90),
+        )
+        np.testing.assert_array_almost_equal(
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern=90, southern=-20, eastern=170, western=-170
+            ),
+            (-170, -20, 170, 90),
+        )
+        np.testing.assert_array_almost_equal(
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern=90, southern=-20, eastern=-170, western=170
+            ),
+            (170, -20, 190, 90),
+        )
+        np.testing.assert_array_almost_equal(
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern=90, southern=-20, eastern=170, western=175
+            ),
+            (175, -20, 530, 90),
+        )
+
+        # some invalid cases
+        with self.assertRaises(TypeError):
+            u_coord.bounding_box_from_cardinal_bounds(
+                southern=-20, eastern=30, western=20
+            )
+        with self.assertRaises(TypeError):
+            u_coord.bounding_box_from_cardinal_bounds([90, -20, 30, 20])
+        with self.assertRaises(TypeError):
+            u_coord.bounding_box_from_cardinal_bounds(90, -20, 30, 20)
+        with self.assertRaises(TypeError):
+            u_coord.bounding_box_from_cardinal_bounds(
+                northern="90", southern=-20, eastern=30, western=20
+            )
+
+
 # Execute Tests
 if __name__ == "__main__":
     TESTS = unittest.TestLoader().loadTestsFromTestCase(TestFunc)
@@ -2302,4 +2376,5 @@ def test_mask_raster_with_geometry(self):
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRasterMeta))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRasterIO))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestDistance))
+    TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestBoundsFromUserInput))
     unittest.TextTestRunner(verbosity=2).run(TESTS)