diff --git a/CHANGELOG.md b/CHANGELOG.md index ee2d3bdf60..7dc6793de0 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/climada/util/coordinates.py b/climada/util/coordinates.py index cec74b512c..351263e629 100644 --- a/climada/util/coordinates.py +++ b/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 diff --git a/climada/util/test/test_coordinates.py b/climada/util/test/test_coordinates.py index 50d5a8073e..86f95f764d 100644 --- a/climada/util/test/test_coordinates.py +++ b/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)