balck

Author: NicolasColombi

climada/hazard/tc_tracks.py

@@ -1620,11 +1620,44 @@ def from_netcdf(cls, folder_name):
                 data.append(track)
         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.
+
+        Parameters:
+        -----------
+        basin : str
+            Basin of generation of the TC
+        v_max : np.array
+            1D vector of maximal wind speed along the track
+
+        Returns:
+        --------
+        Pc : np.array
+            1D vector of central pressure along the track
+        """
+        a = 3.4
+        Pn = np.full(len(v_max), BASIN_ENV_PRESSURE[basin])
+        Pc = Pn - (v_max ** (1000 / 644)) / a
+
+        return Pc
+
+    @staticmethod
+    def compute_radius_max_winds():
+        pass
+
+    @staticmethod
+    def define_category_storm(self):
+        pass
+
     @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
@@ -1633,6 +1666,8 @@ def from_netcdf_fast(cls, folder_name):
         ----------
         folder_name : str
             Folder name from where to read files.
+        storm_id : int
+            Number of the simulated storm
 
         Returns:
         -------
@@ -1648,6 +1683,9 @@ 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())
 
@@ -1660,8 +1698,6 @@ def from_netcdf_fast(cls, folder_name):
                     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)
 
@@ -1673,17 +1709,19 @@ def from_netcdf_fast(cls, folder_name):
                         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()
+                    # Define central pressure
+                    central_pressure = TCTracks.compute_central_pressure(
+                        v_max=max_sustained_wind, basin=track.tc_basins.data.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),
+                                "central_pressure": ("time", central_pressure),
                                 "basin": ("time", basin_vector),
                             },
                             coords={
@@ -1692,12 +1730,10 @@ def from_netcdf_fast(cls, folder_name):
                                 "lon": ("time", lon),
                             },
                             attrs={
-                                "max_sustained_wind_unit": "kn",
-                                "central_pressure_unit": "mb",
-                                "name": track_name,
-                                "sid": track_name,
-                                "orig_event_flag": False,
+                                "max_sustained_wind_unit": "m/s",
+                                "central_pressure_unit": "hPa",
                                 "data_provider": "FAST",
+                                "orig_event_flag": False,
                                 "id_no": id_no,
                                 "category": category,
                             },

climada/hazard/test/test_tc_tracks.py

@@ -63,7 +63,7 @@
         "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),
+            np.full(20, "WP"),
         ),
         "tc_years": (("n_trk",), np.full(20, 2025)),
         "seeds_per_month": (
@@ -641,16 +641,42 @@ 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_from_netcdf_fast(self):
         """test the import of netcdf files from fast model"""
 
-        # create dummy .nc file to read, delate it in the end
+        # 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)
-        # test various instances
 
-        # remove file
+        expected_attributes = {
+            "max_sustained_wind_unit": "m/s",
+            "central_pressure": "hPa",
+            "data_provider": "FAST",
+            "orig_event_flag": False,
+            "id_no": 0,
+            "category": -1,
+        }
+
+        self.assertIsInstance(
+            tc_track, tc.TCTracks, "tc_track is not an instance of TCTracks"
+        )
+        self.assertIsInstance(
+            tc_track.data, list, "tc_track.data is not an instance of list"
+        )
+        self.assertIsInstance(
+            tc_track.data[0],
+            xr.Dataset,
+            "tc_track.data[0] not an instance of xarray.Dataset",
+        )
+        self.assertEqual(len(tc_track.data), 20)
+        self.assertEqual(tc_track.data[0].attrs, expected_attributes)
+        self.assertEqual(tc_track.data[0].environmental_pressure.data[0], 1005)
+        self.assertEqual(list(tc_track.data[0].coords.keys()), ["time", "lat", "lon"])
+
         os.remove(file_path)
 
     def test_to_geodataframe_points(self):