first version of method from_netcdf_fast and set up the test

Author: NicolasColombi

climada/hazard/tc_tracks.py

@@ -162,6 +162,7 @@
     "SI": 1005,
     "WP": 1005,
     "SP": 1004,
+    "AU": 1004,
 }
 """Basin-specific default environmental pressure"""
 
@@ -1619,6 +1620,92 @@ def from_netcdf(cls, folder_name):
                 data.append(track)
         return cls(data)
 
+    @classmethod
+    def from_netcdf_fast(cls, folder_name):
+        """Create new TCTracks object from NetCDF files created with the FAST model
+        of Jonathan Lin.
+
+        GitHub Repository: https://github.com/linjonathan/tropical_cyclone_risk?
+                            tab=readme-ov-file
+        Publication: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2023MS003686
+
+        Parameters:
+        ----------
+        folder_name : str
+            Folder name from where to read files.
+
+        Returns:
+        -------
+        tracks : TCTracks
+            TCTracks obecjt with tracks data from the given directory of NetCDF files.
+        """
+
+        file_tr = get_file_names(folder_name)
+        LOGGER.info("Reading %s files.", len(file_tr))
+        data = []
+        for file in file_tr:
+            if Path(file).suffix != ".nc":
+                continue
+            with xr.open_dataset(file) as ds:
+                for i in ds.n_trk:
+                    # Select track
+                    track = ds.sel(n_trk=i.item())
+
+                    # Define coordinates
+                    lat = track.lat_trks.data
+                    lon = track.lon_trks.data
+                    time = track.time.data
+
+                    # Define variables
+                    time_step_vector = np.full(time.shape[0], track.time.data[1])
+                    max_sustained_wind = track.v_trks.data
+                    basin_vector = np.full(time.shape[0], track.tc_basins.data.item())
+                    central_pressure = np.nan  # work in progress: get them from model
+                    radius_max_wind = np.nan  # work in progress: get them from model
+                    env_pressure = BASIN_ENV_PRESSURE[track.tc_basins.data.item()]
+                    env_pressure_vect = np.full(time.shape[0], env_pressure)
+
+                    # Define hurricaine category and other attributes
+                    max_sustained_wind_kn = np.nanmax(
+                        max_sustained_wind * 1.943844
+                    )  # convert from m/s to knots
+                    category_test = np.full(
+                        len(SAFFIR_SIM_CAT), max_sustained_wind_kn
+                    ) < np.array(SAFFIR_SIM_CAT)
+                    category = np.argmax(category_test) - 1
+                    track_name = track.n_trk.item()
+                    id_no = track.n_trk.item()
+
+                    data.append(
+                        xr.Dataset(
+                            {
+                                "time_step": ("time", time_step_vector),
+                                "max_sustained_wind": ("time", max_sustained_wind),
+                                # "central_pressure": ("time", central_pressure),
+                                # "radius_max_wind": ("time", radius_max_wind),
+                                "environmental_pressure": ("time", env_pressure_vect),
+                                "basin": ("time", basin_vector),
+                            },
+                            coords={
+                                "time": ("time", time),
+                                "lat": ("time", lat),
+                                "lon": ("time", lon),
+                            },
+                            attrs={
+                                "max_sustained_wind_unit": "kn",
+                                "central_pressure_unit": "mb",
+                                "name": track_name,
+                                "sid": track_name,
+                                "orig_event_flag": False,
+                                "data_provider": "FAST",
+                                "id_no": id_no,
+                                "category": category,
+                            },
+                        )
+                    )
+
+        return cls(data)
+
     def write_hdf5(self, file_name, complevel=5):
         """Write TC tracks in NetCDF4-compliant HDF5 format.
 
@@ -2665,20 +2752,20 @@ def ibtracs_fit_param(explained, explanatory, year_range=(1980, 2019), order=1):
         return sm_results
 
 
-def ibtracs_track_agency(ds_sel):
+def ibtracs_track_agency(track):
     """Get preferred IBTrACS agency for each entry in the dataset.
 
     Parameters
     ----------
