Add tropical cyclones basin bounds

CHANGELOG.md

@@ -72,6 +72,7 @@ Removed:
 
 ### Added
 
+- `climada.hazard.tc_tracks.BASINS_BOUNDS` dictionary [#1031](https://github.com/CLIMADA-project/climada_python/pull/1031)
 - `climada.hazard.tc_tracks.TCTracks.subset_years` function [#1023](https://github.com/CLIMADA-project/climada_python/pull/1023)
 - `climada.hazard.tc_tracks.TCTracks.from_FAST` function, add Australia basin (AU) [#993](https://github.com/CLIMADA-project/climada_python/pull/993)
 - Add `osm-flex` package to CLIMADA core [#981](https://github.com/CLIMADA-project/climada_python/pull/981)

climada/hazard/tc_tracks.py

@@ -28,7 +28,8 @@
 import re
 import shutil
 import warnings
-from operator import itemgetter
+from collections import defaultdict
+from enum import Enum
 from pathlib import Path
 from typing import List, Optional
 
@@ -50,7 +51,8 @@
 from matplotlib.collections import LineCollection
 from matplotlib.colors import BoundaryNorm, ListedColormap
 from matplotlib.lines import Line2D
-from shapely.geometry import LineString, MultiLineString, Point
+from shapely.geometry import LineString, MultiLineString, Point, Polygon
+from shapely.ops import unary_union
 from sklearn.metrics import DistanceMetric
 
 import climada.hazard.tc_tracks_synth
@@ -193,7 +195,98 @@
 dataset using STORM. Scientific Data 7(1): 40."""
 
 
+class Basin(Enum):
+    """
+    Store tropical cyclones basin geographical extent.
+    The boundaries of the basin are represented as a polygon (using the `shapely` Polygon object)
+    and follows the definition of the STORM dataset. Important note: tropical cyclone boundaries
+    may vary bewteen datasets. The following boundaries follows the STORM definition:
+    https://www.nature.com/articles/s41597-020-0381-2
+
+    Attributes:
+    ----------
+        *name : str
+            The name of the tropical cyclone basin (e.g., "NA" for North Atlantic).
+        *polygon : Polygon
+            A shapely Polygon object that represents the geographical boundary of the basin.
+
+    """
+
+    NA = Polygon(
+        [
+            (-100, 19),
+            (-94.21951983987083, 17.039584804350312),
+            (-88.75211790888072, 14.837521327451947),
+            (-84.96610530622198, 12.214318798718033),
+            (-84.89823142225451, 12.181148019885352),
+            (-82.59052306410497, 8.777858931465238),
+            (-81.09730008320902, 8.358383265470449),
+            (-79.50226644452471, 9.196860922133856),
+            (-78.58597052442947, 9.213610839871123),
+            (-77.02487377167459, 7.299350879751048),
+            (-77.02487377167459, 5),
+            (0.0, 5.0),
+            (0.0, 60.0),
+            (-100.0, 60.0),
+            (-100, 19),
+        ]
+    )
+
+    EP = Polygon(
+        [
+            (-180.0, 5.0),
+            (-77.02487377167459, 5),
+            (-77.02487377167459, 7.299350879751048),
+            (-78.58597052442947, 9.213610839871123),
+            (-79.50226644452471, 9.196860922133856),
+            (-81.09730008320902, 8.358383265470449),
+            (-82.59052306410497, 8.777858931465238),
+            (-84.89823142225451, 12.181148019885352),
+            (-84.96610530622198, 12.214318798718033),
+            (-88.75211790888072, 14.837521327451947),
+            (-94.21951983987083, 17.039584804350312),
+            (-100, 19),
+            (-100.0, 60.0),
+            (-180.0, 60.0),
+            (-180.0, 5.0),
+        ]
+    )
+
+    WP = Polygon(
+        [(100.0, 5.0), (180.0, 5.0), (180.0, 60.0), (100.0, 60.0), (100.0, 5.0)]
+    )
+
+    NI = Polygon([(30.0, 5.0), (100.0, 5.0), (100.0, 60.0), (30.0, 60.0), (30.0, 5.0)])
+
+    SI = Polygon(
+        [(10.0, -60.0), (135.0, -60.0), (135.0, -5.0), (10.0, -5.0), (10.0, -60.0)]
+    )
+
+    SP = unary_union(
+        [
+            Polygon(  # west side of antimeridian
+                [
+                    (135.0, -60.0),
+                    (180.0, -60.0),
+                    (180.0, -5.0),
+                    (135.0, -5.0),
+                    (135.0, -60.0),
+                ]
+            ),
+            Polygon(  # east side
+                [
+                    (-180.0, -60.0),
+                    (-120.0, -60.0),
+                    (-120.0, -5.0),
+                    (-180.0, -5.0),
+                    (-180.0, -60.0),
+                ]
+            ),
+        ]
+    )
+
+
 class TCTracks:
     """Contains tropical cyclone tracks.
 
     Attributes
@@ -322,7 +415,69 @@
 
         return out
 
+    def subset_by_basin(self):
+        """Subset all tropical cyclones tracks by basin.
+
+        This function iterates through the tropical cyclones in the dataset and assigns each cyclone
+        to a basin based on its geographical location. It checks whether the cyclone's position
+        (latitude and longitude) lies within the boundaries of any of the predefined basins and
+        then groups the cyclones into separate categories for each basin. The resulting dictionary
+        maps each basin's name to a list of tropical cyclones that fall within it.
+
+        Parameters
+        ----------
+        self : TCTtracks object
+            The object instance containing the tropical cyclone data (`self.data`) to be processed.
+
+        Returns
+        -------
+        dict_tc_basins : dict
+            A dictionary where the keys are basin names (e.g., "NA", "EP", "WP", etc.) and the
+            values are instances of the `TCTracks` class containing the tropical cyclones that
+            belong to each basin.
+
+        Example:
+        --------
+        >>> tc = TCTracks.from_ibtracks("")
+        >>> tc_basins = tc.split_by_basin()
+        >>> tc_basins["NA"] # to access tracks in the North Atlantic
+
+        """
+
+        # Initialize a defaultdict to store lists for each basin
+        basins_dict = defaultdict(list)
+        tracks_outside_basin: list = []
+        # Iterate over each tropical cyclone
+        for track in self.data:
+            lat, lon = track.lat.values[0], track.lon.values[0]
+            origin_point = Point(lon, lat)
+            point_in_basin = False
+
+            # Find the basin that contains the point
+            for basin in Basin:
+                if basin.value.contains(origin_point):
+                    basins_dict[basin.name].append(track)
+                    point_in_basin = True
+                    break
+
+            if not point_in_basin:
+                tracks_outside_basin.append(track.id_no)
+
+        if tracks_outside_basin:
+            warnings.warn(
+                f"A total of {len(tracks_outside_basin)} tracks did not originate in any of the \n"
+                f"defined basins. IDs of the tracks outside the basins: {tracks_outside_basin}",
+                UserWarning,
+            )
+
+        # Create a dictionary with TCTracks for each basin
+        dict_tc_basins = {
+            basin_name: TCTracks(tc_list) for basin_name, tc_list in basins_dict.items()
+        }
+
+        return dict_tc_basins
+
     def subset_year(
         self,
         start_date: tuple = (False, False, False),
         end_date: tuple = (False, False, False),

climada/hazard/test/test_tc_tracks.py

@@ -874,6 +874,28 @@ def test_subset_years(self):
         ):
             tc_test.subset_year((2100, False, False), (2150, False, False))
 
+    def test_subset_basin(self):
+        """test the correct splitting of a single tc object into different tc objets by basin"""
+
+        tc_test = tc.TCTracks.from_simulations_emanuel(TEST_TRACK_EMANUEL)
+        tc_test.data[-1].lat[0] = 0  # modify lat of track to exclude it from a basin
+
+        with self.assertWarnsRegex(
+            UserWarning,
+            "A total of 1 tracks did not originate in any of the \n"
+            "defined basins. IDs of the tracks outside the basins: \[4\]",
+        ):
+            dict_basins = tc_test.subset_by_basin()
+
+        self.assertEqual(dict_basins["EP"].data[0].lat[0].item(), 12.553)
+        self.assertEqual(dict_basins["EP"].data[0].lon[0].item(), -109.445)
+        self.assertEqual(dict_basins["SI"].data[0].lat[0].item(), -8.699)
+        self.assertEqual(dict_basins["SI"].data[0].lon[0].item(), 52.761)
+        self.assertEqual(dict_basins["WP"].data[0].lat[0].item(), 8.502)
+        self.assertEqual(dict_basins["WP"].data[0].lon[0].item(), 164.909)
+        self.assertEqual(dict_basins["WP"].data[1].lat[0].item(), 16.234)
+        self.assertEqual(dict_basins["WP"].data[1].lon[0].item(), 116.424)
+
     def test_get_extent(self):
         """Test extent/bounds attributes."""
         storms = ["1988169N14259", "2002073S16161", "2002143S07157"]