add .nc test file and update code

NicolasColombi · NicolasColombi · commit c1f8847dc255 · 2025-01-16T11:30:31.000+01:00
diff --git a/climada/hazard/tc_tracks.py b/climada/hazard/tc_tracks.py
@@ -52,10 +52,9 @@
 from shapely.geometry import LineString, MultiLineString, Point
 from sklearn.metrics import DistanceMetric
 
-import climada.hazard.tc_tracks_synth
 import climada.util.coordinates as u_coord
 import climada.util.plot as u_plot
-from climada.hazard import Centroids
+from climada.hazard import Centroids, tc_tracks, tc_tracks_synth
 
 # climada dependencies
 from climada.util import ureg
@@ -1621,46 +1620,50 @@ def from_netcdf(cls, folder_name):
         return cls(data)
 
     @staticmethod
-    def compute_central_pressure(basin, v_max):
-        """Compute central pressure of tropical cyclone given the maximal
-        wind speed of the storm. Method needed to load tracks from_netcdf_fast.
-
-        Holland, Greg. (2008). A Revised Hurricane Pressure Wind Model.
-        Monthly Weather Review - MON WEATHER REV. 136. 10.1175/2008MWR2395.1.
+    def define_tc_category_fast(max_sust_wind: np.array, units: str = "kn") -> int:
+        """Define category of the tropical cyclone according to Saffir-Simpson scale.
 
         Parameters:
-        -----------
-        basin : str
-            Basin of generation of the TC
-        v_max : np.array
-            1D vector of maximal wind speed along the track
+        ----------
+        max_wind : str
+            Maximal sustained wind speed
+        units: str
+            Wind speed units, kn or m/s. Default: kn
 
         Returns:
-        --------
-        Pc : np.array
-            1D vector of central pressure along the track
+        -------
+        category : int
+           -1: "Tropical Depression",
+            0: "Tropical Storm",
+            1: "Hurricane Cat. 1",
+            2: "Hurricane Cat. 2",
+            3: "Hurricane Cat. 3",
+            4: "Hurricane Cat. 4",
+            5: "Hurricane Cat. 5",
         """
-        a = 3.4
-        Pn = np.full(len(v_max), BASIN_ENV_PRESSURE[basin])
-        Pc = Pn - (v_max ** (1000 / 644)) / a
 
-        return Pc
+        max_sust_wind = max_sust_wind.astype(
+            float
+        )  # avoid casting errors if max_sust_wind is int
+        max_sust_wind *= 1.943844 if units == "m/s" else 1
 
-    @staticmethod
-    def compute_radius_max_winds():
-        pass
+        max_wind = np.nanmax(max_sust_wind)
+        category_test = np.full(len(SAFFIR_SIM_CAT), max_wind) < np.array(
+            SAFFIR_SIM_CAT
+        )
+        category = np.argmax(category_test) - 1
 
-    @staticmethod
-    def define_category_storm(self):
-        pass
+        return category
 
     @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?
+    def from_fast(cls, folder_name: str):
+        """Create a new TCTracks object from NetCDF files generated by the FAST model, modifying
+        the xr.array structure to ensure compatibility with CLIMADA, and calculating the central
+        pressure and radius of maximum wind.
+
+        Model GitHub Repository: https://github.com/linjonathan/tropical_cyclone_risk?
                             tab=readme-ov-file
-        Publication: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2023MS003686
+        Model Publication: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2023MS003686
 
         Parameters:
         ----------
@@ -1672,7 +1675,7 @@ def from_netcdf_fast(cls, folder_name):
         Returns:
         -------
         tracks : TCTracks