-    ds_sel : xarray.Dataset
+    track : xarray.Dataset
         Subselection of original IBTrACS NetCDF dataset.
 
     Returns
     -------
     agency_pref : list of str
         Names of IBTrACS agencies in order of preference.
     track_agency_ix : xarray.DataArray of ints
-        For each entry in `ds_sel`, the agency to use, given as an index into `agency_pref`.
+        For each entry in `track`, the agency to use, given as an index into `agency_pref`.
     """
     agency_pref = ["wmo"] + IBTRACS_AGENCIES
     agency_map = {a.encode("utf-8"): i for i, a in enumerate(agency_pref)}
@@ -2687,11 +2774,11 @@ def ibtracs_track_agency(ds_sel):
     )
     agency_map[b""] = agency_map[b"wmo"]
     agency_fun = lambda x: agency_map[x]
-    if "track_agency" not in ds_sel.data_vars.keys():
-        ds_sel["track_agency"] = ds_sel["wmo_agency"].where(
-            ds_sel["wmo_agency"] != b"", ds_sel["usa_agency"]
+    if "track_agency" not in track.data_vars.keys():
+        track["track_agency"] = track["wmo_agency"].where(
+            track["wmo_agency"] != b"", track["usa_agency"]
         )
-    track_agency_ix = xr.apply_ufunc(agency_fun, ds_sel["track_agency"], vectorize=True)
+    track_agency_ix = xr.apply_ufunc(agency_fun, track["track_agency"], vectorize=True)
     return agency_pref, track_agency_ix
 
 

climada/hazard/test/test_tc_tracks.py

@@ -19,6 +19,7 @@
 Test tc_tracks module.
 """
 
+import os
 import unittest
 from datetime import datetime as dt
 
@@ -48,6 +49,37 @@
 TEST_TRACKS_ANTIMERIDIAN = DATA_DIR.joinpath("tracks-antimeridian")
 TEST_TRACKS_LEGACY_HDF5 = DATA_DIR.joinpath("tctracks_hdf5_legacy.nc")
 
+TEST_TRACKS_FAST_dummy = xr.Dataset(
+    data_vars={
+        "lon_trks": (("n_trk", "time"), np.random.uniform(-180, 180, size=(20, 361))),
+        "lat_trks": (("n_trk", "time"), np.random.uniform(-90, 90, size=(20, 361))),
+        "u250_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "v250_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "u850_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "v850_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "v_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "m_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "vmax_trks": (("n_trk", "time"), np.random.randn(20, 361)),
+        "tc_month": (("n_trk",), np.random.randint(1, 13, size=20)),
+        "tc_basins": (
+            ("n_trk",),
+            np.random.choice(["AU", "EP", "NA", "NI", "SI", "SP", "WP"], size=20),
+        ),
+        "tc_years": (("n_trk",), np.full(20, 2025)),
+        "seeds_per_month": (
+            ("year", "basin", "month"),
+            np.random.randint(0, 5, size=(1, 7, 12)),
+        ),
+    },
+    coords={
+        "n_trk": np.arange(20),
+        "time": np.linspace(0, 1.296e6, 361),
+        "year": [2025],
+        "basin": ["AU", "EP", "NA", "NI", "SI", "SP", "WP"],
+        "month": np.arange(1, 13),
+    },
+)
+
 
 class TestIbtracs(unittest.TestCase):
     """Test reading and model of TC from IBTrACS files"""
@@ -609,6 +641,18 @@ def test_from_simulations_storm(self):
         tc_track = tc.TCTracks.from_simulations_storm(TEST_TRACK_STORM, years=[7])
         self.assertEqual(len(tc_track.data), 0)
 
+    def test_from_netcdf_fast(self):
+        """test the import of netcdf files from fast model"""
+
+        # create dummy .nc file to read, delate it in the end
+        file_path = DATA_DIR.joinpath("fast_test_tracks.nc")
+        TEST_TRACKS_FAST_dummy.to_netcdf(file_path)
+        tc_track = tc.TCTracks.from_netcdf_fast(file_path)
+        # test various instances
+
+        # remove file
+        os.remove(file_path)
+
     def test_to_geodataframe_points(self):
         """Conversion of TCTracks to GeoDataFrame using Points."""
         tc_track = tc.TCTracks.from_processed_ibtracs_csv(TEST_TRACK)