TCTracks: Option to read additional variables from IBTrACS (#728)

Add option to read additional variables from IBTrACS datasets
(e.g. "nature" or "stormtype")

Author: tovogt

CHANGELOG.md

@@ -41,6 +41,7 @@ Removed:
 - `Exposures.affected_total_value` now takes a hazard intensity threshold as argument. Affected values are only those for which at least one event exceeds the threshold. (previously, all exposures points with an assigned centroid were considered affected). By default the centroids are reassigned. [#702](https://github.com/CLIMADA-project/climada_python/pull/702) [#730](https://github.com/CLIMADA-project/climada_python/pull/730)
 - Add option to pass region ID to `LitPop.from_shape` [#720](https://github.com/CLIMADA-project/climada_python/pull/720)
 - Slightly improved performance on `LitPop`-internal computations [#720](https://github.com/CLIMADA-project/climada_python/pull/720)
+- Add option to read additional variables from IBTrACS when using `TCTracks.from_ibtracs_netcdf` [#728](https://github.com/CLIMADA-project/climada_python/pull/728)
 
 ### Fixed
 

climada/hazard/tc_tracks.py

@@ -189,6 +189,9 @@ class TCTracks():
         Computed during processing:
             - on_land (bool for each track position)
             - dist_since_lf (in km)
+        Additional data variables such as "nature" (specifiying, for each track position, whether a
+        system is a disturbance, tropical storm, post-transition extratropical storm etc.) might be
+        included, depending on the data source and on use cases.
     """
     def __init__(self,
                  data: Optional[List[xr.Dataset]] = None,
@@ -328,7 +331,7 @@ def read_ibtracs_netcdf(self, *args, **kwargs):
     def from_ibtracs_netcdf(cls, provider=None, rescale_windspeeds=True, storm_id=None,
                             year_range=None, basin=None, genesis_basin=None,
                             interpolate_missing=True, estimate_missing=False, correct_pres=False,
-                            discard_single_points=True,
+                            discard_single_points=True, additional_variables=None,
                             file_name='IBTrACS.ALL.v04r00.nc'):
         """Create new TCTracks object from IBTrACS databse.
 
@@ -438,6 +441,10 @@ def from_ibtracs_netcdf(cls, provider=None, rescale_windspeeds=True, storm_id=No
         file_name : str, optional
             Name of NetCDF file to be dowloaded or located at climada/data/system.
             Default: 'IBTrACS.ALL.v04r00.nc'
+        additional_variables : list of str, optional
+            If specified, additional IBTrACS data variables are extracted, such as "nature" or
+            "storm_speed". Only variables that are not agency-specific are supported.
+            Default: None.
 
         Returns
         -------
@@ -462,6 +469,9 @@ def from_ibtracs_netcdf(cls, provider=None, rescale_windspeeds=True, storm_id=No
                     f'Error while downloading {IBTRACS_URL}. Try to download it manually and '
                     f'put the file in {ibtracs_path}') from err
 
+        if additional_variables is None:
+            additional_variables = []
+
         ibtracs_ds = xr.open_dataset(ibtracs_path)
         ibtracs_date = ibtracs_ds.attrs["date_created"]
         if (np.datetime64('today') - np.datetime64(ibtracs_date)).item().days > 180:
@@ -576,7 +586,8 @@ def from_ibtracs_netcdf(cls, provider=None, rescale_windspeeds=True, storm_id=No
                             ibtracs_ds[tc_var].sel(storm=nonsingular_mask).interpolate_na(
                                 dim="date_time", method="linear"))
         ibtracs_ds = ibtracs_ds[['sid', 'name', 'basin', 'time', 'valid_t']
-                                + phys_vars + [f'{v}_agency' for v in phys_vars]]
+                                + additional_variables + phys_vars
+                                + [f'{v}_agency' for v in phys_vars]]
 
         if estimate_missing:
             ibtracs_ds['pres'][:] = _estimate_pressure(
@@ -685,28 +696,37 @@ def from_ibtracs_netcdf(cls, provider=None, rescale_windspeeds=True, storm_id=No
                     "{}({})".format(v, track_ds[f'{v}_agency'].astype(str).item())
                     for v in phys_vars)
 
-            all_tracks.append(xr.Dataset({
-                'time_step': ('time', track_ds.time_step.data),
+            data_vars = {
                 'radius_max_wind': ('time', track_ds.rmw.data),
                 'radius_oci': ('time', track_ds.roci.data),
                 'max_sustained_wind': ('time', track_ds.wind.data),
                 'central_pressure': ('time', track_ds.pres.data),
                 'environmental_pressure': ('time', track_ds.poci.data),
-                'basin': ('time', track_ds.basin.data.astype("<U2")),
-            }, coords={
-                'time': track_ds.time.dt.round('s').data,
+            }
+            coords = {
+                'time': ('time', track_ds.time.dt.round('s').data),
                 'lat': ('time', track_ds.lat.data),
                 'lon': ('time', track_ds.lon.data),
-            }, attrs={
+            }
+            attrs = {
                 'max_sustained_wind_unit': 'kn',
                 'central_pressure_unit': 'mb',
-                'name': track_ds.name.astype(str).item(),
-                'sid': track_ds.sid.astype(str).item(),
                 'orig_event_flag': True,
                 'data_provider': provider_str,
-                'id_no': track_ds.id_no.item(),
                 'category': category[i_track],
-            }))
+            }
+            # automatically assign the remaining variables as attributes or data variables
+            for varname in ["time_step", "basin", "name", "sid", "id_no"] + additional_variables:
+                values = track_ds[varname].data
+                if track_ds[varname].dtype.kind == "S":
+                    # This converts the `bytes` (dtype "|S*") in IBTrACS to the more common `str`
+                    # objects (dtype "<U*") that we use in CLIMADA.
+                    values = values.astype(str)
+                if values.ndim == 0:
+                    attrs[varname] = values.item()
+                else:
+                    data_vars[varname] = ('time', values)
+            all_tracks.append(xr.Dataset(data_vars, coords=coords, attrs=attrs))
         if last_perc != 100:
             LOGGER.info("Progress: 100%")
         if len(all_tracks) == 0:
@@ -1382,8 +1402,10 @@ def from_hdf5(cls, file_name):
         for i in track_no:
             track = ds_dict[f'track{i}']
             track.attrs['orig_event_flag'] = bool(track.attrs['orig_event_flag'])
-            # when writing '<U2' and reading in again, xarray reads as dtype 'object'. undo this:
-            track['basin'] = track['basin'].astype('<U2')
+            # when writing '<U*' and reading in again, xarray reads as dtype 'object'. undo this:
+            for varname in track.data_vars:
+                if track[varname].dtype == "object":
+                    track[varname] = track[varname].astype(str)
             data.append(track)
         return cls(data)
 
@@ -1499,7 +1521,9 @@ def _one_interp_data(track, time_step_h, land_geom=None):
             time_step = pd.tseries.frequencies.to_offset(pd.Timedelta(hours=time_step_h)).freqstr
             track_int = track.resample(time=time_step, skipna=True)\
                              .interpolate('linear')
-            track_int['basin'] = track.basin.resample(time=time_step).nearest()
+            for var in track.data_vars:
+                if "time" in track[var].dims and track[var].dtype.kind != "f":
+                    track_int[var] = track[var].resample(time=time_step).nearest()
             track_int['time_step'][:] = time_step_h
             lon_int = lon.resample(time=time_step).interpolate(method)
             lon_int[lon_int > 180] -= 360

climada/hazard/test/test_tc_tracks.py

@@ -256,18 +256,56 @@ def test_ibtracs_discard_single_points(self):
                     break
         self.assertTrue(passed)
 
+    def test_ibtracs_additional_variables(self):
+        """Check additional_variables option"""
+        # this is the complete list (as of May 2023) of variables in IBTrACS that are not
+        # agency-specific and that are not already considered by other parts of
+        # `from_ibtracs_netcdf`:
+        addtl_vars = [
+            'numobs', 'season', 'number', 'subbasin', 'name', 'source_usa', 'source_jma',
+            'source_cma', 'source_hko', 'source_new', 'source_reu', 'source_bom', 'source_nad',
+            'source_wel', 'source_td5', 'source_td6', 'source_ds8', 'source_neu', 'source_mlc',
+            'iso_time', 'nature', 'wmo_wind', 'wmo_pres', 'wmo_agency', 'track_type',
+            'main_track_sid', 'dist2land', 'landfall', 'iflag', 'storm_speed', 'storm_dir',
+        ]
+        tc_track = tc.TCTracks.from_ibtracs_netcdf(
+            storm_id='2017242N16333',
+            additional_variables=addtl_vars,
+        )
+        track_ds = tc_track.get_track()
+        for v in addtl_vars:
+            self.assertIn(v, list(track_ds.data_vars) + list(track_ds.attrs))
+        self.assertEqual(track_ds.attrs["numobs"], 123)
+        self.assertEqual(track_ds.attrs["season"], 2017)
+        self.assertEqual(track_ds["nature"].values[0], "TS")
+        self.assertEqual(track_ds["nature"].values[-1], "DS")
+        self.assertEqual(track_ds["subbasin"].values[0], "NA")
+        self.assertEqual(track_ds["subbasin"].values[58], "CS")
+        self.assertEqual(track_ds["dist2land"].values[0], 1020)
+        self.assertEqual(track_ds["dist2land"].values[-1], 0)
+        self.assertEqual(track_ds["storm_speed"].values[0], 13.0)
+        self.assertEqual(track_ds["storm_speed"].values[5], 11.0)
+        self.assertEqual(track_ds["storm_speed"].values[-1], 8.0)
+
 class TestIO(unittest.TestCase):
     """Test reading of tracks from files of different formats"""
     def test_netcdf_io(self):
         """Test writting and reading netcdf4 TCTracks instances"""
         path = DATA_DIR.joinpath("tc_tracks_nc")
         path.mkdir(exist_ok=True)
         tc_track = tc.TCTracks.from_ibtracs_netcdf(
-            provider='usa', storm_id='1988234N13299', estimate_missing=True)
-        tc_track.write_netcdf(str(path))
+            provider='usa', storm_id='1988234N13299', estimate_missing=True,
+            additional_variables=["numobs", "storm_speed", "nature"],
+        )
+        tc_track.write_netcdf(path)
+        tc_read = tc.TCTracks.from_netcdf(path)
 
-        tc_read = tc.TCTracks.from_netcdf(str(path))
-        self.assertEqual(tc_track.get_track().sid, tc_read.get_track().sid)
+        ds_read = tc_read.get_track()
+        ds_write = tc_track.get_track()
+        self.assertEqual(ds_write.attrs, ds_read.attrs)
+        self.assertEqual(sorted(ds_write.data_vars), sorted(ds_read.data_vars))
+        for v in ds_write.data_vars:
+            np.testing.assert_array_equal(ds_write[v], ds_read[v])
 
     def test_read_legacy_netcdf(self):
         """Test reading from NetCDF files with legacy basin attributes"""
@@ -281,7 +319,7 @@ def test_read_legacy_netcdf(self):
             np.testing.assert_array_equal(tr.basin, "SP")
 
     def test_hdf5_io(self):
-        """Test writting and reading hdf5 TCTracks instances"""
+        """Test writing and reading hdf5 TCTracks instances"""
         path = DATA_DIR.joinpath("tc_tracks.h5")
         tc_track = tc.TCTracks.from_ibtracs_netcdf(
             provider='usa', year_range=(1993, 1994), basin='EP', estimate_missing=True)