-            TCTracks obecjt with tracks data from the given directory of NetCDF files.
+            TCTracks object with tracks data from the given directory of NetCDF files.
         """
 
         file_tr = get_file_names(folder_name)
@@ -1683,45 +1686,53 @@ def from_netcdf_fast(cls, folder_name):
                 continue
             with xr.open_dataset(file) as ds:
                 for i in ds.n_trk:
-                    # if storm_id:
-                    #     i == storm_id
 
                     # Select track
-                    track = ds.sel(n_trk=i.item())
+                    track = ds.sel(n_trk=i)
 
                     # Define coordinates
                     lat = track.lat_trks.data
                     lon = track.lon_trks.data
                     time = track.time.data
 
+                    # Convert time to pandas Datetime "yyyy.mm.dd"
+                    reference_time = (
+                        f"{track.tc_years.item()}-{int(track.tc_month.item())}-01"
+                    )
+                    time = pd.to_datetime(time, unit="s", origin=reference_time).astype(
+                        "datetime64[s]"
+                    )
+
                     # Define variables
                     time_step_vector = np.full(time.shape[0], track.time.data[1])
-                    max_sustained_wind = track.v_trks.data
+                    max_sustained_wind_ms = track.v_trks.data
+                    max_sustained_wind_knots = max_sustained_wind_ms * 1.943844
                     basin_vector = np.full(time.shape[0], track.tc_basins.data.item())
                     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
-                    id_no = track.n_trk.item()
-                    # Define central pressure
-                    central_pressure = TCTracks.compute_central_pressure(
-                        v_max=max_sustained_wind, basin=track.tc_basins.data.item()
+                    cen_pres_missing = np.full(lat.shape, np.nan)
+                    rmw_missing = np.full(lat.shape, np.nan)
+                    cen_pres = tc_tracks._estimate_pressure(
+                        cen_pres_missing, lat, lon, max_sustained_wind_knots
                     )
+                    rmw = tc_tracks.estimate_rmw(rmw_missing, cen_pres)
+
+                    # Define attributes
+                    category = tc_tracks.TCTracks.define_tc_category_fast(
+                        max_sustained_wind_knots
+                    )
+                    id_no = track.n_trk.item()
+                    track_name = f"storm_{id_no}"
 
                     data.append(
                         xr.Dataset(
                             {
                                 "time_step": ("time", time_step_vector),
-                                "max_sustained_wind": ("time", max_sustained_wind),
+                                "max_sustained_wind": ("time", max_sustained_wind_ms),
+                                "central_pressure": ("time", cen_pres),
+                                "radius_max_wind": ("time", rmw),
                                 "environmental_pressure": ("time", env_pressure_vect),
-                                "central_pressure": ("time", central_pressure),
                                 "basin": ("time", basin_vector),
                             },
                             coords={
@@ -1732,8 +1743,10 @@ def from_netcdf_fast(cls, folder_name):
                             attrs={
                                 "max_sustained_wind_unit": "m/s",
                                 "central_pressure_unit": "hPa",
+                                "name": track_name,
+                                "sid": id_no,
+                                "orig_event_flag": True,
                                 "data_provider": "FAST",
-                                "orig_event_flag": False,
                                 "id_no": id_no,
                                 "category": category,
                             },
diff --git a/climada/hazard/test/data/FAST_test_tracks.nc b/climada/hazard/test/data/FAST_test_tracks.nc
diff --git a/climada/hazard/test/test_tc_tracks.py b/climada/hazard/test/test_tc_tracks.py
@@ -45,41 +45,11 @@
 TEST_TRACK_EMANUEL = DATA_DIR.joinpath("emanuel_test_tracks.mat")
 TEST_TRACK_EMANUEL_CORR = DATA_DIR.joinpath("temp_mpircp85cal_full.mat")
 TEST_TRACK_CHAZ = DATA_DIR.joinpath("chaz_test_tracks.nc")
+TEST_TRACK_FAST = DATA_DIR.joinpath("FAST_test_tracks.nc")
 TEST_TRACK_STORM = DATA_DIR.joinpath("storm_test_tracks.txt")
 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.full(20, "WP"),
-        ),
-        "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"""
@@ -663,24 +633,32 @@ 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_compute_central_pressure(self):
-        pass
+    def test_define_tc_category_fast(self):
+        """test that the correct category is assigned to a tc."""
+
+        max_wind = np.array([20, 72, 36, 50])  # knots
+        category1 = tc.TCTracks.define_tc_category_fast(max_wind)
+        category2 = tc.TCTracks.define_tc_category_fast(
+            max_sust_wind=max_wind, units="m/s"
+        )
+        self.assertEqual(category1, 1)
+        self.assertEqual(category2, 5)
 
-    def test_from_netcdf_fast(self):
-        """test the import of netcdf files from fast model"""
+    def test_from_fast(self):
+        """test the correct import of netcdf files from fast model and the conversion to a
+        different xr.array structure compatible with CLIMADA."""
 
-        # create dummy .nc file to be read
-        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)
+        tc_track = tc.TCTracks.from_fast(TEST_TRACK_FAST)
 
         expected_attributes = {
             "max_sustained_wind_unit": "m/s",
-            "central_pressure": "hPa",
+            "central_pressure_unit": "hPa",
+            "name": "storm_0",
+            "sid": 0,
+            "orig_event_flag": True,
             "data_provider": "FAST",
-            "orig_event_flag": False,
             "id_no": 0,
-            "category": -1,
+            "category": 0,
         }
 
         self.assertIsInstance(
@@ -694,13 +672,11 @@ def test_from_netcdf_fast(self):
             xr.Dataset,
             "tc_track.data[0] not an instance of xarray.Dataset",
         )
-        self.assertEqual(len(tc_track.data), 20)
+        self.assertEqual(len(tc_track.data), 5)
         self.assertEqual(tc_track.data[0].attrs, expected_attributes)
-        self.assertEqual(tc_track.data[0].environmental_pressure.data[0], 1005)
+        self.assertEqual(tc_track.data[0].environmental_pressure.data[0], 1010)
         self.assertEqual(list(tc_track.data[0].coords.keys()), ["time", "lat", "lon"])
 
-        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)
