diff --git a/CHANGELOG.md b/CHANGELOG.md index 070763bb16..90d02c85cb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -27,12 +27,14 @@ CLIMADA tutorials. [#872](https://github.com/CLIMADA-project/climada_python/pull - Changed module structure: `climada.hazard.Hazard` has been split into the modules `base`, `io` and `plot` [#871](https://github.com/CLIMADA-project/climada_python/pull/871) - `Impact.from_hdf5` now calls `str` on `event_name` data that is not strings, and issue a warning then [#894](https://github.com/CLIMADA-project/climada_python/pull/894) - `Impact.write_hdf5` now throws an error if `event_name` is does not contain strings exclusively [#894](https://github.com/CLIMADA-project/climada_python/pull/894) +- Split `climada.hazard.trop_cyclone` module into smaller submodules without affecting module usage [#911](https://github.com/CLIMADA-project/climada_python/pull/911) ### Fixed - Avoid an issue where a Hazard subselection would have a fraction matrix with only zeros as entries by throwing an error [#866](https://github.com/CLIMADA-project/climada_python/pull/866) - Allow downgrading the Python bugfix version to improve environment compatibility [#900](https://github.com/CLIMADA-project/climada_python/pull/900) - Fix broken links in `CONTRIBUTING.md` [#900](https://github.com/CLIMADA-project/climada_python/pull/900) +- When writing `TCTracks` to NetCDF, only apply compression to `float` or `int` data types. This fixes a downstream issue, see [climada_petals#135](https://github.com/CLIMADA-project/climada_petals/issues/135) [#911](https://github.com/CLIMADA-project/climada_python/pull/911) ### Added diff --git a/climada/hazard/tc_tracks.py b/climada/hazard/tc_tracks.py index 89a44b3744..7023fe3290 100644 --- a/climada/hazard/tc_tracks.py +++ b/climada/hazard/tc_tracks.py @@ -1383,7 +1383,10 @@ def write_hdf5(self, file_name, complevel=5): df_attrs = pd.DataFrame([t.attrs for t in data], index=ds_combined["storm"].to_series()) ds_combined = xr.merge([ds_combined, df_attrs.to_xarray()]) - encoding = {v: dict(zlib=True, complevel=complevel) for v in ds_combined.data_vars} + encoding = { + v: dict(zlib=_zlib_from_dataarray(ds_combined[v]), complevel=complevel) + for v in ds_combined.data_vars + } LOGGER.info('Writing %d tracks to %s', self.size, file_name) ds_combined.to_netcdf(file_name, encoding=encoding) @@ -2435,3 +2438,18 @@ def set_category(max_sus_wind, wind_unit='kn', saffir_scale=None): return (np.argwhere(max_wind < saffir_scale) - 1)[0][0] except IndexError: return -1 + +def _zlib_from_dataarray(data_var: xr.DataArray) -> bool: + """Return true if data_var is of numerical type, return False otherwise + + Parameters + ---------- + data_var : xr.DataArray + + Returns + ------- + bool + """ + if np.issubdtype(data_var.dtype, float) or np.issubdtype(data_var.dtype, int): + return True + return False diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index b58bf268eb..41ab3a81fb 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -26,25 +26,21 @@ import numpy as np from scipy import sparse -import xarray as xr +import climada.hazard.test as hazard_test from climada.util import ureg from climada.test import get_test_file from climada.hazard.tc_tracks import TCTracks -from climada.hazard.trop_cyclone import ( - TropCyclone, get_close_centroids, _vtrans, _B_holland_1980, _bs_holland_2008, - _v_max_s_holland_2008, _x_holland_2010, _stat_holland_1980, _stat_holland_2010, _stat_er_2011, - tctrack_to_si, MBAR_TO_PA, KM_TO_M, H_TO_S, -) +from climada.hazard.trop_cyclone.trop_cyclone import ( + TropCyclone, ) from climada.hazard.centroids.centr import Centroids -import climada.hazard.test as hazard_test - DATA_DIR = Path(hazard_test.__file__).parent.joinpath('data') TEST_TRACK = DATA_DIR.joinpath("trac_brb_test.csv") TEST_TRACK_SHORT = DATA_DIR.joinpath("trac_short_test.csv") + CENTR_TEST_BRB = Centroids.from_hdf5(get_test_file('centr_test_brb', file_format='hdf5')) @@ -294,236 +290,6 @@ def test_two_files_pass(self): self.assertEqual(tc_haz.fraction.nonzero()[0].size, 0) self.assertEqual(tc_haz.intensity.nonzero()[0].size, 0) -class TestWindfieldHelpers(unittest.TestCase): - """Test helper functions of TC wind field model""" - - def test_get_close_centroids_pass(self): - """Test get_close_centroids function.""" - si_track = xr.Dataset({ - "lat": ("time", np.array([0, -0.5, 0])), - "lon": ("time", np.array([0.9, 2, 3.2])), - }, attrs={"mid_lon": 0.0}) - centroids = np.array([ - [0, -0.2], [0, 0.9], [-1.1, 1.2], [1, 2.1], [0, 4.3], [0.6, 3.8], [0.9, 4.1], - ]) - centroids_close, mask_close, mask_close_alongtrack = ( - get_close_centroids(si_track, centroids, 112.0) - ) - self.assertEqual(centroids_close.shape[0], mask_close.sum()) - self.assertEqual(mask_close_alongtrack.shape[0], si_track.sizes["time"]) - self.assertEqual(mask_close_alongtrack.shape[1], centroids_close.shape[0]) - np.testing.assert_equal(mask_close_alongtrack.any(axis=0), True) - np.testing.assert_equal(mask_close, np.array( - [False, True, True, False, False, True, False] - )) - np.testing.assert_equal(mask_close_alongtrack, np.array([ - [True, False, False], - [False, True, False], - [False, False, True], - ])) - np.testing.assert_equal(centroids_close, centroids[mask_close]) - - # example where antimeridian is crossed - si_track = xr.Dataset({ - "lat": ("time", np.linspace(-10, 10, 11)), - "lon": ("time", np.linspace(170, 200, 11)), - }, attrs={"mid_lon": 180.0}) - centroids = np.array([[-11, 169], [-7, 176], [4, -170], [10, 170], [-10, -160]]) - centroids_close, mask_close, mask_close_alongtrack = ( - get_close_centroids(si_track, centroids, 600.0) - ) - self.assertEqual(centroids_close.shape[0], mask_close.sum()) - self.assertEqual(mask_close_alongtrack.shape[0], si_track.sizes["time"]) - self.assertEqual(mask_close_alongtrack.shape[1], centroids_close.shape[0]) - np.testing.assert_equal(mask_close_alongtrack.any(axis=0), True) - np.testing.assert_equal(mask_close, np.array([True, True, True, False, False])) - np.testing.assert_equal(centroids_close, np.array([ - # the longitudinal coordinate of the third centroid is normalized - [-11, 169], [-7, 176], [4, 190], - ])) - - def test_B_holland_1980_pass(self): - """Test _B_holland_1980 function.""" - si_track = xr.Dataset({ - "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), - "vgrad": ("time", [35, 40]), - "rho_air": ("time", [1.15, 1.15]) - }) - _B_holland_1980(si_track) - np.testing.assert_array_almost_equal(si_track["hol_b"], [2.5, 1.667213]) - - si_track = xr.Dataset({ - "pdelta": ("time", MBAR_TO_PA * np.array([4.74, 15, 30, 40])), - "vmax": ("time", [np.nan, 22.5, 25.4, 42.5]), - "rho_air": ("time", [1.2, 1.2, 1.2, 1.2]) - }) - _B_holland_1980(si_track, gradient_to_surface_winds=0.9) - np.testing.assert_allclose(si_track["hol_b"], [np.nan, 1.101, 0.810, 1.473], atol=1e-3) - - def test_bs_holland_2008_pass(self): - """Test _bs_holland_2008 function. Compare to MATLAB reference.""" - si_track = xr.Dataset({ - "tstep": ("time", H_TO_S * np.array([1.0, 1.0, 1.0])), - "lat": ("time", [12.299999504631234, 12.299999504631343, 12.299999279463769]), - "pdelta": ("time", MBAR_TO_PA * np.array([4.74, 4.73, 4.73])), - "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 1005.2633, 1005.2682])), - "vtrans_norm": ("time", [np.nan, 5.241999541820597, 5.123882725120426]), - }) - _bs_holland_2008(si_track) - np.testing.assert_allclose( - si_track["hol_b"], [np.nan, 1.27, 1.27], atol=1e-2 - ) - - def test_v_max_s_holland_2008_pass(self): - """Test _v_max_s_holland_2008 function.""" - # Numbers analogous to test_B_holland_1980_pass - si_track = xr.Dataset({ - "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), - "hol_b": ("time", [2.5, 1.67]), - "rho_air": ("time", [1.15, 1.15]), - }) - _v_max_s_holland_2008(si_track) - np.testing.assert_array_almost_equal(si_track["vmax"], [34.635341, 40.033421]) - - def test_holland_2010_pass(self): - """Test Holland et al. 2010 wind field model.""" - # The parameter "x" is designed to be exactly 0.5 inside the radius of max wind (RMW) and - # to increase or decrease linearly outside of it in radial direction. - # - # An increase (decrease) of "x" outside of the RMW is for cases where the max wind is very - # high (low), but the RMW is still comparably large (small). This means, wind speeds need - # to decay very sharply (only moderately) outside of the RMW to reach the low prescribed - # peripheral wind speeds. - # - # The "hol_b" parameter tunes the meaning of a "comparably" large or small RMW. - si_track = xr.Dataset({ - # four test cases: - # - low vmax, moderate RMW: x decreases moderately - # - large hol_b: x decreases sharply - # - very low vmax: x decreases so much, it needs to be clipped at 0 - # - large vmax, large RMW: x increases - "rad": ("time", KM_TO_M * np.array([75, 75, 75, 90])), - "vmax": ("time", [35.0, 35.0, 16.0, 90.0]), - "hol_b": ("time", [1.75, 2.5, 1.9, 1.6]), - }) - d_centr = KM_TO_M * np.array([ - # first column is for locations within the storm eye - # second column is for locations at or close to the radius of max wind - # third column is for locations outside the storm eye - # fourth column is for locations exactly at the peripheral radius - # fifth column is for locations outside the peripheral radius - [0., 75, 220, 300, 490], - [30, 74, 170, 300, 501], - [21, 76, 230, 300, 431], - [32, 91, 270, 300, 452], - ], dtype=float) - close_centr = np.array([ - # note that we set one of these to "False" for testing - [True, True, True, True, True], - [True, True, True, True, False], - [True, True, True, True, True], - [True, True, True, True, True], - ], dtype=bool) - hol_x = _x_holland_2010(si_track, d_centr, close_centr) - np.testing.assert_array_almost_equal(hol_x, [ - [0.5, 0.500000, 0.485077, 0.476844, 0.457291], - [0.5, 0.500000, 0.410997, 0.400000, 0.000000], - [0.5, 0.497620, 0.400000, 0.400000, 0.400000], - [0.5, 0.505022, 1.403952, 1.554611, 2.317948], - ]) - - v_ang_norm = _stat_holland_2010(si_track, d_centr, close_centr, hol_x) - np.testing.assert_allclose(v_ang_norm, [ - # first column: converge to 0 when approaching storm eye - # second column: vmax at RMW - # fourth column: peripheral speed (17 m/s) at peripheral radius (unless x is clipped!) - [ 0.000000, 35.000000, 21.181497, 17.000000, 12.1034610], - [ 1.296480, 34.990037, 21.593755, 12.891313, 0.0000000], - [ 0.321952, 15.997500, 9.712006, 8.087240, 6.2289690], - [24.823469, 89.992938, 24.381965, 17.000000, 1.9292020], - ], atol=1e-6) - - def test_stat_holland_1980(self): - """Test _stat_holland_1980 function. Compare to MATLAB reference.""" - d_centr = KM_TO_M * np.array([ - [299.4501244109841, 291.0737897183741, 292.5441003235722, 40.665454622610511], - [293.6067129546862, 1000.0, 298.2652319413182, 70.0], - ]) - si_track = xr.Dataset({ - "rad": ("time", KM_TO_M * np.array([40.665454622610511, 75.547902916671745])), - "hol_b": ("time", [1.486076257880692, 1.265551666104679]), - "pdelta": ("time", MBAR_TO_PA * np.array([39.12, 4.73])), - "lat": ("time", [-14.089110370469488, 12.299999279463769]), - "cp": ("time", [3.54921922e-05, 3.10598285e-05]), - "rho_air": ("time", [1.15, 1.15]), - }) - mask = np.array([[True, True, True, True], [True, False, True, True]], dtype=bool) - v_ang_norm = _stat_holland_1980(si_track, d_centr, mask) - np.testing.assert_allclose( - v_ang_norm, [[11.28, 11.68, 11.61, 42.41], [5.38, 0, 5.28, 12.76]], atol=1e-2, - ) - - # without Coriolis force, values are higher, esp. far away from the center: - v_ang_norm = _stat_holland_1980(si_track, d_centr, mask, cyclostrophic=True) - np.testing.assert_allclose( - v_ang_norm, [[15.72, 16.04, 15.98, 43.13], [8.84, 0, 8.76, 13.81]], atol=1e-2, - ) - - d_centr = np.array([[], []]) - mask = np.ones_like(d_centr, dtype=bool) - v_ang_norm = _stat_holland_1980(si_track, d_centr, mask) - np.testing.assert_array_equal(v_ang_norm, np.array([[], []])) - - def test_er_2011_pass(self): - """Test Emanuel and Rotunno 2011 wind field model.""" - # test at centroids within and outside of radius of max wind - d_centr = KM_TO_M * np.array([[35, 70, 75, 220], [30, 150, 1000, 300]], dtype=float) - si_track = xr.Dataset({ - "rad": ("time", KM_TO_M * np.array([75.0, 40.0])), - "vmax": ("time", [35.0, 40.0]), - "lat": ("time", [20.0, 27.0]), - "cp": ("time", [4.98665369e-05, 6.61918149e-05]), - }) - mask = np.array([[True, True, True, True], [True, False, True, True]], dtype=bool) - v_ang_norm = _stat_er_2011(si_track, d_centr, mask) - np.testing.assert_array_almost_equal(v_ang_norm, - [[28.258025, 36.782418, 36.869995, 22.521237], - [39.670883, 0, 3.300626, 10.827206]]) - - def test_vtrans_pass(self): - """Test _vtrans function. Compare to MATLAB reference.""" - tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK) - tc_track.equal_timestep() - track_ds = tc_track.data[0] - - si_track = tctrack_to_si(track_ds) - _vtrans(si_track) - - to_kn = (1.0 * ureg.meter / ureg.second).to(ureg.knot).magnitude - - self.assertEqual(si_track["vtrans_norm"].size, track_ds["time"].size) - self.assertEqual(si_track["vtrans_norm"].values[0], 0) - self.assertAlmostEqual(si_track["vtrans_norm"].values[1] * to_kn, 10.191466246) - - def testtctrack_to_si(self): - """ Test tctrack_to_si should create the same vmax output independent of the input unit """ - tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK_SHORT).data[0] - - tc_track_kmph = tc_track.copy(deep=True) - tc_track_kmph['max_sustained_wind'] *= ( - (1.0 * ureg.knot).to(ureg.km / ureg.hour).magnitude - ) - tc_track_kmph.attrs['max_sustained_wind_unit'] = 'km/h' - - si_track = tctrack_to_si(tc_track) - si_track_from_kmph = tctrack_to_si(tc_track_kmph) - - np.testing.assert_array_almost_equal(si_track["vmax"], si_track_from_kmph["vmax"]) - - tc_track.attrs['max_sustained_wind_unit'] = 'elbows/fortnight' - with self.assertRaises(ValueError): - tctrack_to_si(tc_track) - class TestClimateSce(unittest.TestCase): def test_apply_criterion_track(self): @@ -682,7 +448,6 @@ def tearDown(self): if __name__ == "__main__": TESTS = unittest.TestLoader().loadTestsFromTestCase(TestReader) - TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestWindfieldHelpers)) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestClimateSce)) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestDumpReloadCycle)) unittest.TextTestRunner(verbosity=2).run(TESTS) diff --git a/climada/hazard/test/test_trop_cyclone_windfields.py b/climada/hazard/test/test_trop_cyclone_windfields.py new file mode 100644 index 0000000000..f91ac075ab --- /dev/null +++ b/climada/hazard/test/test_trop_cyclone_windfields.py @@ -0,0 +1,247 @@ +import unittest + +import numpy as np +import xarray as xr + +from climada.hazard import TCTracks +from climada.hazard.test.test_trop_cyclone import TEST_TRACK, TEST_TRACK_SHORT +from climada.hazard.trop_cyclone.trop_cyclone_windfields import get_close_centroids, MBAR_TO_PA, _B_holland_1980, H_TO_S, \ + _bs_holland_2008, _v_max_s_holland_2008, KM_TO_M, _x_holland_2010, _stat_holland_2010, _stat_holland_1980, \ + _stat_er_2011, tctrack_to_si, _vtrans +from climada.util import ureg + + +class TestWindfieldHelpers(unittest.TestCase): + """Test helper functions of TC wind field model""" + + def test_get_close_centroids_pass(self): + """Test get_close_centroids function.""" + si_track = xr.Dataset({ + "lat": ("time", np.array([0, -0.5, 0])), + "lon": ("time", np.array([0.9, 2, 3.2])), + }, attrs={"mid_lon": 0.0}) + centroids = np.array([ + [0, -0.2], [0, 0.9], [-1.1, 1.2], [1, 2.1], [0, 4.3], [0.6, 3.8], [0.9, 4.1], + ]) + centroids_close, mask_close, mask_close_alongtrack = ( + get_close_centroids(si_track, centroids, 112.0) + ) + self.assertEqual(centroids_close.shape[0], mask_close.sum()) + self.assertEqual(mask_close_alongtrack.shape[0], si_track.sizes["time"]) + self.assertEqual(mask_close_alongtrack.shape[1], centroids_close.shape[0]) + np.testing.assert_equal(mask_close_alongtrack.any(axis=0), True) + np.testing.assert_equal(mask_close, np.array( + [False, True, True, False, False, True, False] + )) + np.testing.assert_equal(mask_close_alongtrack, np.array([ + [True, False, False], + [False, True, False], + [False, False, True], + ])) + np.testing.assert_equal(centroids_close, centroids[mask_close]) + + # example where antimeridian is crossed + si_track = xr.Dataset({ + "lat": ("time", np.linspace(-10, 10, 11)), + "lon": ("time", np.linspace(170, 200, 11)), + }, attrs={"mid_lon": 180.0}) + centroids = np.array([[-11, 169], [-7, 176], [4, -170], [10, 170], [-10, -160]]) + centroids_close, mask_close, mask_close_alongtrack = ( + get_close_centroids(si_track, centroids, 600.0) + ) + self.assertEqual(centroids_close.shape[0], mask_close.sum()) + self.assertEqual(mask_close_alongtrack.shape[0], si_track.sizes["time"]) + self.assertEqual(mask_close_alongtrack.shape[1], centroids_close.shape[0]) + np.testing.assert_equal(mask_close_alongtrack.any(axis=0), True) + np.testing.assert_equal(mask_close, np.array([True, True, True, False, False])) + np.testing.assert_equal(centroids_close, np.array([ + # the longitudinal coordinate of the third centroid is normalized + [-11, 169], [-7, 176], [4, 190], + ])) + + def test_B_holland_1980_pass(self): + """Test _B_holland_1980 function.""" + si_track = xr.Dataset({ + "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), + "vgrad": ("time", [35, 40]), + "rho_air": ("time", [1.15, 1.15]) + }) + _B_holland_1980(si_track) + np.testing.assert_array_almost_equal(si_track["hol_b"], [2.5, 1.667213]) + + si_track = xr.Dataset({ + "pdelta": ("time", MBAR_TO_PA * np.array([4.74, 15, 30, 40])), + "vmax": ("time", [np.nan, 22.5, 25.4, 42.5]), + "rho_air": ("time", [1.2, 1.2, 1.2, 1.2]) + }) + _B_holland_1980(si_track, gradient_to_surface_winds=0.9) + np.testing.assert_allclose(si_track["hol_b"], [np.nan, 1.101, 0.810, 1.473], atol=1e-3) + + def test_bs_holland_2008_pass(self): + """Test _bs_holland_2008 function. Compare to MATLAB reference.""" + si_track = xr.Dataset({ + "tstep": ("time", H_TO_S * np.array([1.0, 1.0, 1.0])), + "lat": ("time", [12.299999504631234, 12.299999504631343, 12.299999279463769]), + "pdelta": ("time", MBAR_TO_PA * np.array([4.74, 4.73, 4.73])), + "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 1005.2633, 1005.2682])), + "vtrans_norm": ("time", [np.nan, 5.241999541820597, 5.123882725120426]), + }) + _bs_holland_2008(si_track) + np.testing.assert_allclose( + si_track["hol_b"], [np.nan, 1.27, 1.27], atol=1e-2 + ) + + def test_v_max_s_holland_2008_pass(self): + """Test _v_max_s_holland_2008 function.""" + # Numbers analogous to test_B_holland_1980_pass + si_track = xr.Dataset({ + "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), + "hol_b": ("time", [2.5, 1.67]), + "rho_air": ("time", [1.15, 1.15]), + }) + _v_max_s_holland_2008(si_track) + np.testing.assert_array_almost_equal(si_track["vmax"], [34.635341, 40.033421]) + + def test_holland_2010_pass(self): + """Test Holland et al. 2010 wind field model.""" + # The parameter "x" is designed to be exactly 0.5 inside the radius of max wind (RMW) and + # to increase or decrease linearly outside of it in radial direction. + # + # An increase (decrease) of "x" outside of the RMW is for cases where the max wind is very + # high (low), but the RMW is still comparably large (small). This means, wind speeds need + # to decay very sharply (only moderately) outside of the RMW to reach the low prescribed + # peripheral wind speeds. + # + # The "hol_b" parameter tunes the meaning of a "comparably" large or small RMW. + si_track = xr.Dataset({ + # four test cases: + # - low vmax, moderate RMW: x decreases moderately + # - large hol_b: x decreases sharply + # - very low vmax: x decreases so much, it needs to be clipped at 0 + # - large vmax, large RMW: x increases + "rad": ("time", KM_TO_M * np.array([75, 75, 75, 90])), + "vmax": ("time", [35.0, 35.0, 16.0, 90.0]), + "hol_b": ("time", [1.75, 2.5, 1.9, 1.6]), + }) + d_centr = KM_TO_M * np.array([ + # first column is for locations within the storm eye + # second column is for locations at or close to the radius of max wind + # third column is for locations outside the storm eye + # fourth column is for locations exactly at the peripheral radius + # fifth column is for locations outside the peripheral radius + [0., 75, 220, 300, 490], + [30, 74, 170, 300, 501], + [21, 76, 230, 300, 431], + [32, 91, 270, 300, 452], + ], dtype=float) + close_centr = np.array([ + # note that we set one of these to "False" for testing + [True, True, True, True, True], + [True, True, True, True, False], + [True, True, True, True, True], + [True, True, True, True, True], + ], dtype=bool) + hol_x = _x_holland_2010(si_track, d_centr, close_centr) + np.testing.assert_array_almost_equal(hol_x, [ + [0.5, 0.500000, 0.485077, 0.476844, 0.457291], + [0.5, 0.500000, 0.410997, 0.400000, 0.000000], + [0.5, 0.497620, 0.400000, 0.400000, 0.400000], + [0.5, 0.505022, 1.403952, 1.554611, 2.317948], + ]) + + v_ang_norm = _stat_holland_2010(si_track, d_centr, close_centr, hol_x) + np.testing.assert_allclose(v_ang_norm, [ + # first column: converge to 0 when approaching storm eye + # second column: vmax at RMW + # fourth column: peripheral speed (17 m/s) at peripheral radius (unless x is clipped!) + [ 0.000000, 35.000000, 21.181497, 17.000000, 12.1034610], + [ 1.296480, 34.990037, 21.593755, 12.891313, 0.0000000], + [ 0.321952, 15.997500, 9.712006, 8.087240, 6.2289690], + [24.823469, 89.992938, 24.381965, 17.000000, 1.9292020], + ], atol=1e-6) + + def test_stat_holland_1980(self): + """Test _stat_holland_1980 function. Compare to MATLAB reference.""" + d_centr = KM_TO_M * np.array([ + [299.4501244109841, 291.0737897183741, 292.5441003235722, 40.665454622610511], + [293.6067129546862, 1000.0, 298.2652319413182, 70.0], + ]) + si_track = xr.Dataset({ + "rad": ("time", KM_TO_M * np.array([40.665454622610511, 75.547902916671745])), + "hol_b": ("time", [1.486076257880692, 1.265551666104679]), + "pdelta": ("time", MBAR_TO_PA * np.array([39.12, 4.73])), + "lat": ("time", [-14.089110370469488, 12.299999279463769]), + "cp": ("time", [3.54921922e-05, 3.10598285e-05]), + "rho_air": ("time", [1.15, 1.15]), + }) + mask = np.array([[True, True, True, True], [True, False, True, True]], dtype=bool) + v_ang_norm = _stat_holland_1980(si_track, d_centr, mask) + np.testing.assert_allclose( + v_ang_norm, [[11.28, 11.68, 11.61, 42.41], [5.38, 0, 5.28, 12.76]], atol=1e-2, + ) + + # without Coriolis force, values are higher, esp. far away from the center: + v_ang_norm = _stat_holland_1980(si_track, d_centr, mask, cyclostrophic=True) + np.testing.assert_allclose( + v_ang_norm, [[15.72, 16.04, 15.98, 43.13], [8.84, 0, 8.76, 13.81]], atol=1e-2, + ) + + d_centr = np.array([[], []]) + mask = np.ones_like(d_centr, dtype=bool) + v_ang_norm = _stat_holland_1980(si_track, d_centr, mask) + np.testing.assert_array_equal(v_ang_norm, np.array([[], []])) + + def test_er_2011_pass(self): + """Test Emanuel and Rotunno 2011 wind field model.""" + # test at centroids within and outside of radius of max wind + d_centr = KM_TO_M * np.array([[35, 70, 75, 220], [30, 150, 1000, 300]], dtype=float) + si_track = xr.Dataset({ + "rad": ("time", KM_TO_M * np.array([75.0, 40.0])), + "vmax": ("time", [35.0, 40.0]), + "lat": ("time", [20.0, 27.0]), + "cp": ("time", [4.98665369e-05, 6.61918149e-05]), + }) + mask = np.array([[True, True, True, True], [True, False, True, True]], dtype=bool) + v_ang_norm = _stat_er_2011(si_track, d_centr, mask) + np.testing.assert_array_almost_equal(v_ang_norm, + [[28.258025, 36.782418, 36.869995, 22.521237], + [39.670883, 0, 3.300626, 10.827206]]) + + def test_vtrans_pass(self): + """Test _vtrans function. Compare to MATLAB reference.""" + tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK) + tc_track.equal_timestep() + track_ds = tc_track.data[0] + + si_track = tctrack_to_si(track_ds) + _vtrans(si_track) + + to_kn = (1.0 * ureg.meter / ureg.second).to(ureg.knot).magnitude + + self.assertEqual(si_track["vtrans_norm"].size, track_ds["time"].size) + self.assertEqual(si_track["vtrans_norm"].values[0], 0) + self.assertAlmostEqual(si_track["vtrans_norm"].values[1] * to_kn, 10.191466246) + + def testtctrack_to_si(self): + """ Test tctrack_to_si should create the same vmax output independent of the input unit """ + tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK_SHORT).data[0] + + tc_track_kmph = tc_track.copy(deep=True) + tc_track_kmph['max_sustained_wind'] *= ( + (1.0 * ureg.knot).to(ureg.km / ureg.hour).magnitude + ) + tc_track_kmph.attrs['max_sustained_wind_unit'] = 'km/h' + + si_track = tctrack_to_si(tc_track) + si_track_from_kmph = tctrack_to_si(tc_track_kmph) + + np.testing.assert_array_almost_equal(si_track["vmax"], si_track_from_kmph["vmax"]) + + tc_track.attrs['max_sustained_wind_unit'] = 'elbows/fortnight' + with self.assertRaises(ValueError): + tctrack_to_si(tc_track) + + +if __name__ == "__main__": + TESTS = unittest.TestLoader().loadTestsFromTestCase(TestWindfieldHelpers) + unittest.TextTestRunner(verbosity=2).run(TESTS) diff --git a/climada/hazard/trop_cyclone/__init__.py b/climada/hazard/trop_cyclone/__init__.py new file mode 100644 index 0000000000..452bf43646 --- /dev/null +++ b/climada/hazard/trop_cyclone/__init__.py @@ -0,0 +1,5 @@ +from climada.hazard.trop_cyclone.trop_cyclone import * +from climada.hazard.trop_cyclone.trop_cyclone_windfields import compute_windfields_sparse, compute_angular_windspeeds, tctrack_to_si, \ + get_close_centroids, KN_TO_MS, KM_TO_M, KM_TO_M, H_TO_S, NM_TO_KM, KMH_TO_MS, MBAR_TO_PA, \ + DEF_MAX_DIST_EYE_KM, DEF_INTENSITY_THRES, DEF_MAX_MEMORY_GB, MODEL_VANG, DEF_RHO_AIR, DEF_GRADIENT_TO_SURFACE_WINDS, \ + T_ICE_K, V_ANG_EARTH diff --git a/climada/hazard/trop_cyclone/trop_cyclone.py b/climada/hazard/trop_cyclone/trop_cyclone.py new file mode 100644 index 0000000000..6fc8347e97 --- /dev/null +++ b/climada/hazard/trop_cyclone/trop_cyclone.py @@ -0,0 +1,696 @@ +""" +This file is part of CLIMADA. + +Copyright (C) 2017 ETH Zurich, CLIMADA contributors listed in AUTHORS. + +CLIMADA is free software: you can redistribute it and/or modify it under the +terms of the GNU General Public License as published by the Free +Software Foundation, version 3. + +CLIMADA is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A +PARTICULAR PURPOSE. See the GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with CLIMADA. If not, see . + +--- + +Define TC wind hazard (TropCyclone class). +""" + +__all__ = ['TropCyclone'] + +import copy +import datetime as dt +import itertools +import logging +import time +from typing import Optional, Tuple, List + +import numpy as np +from scipy import sparse +import matplotlib.animation as animation +from tqdm import tqdm +import pathos.pools +import xarray as xr + +from climada.hazard.base import Hazard +from climada.hazard.tc_tracks import TCTracks +from climada.hazard.tc_clim_change import get_knutson_criterion, calc_scale_knutson +from climada.hazard.centroids.centr import Centroids +import climada.util.constants as u_const +import climada.util.coordinates as u_coord +import climada.util.plot as u_plot + +from .trop_cyclone_windfields import DEF_MAX_DIST_EYE_KM, DEF_INTENSITY_THRES, \ + DEF_MAX_MEMORY_GB, compute_windfields_sparse + +LOGGER = logging.getLogger(__name__) + +HAZ_TYPE = 'TC' +"""Hazard type acronym for Tropical Cyclone""" + + +class TropCyclone(Hazard): + """ + Contains tropical cyclone events. + + Attributes + ---------- + category : np.ndarray of ints + for every event, the TC category using the Saffir-Simpson scale: + + * -1 tropical depression + * 0 tropical storm + * 1 Hurrican category 1 + * 2 Hurrican category 2 + * 3 Hurrican category 3 + * 4 Hurrican category 4 + * 5 Hurrican category 5 + basin : list of str + Basin where every event starts: + + * 'NA' North Atlantic + * 'EP' Eastern North Pacific + * 'WP' Western North Pacific + * 'NI' North Indian + * 'SI' South Indian + * 'SP' Southern Pacific + * 'SA' South Atlantic + windfields : list of csr_matrix + For each event, the full velocity vectors at each centroid and track position in a sparse + matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full ndarray of + shape (npositions, ncentroids, 2). + """ + intensity_thres = DEF_INTENSITY_THRES + """intensity threshold for storage in m/s""" + + vars_opt = Hazard.vars_opt.union({'category'}) + """Name of the variables that are not needed to compute the impact.""" + + def __init__( + self, + category: Optional[np.ndarray] = None, + basin: Optional[List] = None, + windfields: Optional[List[sparse.csr_matrix]] = None, + **kwargs, + ): + """Initialize values. + + Parameters + ---------- + category : np.ndarray of int, optional + For every event, the TC category using the Saffir-Simpson scale: + -1 tropical depression + 0 tropical storm + 1 Hurrican category 1 + 2 Hurrican category 2 + 3 Hurrican category 3 + 4 Hurrican category 4 + 5 Hurrican category 5 + basin : list of str, optional + Basin where every event starts: + 'NA' North Atlantic + 'EP' Eastern North Pacific + 'WP' Western North Pacific + 'NI' North Indian + 'SI' South Indian + 'SP' Southern Pacific + 'SA' South Atlantic + windfields : list of csr_matrix, optional + For each event, the full velocity vectors at each centroid and track position in a + sparse matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full + ndarray of shape (npositions, ncentroids, 2). + **kwargs : Hazard properties, optional + All other keyword arguments are passed to the Hazard constructor. + """ + kwargs.setdefault('haz_type', HAZ_TYPE) + Hazard.__init__(self, **kwargs) + self.category = category if category is not None else np.array([], int) + self.basin = basin if basin is not None else [] + self.windfields = windfields if windfields is not None else [] + + def set_from_tracks(self, *args, **kwargs): + """This function is deprecated, use TropCyclone.from_tracks instead.""" + LOGGER.warning("The use of TropCyclone.set_from_tracks is deprecated." + "Use TropCyclone.from_tracks instead.") + if "intensity_thres" not in kwargs: + # some users modify the threshold attribute before calling `set_from_tracks` + kwargs["intensity_thres"] = self.intensity_thres + if self.pool is not None and 'pool' not in kwargs: + kwargs['pool'] = self.pool + self.__dict__ = TropCyclone.from_tracks(*args, **kwargs).__dict__ + + @classmethod + def from_tracks( + cls, + tracks: TCTracks, + centroids: Centroids, + pool: Optional[pathos.pools.ProcessPool] = None, + model: str = 'H08', + model_kwargs: Optional[dict] = None, + ignore_distance_to_coast: bool = False, + store_windfields: bool = False, + metric: str = "equirect", + intensity_thres: float = DEF_INTENSITY_THRES, + max_latitude: float = 61, + max_dist_inland_km: float = 1000, + max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, + max_memory_gb: float = DEF_MAX_MEMORY_GB, + ): + """ + Create new TropCyclone instance that contains windfields from the specified tracks. + + This function sets the ``intensity`` attribute to contain, for each centroid, + the maximum wind speed (1-minute sustained winds at 10 meters above ground) + experienced over the whole period of each TC event in m/s. The wind speed is set + to 0 if it doesn't exceed the threshold ``intensity_thres``. + + The ``category`` attribute is set to the value of the ``category``-attribute + of each of the given track data sets. + + The ``basin`` attribute is set to the genesis basin for each event, which + is the first value of the ``basin``-variable in each of the given track data sets. + + Optionally, the time dependent, vectorial winds can be stored using the + ``store_windfields`` function parameter (see below). + + Parameters + ---------- + tracks : climada.hazard.TCTracks + Tracks of storm events. + centroids : Centroids, optional + Centroids where to model TC. Default: global centroids at 360 arc-seconds + resolution. + pool : pathos.pool, optional + Pool that will be used for parallel computation of wind fields. Default: + None + model : str, optional + Parametric wind field model to use. Default: "H08". + + * ``"H1980"`` (the prominent Holland 1980 model) from the paper: + Holland, G.J. (1980): An Analytic Model of the Wind and Pressure + Profiles in Hurricanes. Monthly Weather Review 108(8): 1212–1218. + ``https://doi.org/10.1175/1520-0493(1980)108<1212:AAMOTW>2.0.CO;2`` + * ``"H08"`` (Holland 1980 with b-value from Holland 2008) from the paper: + Holland, G. (2008). A revised hurricane pressure-wind model. Monthly + Weather Review, 136(9), 3432–3445. + https://doi.org/10.1175/2008MWR2395.1 + * ``"H10"`` (Holland et al. 2010) from the paper: + Holland et al. (2010): A Revised Model for Radial Profiles of + Hurricane Winds. Monthly Weather Review 138(12): 4393–4401. + https://doi.org/10.1175/2010MWR3317.1 + * ``"ER11"`` (Emanuel and Rotunno 2011) from the paper: + Emanuel, K., Rotunno, R. (2011): Self-Stratification of Tropical + Cyclone Outflow. Part I: Implications for Storm Structure. Journal + of the Atmospheric Sciences 68(10): 2236–2249. + https://dx.doi.org/10.1175/JAS-D-10-05024.1 + model_kwargs : dict, optional + If given, forward these kwargs to the selected wind model. None of the + parameters is currently supported by the ER11 model. Default: None. + The Holland models support the following parameters, in alphabetical order: + + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use. In H1980, the wind + profile is computed on the gradient level, and wind speeds are converted + to the surface level using this factor. In H08 and H10, the wind profile + is computed on the surface level, but the clipping interval of the + B-value depends on this factor. Default: 0.9 + rho_air_const : float or None, optional + The constant value for air density (in kg/m³) to assume in the formulas + from Holland 1980. By default, the constant value suggested in Holland + 1980 is used. If set to None, the air density is computed from pressure + following equation (9) in Holland et al. 2010. Default: 1.15 + vmax_from_cen : boolean, optional + Only used in H10. If True, replace the recorded value of vmax along the + track by an estimate from pressure, following equation (8) in Holland et + al. 2010. Default: True + vmax_in_brackets : bool, optional + Only used in H10. Specifies which of the two formulas in equation (6) of + Holland et al. 2010 to use. If False, the formula with vmax outside of + the brackets is used. Note that, a side-effect of the formula with vmax + inside of the brackets is that the wind speed maximum is attained a bit + farther away from the center than according to the recorded radius of + maximum winds (RMW). Default: False + + ignore_distance_to_coast : boolean, optional + If True, centroids far from coast are not ignored. + If False, the centroids' distances to the coast are calculated with the + `Centroids.get_dist_coast()` method (unless there is "dist_coast" column in + the centroids' GeoDataFrame) and centroids far from coast are ignored. + Default: False. + store_windfields : boolean, optional + If True, the Hazard object gets a list ``windfields`` of sparse matrices. + For each track, the full velocity vectors at each centroid and track + position are stored in a sparse matrix of shape (npositions, + ncentroids * 2) that can be reshaped to a full ndarray of shape (npositions, + ncentroids, 2). Default: False. + metric : str, optional + Specify an approximation method to use for earth distances: + + * "equirect": Distance according to sinusoidal projection. Fast, but + inaccurate for large distances and high latitudes. + * "geosphere": Exact spherical distance. Much more accurate at all + distances, but slow. + + Default: "equirect". + intensity_thres : float, optional + Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 + max_latitude : float, optional + No wind speed calculation is done for centroids with latitude larger than + this parameter. Default: 61 + max_dist_inland_km : float, optional + No wind speed calculation is done for centroids with a distance (in km) to + the coast larger than this parameter. Default: 1000 + max_dist_eye_km : float, optional + No wind speed calculation is done for centroids with a distance (in km) to + the TC center ("eye") larger than this parameter. Default: 300 + max_memory_gb : float, optional + To avoid memory issues, the computation is done for chunks of the track + sequentially. The chunk size is determined depending on the available memory + (in GB). Note that this limit applies to each thread separately if a + ``pool`` is used. Default: 8 + + Raises + ------ + ValueError + + Returns + ------- + TropCyclone + """ + num_tracks = tracks.size + + if ignore_distance_to_coast: + # Select centroids with lat <= max_latitude + [idx_centr_filter] = (np.abs(centroids.lat) <= max_latitude).nonzero() + else: + # Select centroids which are inside max_dist_inland_km and lat <= max_latitude + if 'dist_coast' not in centroids.gdf.columns: + dist_coast = centroids.get_dist_coast() + else: + dist_coast = centroids.gdf.dist_coast.values + [idx_centr_filter] = ( + (dist_coast <= max_dist_inland_km * 1000) + & (np.abs(centroids.lat) <= max_latitude) + ).nonzero() + + # Filter early with a larger threshold, but inaccurate (lat/lon) distances. + # Later, there will be another filtering step with more accurate distances in km. + max_dist_eye_deg = max_dist_eye_km / ( + u_const.ONE_LAT_KM * np.cos(np.radians(max_latitude)) + ) + + # Restrict to coastal centroids within reach of any of the tracks + t_lon_min, t_lat_min, t_lon_max, t_lat_max = tracks.get_bounds(deg_buffer=max_dist_eye_deg) + t_mid_lon = 0.5 * (t_lon_min + t_lon_max) + filtered_centroids = centroids.coord[idx_centr_filter] + u_coord.lon_normalize(filtered_centroids[:, 1], center=t_mid_lon) + idx_centr_filter = idx_centr_filter[ + (t_lon_min <= filtered_centroids[:, 1]) + & (filtered_centroids[:, 1] <= t_lon_max) + & (t_lat_min <= filtered_centroids[:, 0]) + & (filtered_centroids[:, 0] <= t_lat_max) + ] + + # prepare keyword arguments to pass to `from_single_track` + kwargs_from_single_track = dict( + centroids=centroids, + idx_centr_filter=idx_centr_filter, + model=model, + model_kwargs=model_kwargs, + store_windfields=store_windfields, + metric=metric, + intensity_thres=intensity_thres, + max_dist_eye_km=max_dist_eye_km, + max_memory_gb=max_memory_gb, + ) + + LOGGER.info( + 'Mapping %d tracks to %d coastal centroids.', num_tracks, idx_centr_filter.size, + ) + if pool: + chunksize = max(min(num_tracks // pool.ncpus, 1000), 1) + kwargs_repeated = [ + itertools.repeat(val, num_tracks) + for val in kwargs_from_single_track.values() + ] + tc_haz_list = pool.map( + cls.from_single_track, + tracks.data, + *kwargs_repeated, + chunksize=chunksize, + ) + else: + last_perc = 0 + tc_haz_list = [] + for track in tracks.data: + perc = 100 * len(tc_haz_list) / len(tracks.data) + if perc - last_perc >= 10: + LOGGER.info("Progress: %d%%", perc) + last_perc = perc + tc_haz_list.append( + cls.from_single_track(track, **kwargs_from_single_track) + ) + if last_perc < 100: + LOGGER.info("Progress: 100%") + + LOGGER.debug('Concatenate events.') + haz = cls.concat(tc_haz_list) + haz.pool = pool + haz.intensity_thres = intensity_thres + LOGGER.debug('Compute frequency.') + haz.frequency_from_tracks(tracks.data) + return haz + + def apply_climate_scenario_knu( + self, + ref_year: int = 2050, + rcp_scenario: int = 45 + ): + """ + From current TC hazard instance, return new hazard set with + future events for a given RCP scenario and year based on the + parametrized values derived from Table 3 in Knutson et al 2015. + https://doi.org/10.1175/JCLI-D-15-0129.1 . The scaling for different + years and RCP scenarios is obtained by linear interpolation. + + Note: The parametrized values are derived from the overall changes + in statistical ensemble of tracks. Hence, this method should only be + applied to sufficiently large tropical cyclone event sets that + approximate the reference years 1981 - 2008 used in Knutson et. al. + + The frequency and intensity changes are applied independently from + one another. The mean intensity factors can thus slightly deviate + from the Knutson value (deviation was found to be less than 1% + for default IBTrACS event sets 1980-2020 for each basin). + + Parameters + ---------- + ref_year : int + year between 2000 ad 2100. Default: 2050 + rcp_scenario : int + 26 for RCP 2.6, 45 for RCP 4.5, 60 for RCP 6.0 and 85 for RCP 8.5. + The default is 45. + + Returns + ------- + haz_cc : climada.hazard.TropCyclone + Tropical cyclone with frequencies and intensity scaled according + to the Knutson criterion for the given year and RCP. Returns + a new instance of climada.hazard.TropCyclone, self is not + modified. + """ + chg_int_freq = get_knutson_criterion() + scale_rcp_year = calc_scale_knutson(ref_year, rcp_scenario) + haz_cc = self._apply_knutson_criterion(chg_int_freq, scale_rcp_year) + return haz_cc + + def set_climate_scenario_knu(self, *args, **kwargs): + """This function is deprecated, use TropCyclone.apply_climate_scenario_knu instead.""" + LOGGER.warning("The use of TropCyclone.set_climate_scenario_knu is deprecated." + "Use TropCyclone.apply_climate_scenario_knu instead.") + return self.apply_climate_scenario_knu(*args, **kwargs) + + @classmethod + def video_intensity( + cls, + track_name: str, + tracks: TCTracks, + centroids: Centroids, + file_name: Optional[str] = None, + writer: animation = animation.PillowWriter(bitrate=500), + figsize: Tuple[float, float] = (9, 13), + adapt_fontsize: bool = True, + **kwargs + ): + """ + Generate video of TC wind fields node by node and returns its + corresponding TropCyclone instances and track pieces. + + Parameters + ---------- + track_name : str + name of the track contained in tracks to record + tracks : climada.hazard.TCTracks + tropical cyclone tracks + centroids : climada.hazard.Centroids + centroids where wind fields are mapped + file_name : str, optional + file name to save video (including full path and file extension) + writer : matplotlib.animation.*, optional + video writer. Default is pillow with bitrate=500 + figsize : tuple, optional + figure size for plt.subplots + adapt_fontsize : bool, optional + If set to true, the size of the fonts will be adapted to the size of the figure. + Otherwise the default matplotlib font size is used. Default is True. + kwargs : optional + arguments for pcolormesh matplotlib function used in event plots + + Returns + ------- + tc_list, tc_coord : list(TropCyclone), list(np.ndarray) + + Raises + ------ + ValueError + + """ + # initialization + track = tracks.get_track(track_name) + if not track: + raise ValueError(f'{track_name} not found in track data.') + idx_plt = np.argwhere( + (track.lon.values < centroids.total_bounds[2] + 1) + & (centroids.total_bounds[0] - 1 < track.lon.values) + & (track.lat.values < centroids.total_bounds[3] + 1) + & (centroids.total_bounds[1] - 1 < track.lat.values) + ).reshape(-1) + + tc_list = [] + tr_coord = {'lat': [], 'lon': []} + for node in range(idx_plt.size - 2): + tr_piece = track.sel( + time=slice(track.time.values[idx_plt[node]], + track.time.values[idx_plt[node + 2]])) + tr_piece.attrs['n_nodes'] = 2 # plot only one node + tr_sel = TCTracks() + tr_sel.append(tr_piece) + tr_coord['lat'].append(tr_sel.data[0].lat.values[:-1]) + tr_coord['lon'].append(tr_sel.data[0].lon.values[:-1]) + + tc_tmp = cls.from_tracks(tr_sel, centroids=centroids) + tc_tmp.event_name = [ + track.name + ' ' + time.strftime( + "%d %h %Y %H:%M", + time.gmtime(tr_sel.data[0].time[1].values.astype(int) + / 1000000000) + ) + ] + tc_list.append(tc_tmp) + + if 'cmap' not in kwargs: + kwargs['cmap'] = 'Greys' + if 'vmin' not in kwargs: + kwargs['vmin'] = np.array([tc_.intensity.min() for tc_ in tc_list]).min() + if 'vmax' not in kwargs: + kwargs['vmax'] = np.array([tc_.intensity.max() for tc_ in tc_list]).max() + + def run(node): + tc_list[node].plot_intensity(1, axis=axis, **kwargs) + axis.plot(tr_coord['lon'][node], tr_coord['lat'][node], 'k') + axis.set_title(tc_list[node].event_name[0]) + pbar.update() + + if file_name: + LOGGER.info('Generating video %s', file_name) + fig, axis, _fontsize = u_plot.make_map(figsize=figsize, adapt_fontsize=adapt_fontsize) + pbar = tqdm(total=idx_plt.size - 2) + ani = animation.FuncAnimation(fig, run, frames=idx_plt.size - 2, + interval=500, blit=False) + fig.tight_layout() + ani.save(file_name, writer=writer) + pbar.close() + return tc_list, tr_coord + + def frequency_from_tracks(self, tracks: List): + """ + Set hazard frequency from tracks data. + + Parameters + ---------- + tracks : list of xarray.Dataset + """ + if not tracks: + return + year_max = np.amax([t.time.dt.year.values.max() for t in tracks]) + year_min = np.amin([t.time.dt.year.values.min() for t in tracks]) + year_delta = year_max - year_min + 1 + num_orig = np.count_nonzero(self.orig) + ens_size = (self.event_id.size / num_orig) if num_orig > 0 else 1 + self.frequency = np.ones(self.event_id.size) / (year_delta * ens_size) + + @classmethod + def from_single_track( + cls, + track: xr.Dataset, + centroids: Centroids, + idx_centr_filter: np.ndarray, + model: str = 'H08', + model_kwargs: Optional[dict] = None, + store_windfields: bool = False, + metric: str = "equirect", + intensity_thres: float = DEF_INTENSITY_THRES, + max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, + max_memory_gb: float = DEF_MAX_MEMORY_GB, + ): + """ + Generate windfield hazard from a single track dataset + + Parameters + ---------- + track : xr.Dataset + Single tropical cyclone track. + centroids : Centroids + Centroids instance. + idx_centr_filter : np.ndarray + Indices of centroids to restrict to (e.g. sufficiently close to coast). + model : str, optional + Parametric wind field model, one of "H1980" (the prominent Holland 1980 model), + "H08" (Holland 1980 with b-value from Holland 2008), "H10" (Holland et al. 2010), or + "ER11" (Emanuel and Rotunno 2011). + Default: "H08". + model_kwargs: dict, optional + If given, forward these kwargs to the selected model. Default: None + store_windfields : boolean, optional + If True, store windfields. Default: False. + metric : str, optional + Specify an approximation method to use for earth distances: "equirect" (faster) or + "geosphere" (more accurate). See ``dist_approx`` function in + ``climada.util.coordinates``. + Default: "equirect". + intensity_thres : float, optional + Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 + max_dist_eye_km : float, optional + No wind speed calculation is done for centroids with a distance (in km) to the TC + center ("eye") larger than this parameter. Default: 300 + max_memory_gb : float, optional + To avoid memory issues, the computation is done for chunks of the track sequentially. + The chunk size is determined depending on the available memory (in GB). Default: 8 + + Raises + ------ + ValueError, KeyError + + Returns + ------- + haz : TropCyclone + """ + intensity_sparse, windfields_sparse = compute_windfields_sparse( + track=track, + centroids=centroids, + idx_centr_filter=idx_centr_filter, + model=model, + model_kwargs=model_kwargs, + store_windfields=store_windfields, + metric=metric, + intensity_thres=intensity_thres, + max_dist_eye_km=max_dist_eye_km, + max_memory_gb=max_memory_gb, + ) + + new_haz = cls(haz_type=HAZ_TYPE) + new_haz.intensity_thres = intensity_thres + new_haz.intensity = intensity_sparse + if store_windfields: + new_haz.windfields = [windfields_sparse] + new_haz.units = 'm/s' + new_haz.centroids = centroids + new_haz.event_id = np.array([1]) + new_haz.frequency = np.array([1]) + new_haz.event_name = [track.sid] + new_haz.fraction = sparse.csr_matrix(new_haz.intensity.shape) + # store first day of track as date + new_haz.date = np.array([ + dt.datetime(track.time.dt.year.values[0], + track.time.dt.month.values[0], + track.time.dt.day.values[0]).toordinal() + ]) + new_haz.orig = np.array([track.orig_event_flag]) + new_haz.category = np.array([track.category]) + # users that pickle TCTracks objects might still have data with the legacy basin attribute, + # so we have to deal with it here + new_haz.basin = [track.basin if isinstance(track.basin, str) + else str(track.basin.values[0])] + return new_haz + + def _apply_knutson_criterion( + self, + chg_int_freq: List, + scaling_rcp_year: float + ): + """ + Apply changes to intensities and cumulative frequencies. + + Parameters + ---------- + chg_int_freq : list(dict)) + list of criteria from climada.hazard.tc_clim_change + scaling_rcp_year : float + scale parameter because of chosen year and RCP + + Returns + ------- + tc_cc : climada.hazard.TropCyclone + Tropical cyclone with frequency and intensity scaled inspired by + the Knutson criterion. Returns a new instance of TropCyclone. + """ + + tc_cc = copy.deepcopy(self) + + # Criterion per basin + for basin in np.unique(tc_cc.basin): + + bas_sel = np.array(tc_cc.basin) == basin + + # Apply intensity change + inten_chg = [chg + for chg in chg_int_freq + if (chg['variable'] == 'intensity' and + chg['basin'] == basin) + ] + for chg in inten_chg: + sel_cat_chg = np.isin(tc_cc.category, chg['category']) & bas_sel + inten_scaling = 1 + (chg['change'] - 1) * scaling_rcp_year + tc_cc.intensity = sparse.diags( + np.where(sel_cat_chg, inten_scaling, 1) + ).dot(tc_cc.intensity) + + # Apply frequency change + freq_chg = [chg + for chg in chg_int_freq + if (chg['variable'] == 'frequency' and + chg['basin'] == basin) + ] + freq_chg.sort(reverse=False, key=lambda x: len(x['category'])) + + # Scale frequencies by category + cat_larger_list = [] + for chg in freq_chg: + cat_chg_list = [cat + for cat in chg['category'] + if cat not in cat_larger_list + ] + sel_cat_chg = np.isin(tc_cc.category, cat_chg_list) & bas_sel + if sel_cat_chg.any(): + freq_scaling = 1 + (chg['change'] - 1) * scaling_rcp_year + tc_cc.frequency[sel_cat_chg] *= freq_scaling + cat_larger_list += cat_chg_list + + if (tc_cc.frequency < 0).any(): + raise ValueError("The application of the given climate scenario" + "resulted in at least one negative frequency.") + + return tc_cc diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone/trop_cyclone_windfields.py similarity index 67% rename from climada/hazard/trop_cyclone.py rename to climada/hazard/trop_cyclone/trop_cyclone_windfields.py index b1bf7a72c6..e82c0b11e5 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone/trop_cyclone_windfields.py @@ -16,38 +16,28 @@ --- -Define TC wind hazard (TropCyclone class). +Compute Tropical Cyclone windfields (see compute_windfields_sparse function). """ -__all__ = ['TropCyclone'] - -import copy -import datetime as dt -import itertools import logging -import time -from typing import Optional, Tuple, List, Union +from typing import Optional, Union, Tuple import numpy as np -from scipy import sparse -import matplotlib.animation as animation -from tqdm import tqdm -import pathos.pools import xarray as xr +from scipy import sparse -from climada.hazard.base import Hazard -from climada.hazard.tc_tracks import TCTracks, estimate_rmw -from climada.hazard.tc_clim_change import get_knutson_criterion, calc_scale_knutson -from climada.hazard.centroids.centr import Centroids -from climada.util import ureg -import climada.util.constants as u_const -import climada.util.coordinates as u_coord -import climada.util.plot as u_plot +from climada.hazard import Centroids +from climada.hazard.tc_tracks import estimate_rmw +from climada.util import ureg, coordinates as u_coord, constants as u_const LOGGER = logging.getLogger(__name__) -HAZ_TYPE = 'TC' -"""Hazard type acronym for Tropical Cyclone""" +NM_TO_KM = (1.0 * ureg.nautical_mile).to(ureg.kilometer).magnitude +KMH_TO_MS = (1.0 * ureg.km / ureg.hour).to(ureg.meter / ureg.second).magnitude +KM_TO_M = (1.0 * ureg.kilometer).to(ureg.meter).magnitude +H_TO_S = (1.0 * ureg.hours).to(ureg.seconds).magnitude +MBAR_TO_PA = (1.0 * ureg.millibar).to(ureg.pascal).magnitude +"""Unit conversion factors for JIT functions that can't use ureg""" DEF_MAX_DIST_EYE_KM = 300 """Default value for the maximum distance (in km) of a centroid to the TC center at which wind @@ -76,1407 +66,744 @@ other regions and levels, values of 0.8 or even 0.75 might be justified. """ -KMH_TO_MS = (1.0 * ureg.km / ureg.hour).to(ureg.meter / ureg.second).magnitude -KN_TO_MS = (1.0 * ureg.knot).to(ureg.meter / ureg.second).magnitude -NM_TO_KM = (1.0 * ureg.nautical_mile).to(ureg.kilometer).magnitude -KM_TO_M = (1.0 * ureg.kilometer).to(ureg.meter).magnitude -H_TO_S = (1.0 * ureg.hours).to(ureg.seconds).magnitude -MBAR_TO_PA = (1.0 * ureg.millibar).to(ureg.pascal).magnitude -"""Unit conversion factors for JIT functions that can't use ureg""" - T_ICE_K = 273.16 """Freezing temperatur of water (in K), for conversion between K and °C""" V_ANG_EARTH = 7.29e-5 """Earth angular velocity (in radians per second)""" -class TropCyclone(Hazard): - """ - Contains tropical cyclone events. +def _vgrad(si_track, gradient_to_surface_winds): + """Gradient wind speeds (in m/s) without translational influence at each track node - Attributes + The track dataset is modified in place, with the "vgrad" data variable added. + + Parameters ---------- - category : np.ndarray of ints - for every event, the TC category using the Saffir-Simpson scale: - - * -1 tropical depression - * 0 tropical storm - * 1 Hurrican category 1 - * 2 Hurrican category 2 - * 3 Hurrican category 3 - * 4 Hurrican category 4 - * 5 Hurrican category 5 - basin : list of str - Basin where every event starts: - - * 'NA' North Atlantic - * 'EP' Eastern North Pacific - * 'WP' Western North Pacific - * 'NI' North Indian - * 'SI' South Indian - * 'SP' Southern Pacific - * 'SA' South Atlantic - windfields : list of csr_matrix - For each event, the full velocity vectors at each centroid and track position in a sparse - matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full ndarray of - shape (npositions, ncentroids, 2). + si_track : xr.Dataset + Track information as returned by `tctrack_to_si`. The data variables used by this function + are "vmax" and "vtrans_norm". The result is stored in place as new data variable "vgrad". + gradient_to_surface_winds : float + The gradient-to-surface wind reduction factor to use. """ - intensity_thres = DEF_INTENSITY_THRES - """intensity threshold for storage in m/s""" - - vars_opt = Hazard.vars_opt.union({'category'}) - """Name of the variables that are not needed to compute the impact.""" - - def __init__( - self, - category: Optional[np.ndarray] = None, - basin: Optional[List] = None, - windfields: Optional[List[sparse.csr_matrix]] = None, - **kwargs, - ): - """Initialize values. - - Parameters - ---------- - category : np.ndarray of int, optional - For every event, the TC category using the Saffir-Simpson scale: - -1 tropical depression - 0 tropical storm - 1 Hurrican category 1 - 2 Hurrican category 2 - 3 Hurrican category 3 - 4 Hurrican category 4 - 5 Hurrican category 5 - basin : list of str, optional - Basin where every event starts: - 'NA' North Atlantic - 'EP' Eastern North Pacific - 'WP' Western North Pacific - 'NI' North Indian - 'SI' South Indian - 'SP' Southern Pacific - 'SA' South Atlantic - windfields : list of csr_matrix, optional - For each event, the full velocity vectors at each centroid and track position in a - sparse matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full - ndarray of shape (npositions, ncentroids, 2). - **kwargs : Hazard properties, optional - All other keyword arguments are passed to the Hazard constructor. - """ - kwargs.setdefault('haz_type', HAZ_TYPE) - Hazard.__init__(self, **kwargs) - self.category = category if category is not None else np.array([], int) - self.basin = basin if basin is not None else [] - self.windfields = windfields if windfields is not None else [] - - def set_from_tracks(self, *args, **kwargs): - """This function is deprecated, use TropCyclone.from_tracks instead.""" - LOGGER.warning("The use of TropCyclone.set_from_tracks is deprecated." - "Use TropCyclone.from_tracks instead.") - if "intensity_thres" not in kwargs: - # some users modify the threshold attribute before calling `set_from_tracks` - kwargs["intensity_thres"] = self.intensity_thres - if self.pool is not None and 'pool' not in kwargs: - kwargs['pool'] = self.pool - self.__dict__ = TropCyclone.from_tracks(*args, **kwargs).__dict__ - - @classmethod - def from_tracks( - cls, - tracks: TCTracks, - centroids: Centroids, - pool: Optional[pathos.pools.ProcessPool] = None, - model: str = 'H08', - model_kwargs: Optional[dict] = None, - ignore_distance_to_coast: bool = False, - store_windfields: bool = False, - metric: str = "equirect", - intensity_thres: float = DEF_INTENSITY_THRES, - max_latitude: float = 61, - max_dist_inland_km: float = 1000, - max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, - max_memory_gb: float = DEF_MAX_MEMORY_GB, - ): - """ - Create new TropCyclone instance that contains windfields from the specified tracks. - - This function sets the ``intensity`` attribute to contain, for each centroid, - the maximum wind speed (1-minute sustained winds at 10 meters above ground) - experienced over the whole period of each TC event in m/s. The wind speed is set - to 0 if it doesn't exceed the threshold ``intensity_thres``. - - The ``category`` attribute is set to the value of the ``category``-attribute - of each of the given track data sets. - - The ``basin`` attribute is set to the genesis basin for each event, which - is the first value of the ``basin``-variable in each of the given track data sets. - - Optionally, the time dependent, vectorial winds can be stored using the - ``store_windfields`` function parameter (see below). - - Parameters - ---------- - tracks : climada.hazard.TCTracks - Tracks of storm events. - centroids : Centroids, optional - Centroids where to model TC. Default: global centroids at 360 arc-seconds - resolution. - pool : pathos.pool, optional - Pool that will be used for parallel computation of wind fields. Default: - None - model : str, optional - Parametric wind field model to use. Default: "H08". - - * ``"H1980"`` (the prominent Holland 1980 model) from the paper: - Holland, G.J. (1980): An Analytic Model of the Wind and Pressure - Profiles in Hurricanes. Monthly Weather Review 108(8): 1212–1218. - ``https://doi.org/10.1175/1520-0493(1980)108<1212:AAMOTW>2.0.CO;2`` - * ``"H08"`` (Holland 1980 with b-value from Holland 2008) from the paper: - Holland, G. (2008). A revised hurricane pressure-wind model. Monthly - Weather Review, 136(9), 3432–3445. - https://doi.org/10.1175/2008MWR2395.1 - * ``"H10"`` (Holland et al. 2010) from the paper: - Holland et al. (2010): A Revised Model for Radial Profiles of - Hurricane Winds. Monthly Weather Review 138(12): 4393–4401. - https://doi.org/10.1175/2010MWR3317.1 - * ``"ER11"`` (Emanuel and Rotunno 2011) from the paper: - Emanuel, K., Rotunno, R. (2011): Self-Stratification of Tropical - Cyclone Outflow. Part I: Implications for Storm Structure. Journal - of the Atmospheric Sciences 68(10): 2236–2249. - https://dx.doi.org/10.1175/JAS-D-10-05024.1 - model_kwargs : dict, optional - If given, forward these kwargs to the selected wind model. None of the - parameters is currently supported by the ER11 model. Default: None. - The Holland models support the following parameters, in alphabetical order: - - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use. In H1980, the wind - profile is computed on the gradient level, and wind speeds are converted - to the surface level using this factor. In H08 and H10, the wind profile - is computed on the surface level, but the clipping interval of the - B-value depends on this factor. Default: 0.9 - rho_air_const : float or None, optional - The constant value for air density (in kg/m³) to assume in the formulas - from Holland 1980. By default, the constant value suggested in Holland - 1980 is used. If set to None, the air density is computed from pressure - following equation (9) in Holland et al. 2010. Default: 1.15 - vmax_from_cen : boolean, optional - Only used in H10. If True, replace the recorded value of vmax along the - track by an estimate from pressure, following equation (8) in Holland et - al. 2010. Default: True - vmax_in_brackets : bool, optional - Only used in H10. Specifies which of the two formulas in equation (6) of - Holland et al. 2010 to use. If False, the formula with vmax outside of - the brackets is used. Note that, a side-effect of the formula with vmax - inside of the brackets is that the wind speed maximum is attained a bit - farther away from the center than according to the recorded radius of - maximum winds (RMW). Default: False - - ignore_distance_to_coast : boolean, optional - If True, centroids far from coast are not ignored. - If False, the centroids' distances to the coast are calculated with the - `Centroids.get_dist_coast()` method (unless there is "dist_coast" column in - the centroids' GeoDataFrame) and centroids far from coast are ignored. - Default: False. - store_windfields : boolean, optional - If True, the Hazard object gets a list ``windfields`` of sparse matrices. - For each track, the full velocity vectors at each centroid and track - position are stored in a sparse matrix of shape (npositions, - ncentroids * 2) that can be reshaped to a full ndarray of shape (npositions, - ncentroids, 2). Default: False. - metric : str, optional - Specify an approximation method to use for earth distances: - - * "equirect": Distance according to sinusoidal projection. Fast, but - inaccurate for large distances and high latitudes. - * "geosphere": Exact spherical distance. Much more accurate at all - distances, but slow. - - Default: "equirect". - intensity_thres : float, optional - Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 - max_latitude : float, optional - No wind speed calculation is done for centroids with latitude larger than - this parameter. Default: 61 - max_dist_inland_km : float, optional - No wind speed calculation is done for centroids with a distance (in km) to - the coast larger than this parameter. Default: 1000 - max_dist_eye_km : float, optional - No wind speed calculation is done for centroids with a distance (in km) to - the TC center ("eye") larger than this parameter. Default: 300 - max_memory_gb : float, optional - To avoid memory issues, the computation is done for chunks of the track - sequentially. The chunk size is determined depending on the available memory - (in GB). Note that this limit applies to each thread separately if a - ``pool`` is used. Default: 8 - - Raises - ------ - ValueError - - Returns - ------- - TropCyclone - """ - num_tracks = tracks.size - - if ignore_distance_to_coast: - # Select centroids with lat <= max_latitude - [idx_centr_filter] = (np.abs(centroids.lat) <= max_latitude).nonzero() - else: - # Select centroids which are inside max_dist_inland_km and lat <= max_latitude - if 'dist_coast' not in centroids.gdf.columns: - dist_coast = centroids.get_dist_coast() - else: - dist_coast = centroids.gdf.dist_coast.values - [idx_centr_filter] = ( - (dist_coast <= max_dist_inland_km * 1000) - & (np.abs(centroids.lat) <= max_latitude) - ).nonzero() - - # Filter early with a larger threshold, but inaccurate (lat/lon) distances. - # Later, there will be another filtering step with more accurate distances in km. - max_dist_eye_deg = max_dist_eye_km / ( - u_const.ONE_LAT_KM * np.cos(np.radians(max_latitude)) - ) - - # Restrict to coastal centroids within reach of any of the tracks - t_lon_min, t_lat_min, t_lon_max, t_lat_max = tracks.get_bounds(deg_buffer=max_dist_eye_deg) - t_mid_lon = 0.5 * (t_lon_min + t_lon_max) - filtered_centroids = centroids.coord[idx_centr_filter] - u_coord.lon_normalize(filtered_centroids[:, 1], center=t_mid_lon) - idx_centr_filter = idx_centr_filter[ - (t_lon_min <= filtered_centroids[:, 1]) - & (filtered_centroids[:, 1] <= t_lon_max) - & (t_lat_min <= filtered_centroids[:, 0]) - & (filtered_centroids[:, 0] <= t_lat_max) - ] - - # prepare keyword arguments to pass to `from_single_track` - kwargs_from_single_track = dict( - centroids=centroids, - idx_centr_filter=idx_centr_filter, - model=model, - model_kwargs=model_kwargs, - store_windfields=store_windfields, - metric=metric, - intensity_thres=intensity_thres, - max_dist_eye_km=max_dist_eye_km, - max_memory_gb=max_memory_gb, - ) + si_track["vgrad"] = ( + np.fmax(0, si_track["vmax"] - si_track["vtrans_norm"]) / gradient_to_surface_winds + ) - LOGGER.info( - 'Mapping %d tracks to %d coastal centroids.', num_tracks, idx_centr_filter.size, - ) - if pool: - chunksize = max(min(num_tracks // pool.ncpus, 1000), 1) - kwargs_repeated = [ - itertools.repeat(val, num_tracks) - for val in kwargs_from_single_track.values() - ] - tc_haz_list = pool.map( - cls.from_single_track, - tracks.data, - *kwargs_repeated, - chunksize=chunksize, - ) - else: - last_perc = 0 - tc_haz_list = [] - for track in tracks.data: - perc = 100 * len(tc_haz_list) / len(tracks.data) - if perc - last_perc >= 10: - LOGGER.info("Progress: %d%%", perc) - last_perc = perc - tc_haz_list.append( - cls.from_single_track(track, **kwargs_from_single_track) - ) - if last_perc < 100: - LOGGER.info("Progress: 100%") - - LOGGER.debug('Concatenate events.') - haz = cls.concat(tc_haz_list) - haz.pool = pool - haz.intensity_thres = intensity_thres - LOGGER.debug('Compute frequency.') - haz.frequency_from_tracks(tracks.data) - return haz - - def apply_climate_scenario_knu( - self, - ref_year: int = 2050, - rcp_scenario: int = 45 - ): - """ - From current TC hazard instance, return new hazard set with - future events for a given RCP scenario and year based on the - parametrized values derived from Table 3 in Knutson et al 2015. - https://doi.org/10.1175/JCLI-D-15-0129.1 . The scaling for different - years and RCP scenarios is obtained by linear interpolation. - - Note: The parametrized values are derived from the overall changes - in statistical ensemble of tracks. Hence, this method should only be - applied to sufficiently large tropical cyclone event sets that - approximate the reference years 1981 - 2008 used in Knutson et. al. - - The frequency and intensity changes are applied independently from - one another. The mean intensity factors can thus slightly deviate - from the Knutson value (deviation was found to be less than 1% - for default IBTrACS event sets 1980-2020 for each basin). - - Parameters - ---------- - ref_year : int - year between 2000 ad 2100. Default: 2050 - rcp_scenario : int - 26 for RCP 2.6, 45 for RCP 4.5, 60 for RCP 6.0 and 85 for RCP 8.5. - The default is 45. - - Returns - ------- - haz_cc : climada.hazard.TropCyclone - Tropical cyclone with frequencies and intensity scaled according - to the Knutson criterion for the given year and RCP. Returns - a new instance of climada.hazard.TropCyclone, self is not - modified. - """ - chg_int_freq = get_knutson_criterion() - scale_rcp_year = calc_scale_knutson(ref_year, rcp_scenario) - haz_cc = self._apply_knutson_criterion(chg_int_freq, scale_rcp_year) - return haz_cc - - def set_climate_scenario_knu(self, *args, **kwargs): - """This function is deprecated, use TropCyclone.apply_climate_scenario_knu instead.""" - LOGGER.warning("The use of TropCyclone.set_climate_scenario_knu is deprecated." - "Use TropCyclone.apply_climate_scenario_knu instead.") - return self.apply_climate_scenario_knu(*args, **kwargs) - - @classmethod - def video_intensity( - cls, - track_name: str, - tracks: TCTracks, - centroids: Centroids, - file_name: Optional[str] = None, - writer: animation = animation.PillowWriter(bitrate=500), - figsize: Tuple[float, float] = (9, 13), - adapt_fontsize: bool = True, - **kwargs - ): - """ - Generate video of TC wind fields node by node and returns its - corresponding TropCyclone instances and track pieces. - - Parameters - ---------- - track_name : str - name of the track contained in tracks to record - tracks : climada.hazard.TCTracks - tropical cyclone tracks - centroids : climada.hazard.Centroids - centroids where wind fields are mapped - file_name : str, optional - file name to save video (including full path and file extension) - writer : matplotlib.animation.*, optional - video writer. Default is pillow with bitrate=500 - figsize : tuple, optional - figure size for plt.subplots - adapt_fontsize : bool, optional - If set to true, the size of the fonts will be adapted to the size of the figure. - Otherwise the default matplotlib font size is used. Default is True. - kwargs : optional - arguments for pcolormesh matplotlib function used in event plots - - Returns - ------- - tc_list, tc_coord : list(TropCyclone), list(np.ndarray) - - Raises - ------ - ValueError - - """ - # initialization - track = tracks.get_track(track_name) - if not track: - raise ValueError(f'{track_name} not found in track data.') - idx_plt = np.argwhere( - (track.lon.values < centroids.total_bounds[2] + 1) - & (centroids.total_bounds[0] - 1 < track.lon.values) - & (track.lat.values < centroids.total_bounds[3] + 1) - & (centroids.total_bounds[1] - 1 < track.lat.values) - ).reshape(-1) - - tc_list = [] - tr_coord = {'lat': [], 'lon': []} - for node in range(idx_plt.size - 2): - tr_piece = track.sel( - time=slice(track.time.values[idx_plt[node]], - track.time.values[idx_plt[node + 2]])) - tr_piece.attrs['n_nodes'] = 2 # plot only one node - tr_sel = TCTracks() - tr_sel.append(tr_piece) - tr_coord['lat'].append(tr_sel.data[0].lat.values[:-1]) - tr_coord['lon'].append(tr_sel.data[0].lon.values[:-1]) - - tc_tmp = cls.from_tracks(tr_sel, centroids=centroids) - tc_tmp.event_name = [ - track.name + ' ' + time.strftime( - "%d %h %Y %H:%M", - time.gmtime(tr_sel.data[0].time[1].values.astype(int) - / 1000000000) - ) - ] - tc_list.append(tc_tmp) - - if 'cmap' not in kwargs: - kwargs['cmap'] = 'Greys' - if 'vmin' not in kwargs: - kwargs['vmin'] = np.array([tc_.intensity.min() for tc_ in tc_list]).min() - if 'vmax' not in kwargs: - kwargs['vmax'] = np.array([tc_.intensity.max() for tc_ in tc_list]).max() - - def run(node): - tc_list[node].plot_intensity(1, axis=axis, **kwargs) - axis.plot(tr_coord['lon'][node], tr_coord['lat'][node], 'k') - axis.set_title(tc_list[node].event_name[0]) - pbar.update() - - if file_name: - LOGGER.info('Generating video %s', file_name) - fig, axis, _fontsize = u_plot.make_map(figsize=figsize, adapt_fontsize=adapt_fontsize) - pbar = tqdm(total=idx_plt.size - 2) - ani = animation.FuncAnimation(fig, run, frames=idx_plt.size - 2, - interval=500, blit=False) - fig.tight_layout() - ani.save(file_name, writer=writer) - pbar.close() - return tc_list, tr_coord - - def frequency_from_tracks(self, tracks: List): - """ - Set hazard frequency from tracks data. - - Parameters - ---------- - tracks : list of xarray.Dataset - """ - if not tracks: - return - year_max = np.amax([t.time.dt.year.values.max() for t in tracks]) - year_min = np.amin([t.time.dt.year.values.min() for t in tracks]) - year_delta = year_max - year_min + 1 - num_orig = np.count_nonzero(self.orig) - ens_size = (self.event_id.size / num_orig) if num_orig > 0 else 1 - self.frequency = np.ones(self.event_id.size) / (year_delta * ens_size) - - @classmethod - def from_single_track( - cls, - track: xr.Dataset, - centroids: Centroids, - idx_centr_filter: np.ndarray, - model: str = 'H08', - model_kwargs: Optional[dict] = None, - store_windfields: bool = False, - metric: str = "equirect", - intensity_thres: float = DEF_INTENSITY_THRES, - max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, - max_memory_gb: float = DEF_MAX_MEMORY_GB, - ): - """ - Generate windfield hazard from a single track dataset - - Parameters - ---------- - track : xr.Dataset - Single tropical cyclone track. - centroids : Centroids - Centroids instance. - idx_centr_filter : np.ndarray - Indices of centroids to restrict to (e.g. sufficiently close to coast). - model : str, optional - Parametric wind field model, one of "H1980" (the prominent Holland 1980 model), - "H08" (Holland 1980 with b-value from Holland 2008), "H10" (Holland et al. 2010), or - "ER11" (Emanuel and Rotunno 2011). - Default: "H08". - model_kwargs: dict, optional - If given, forward these kwargs to the selected model. Default: None - store_windfields : boolean, optional - If True, store windfields. Default: False. - metric : str, optional - Specify an approximation method to use for earth distances: "equirect" (faster) or - "geosphere" (more accurate). See ``dist_approx`` function in - ``climada.util.coordinates``. - Default: "equirect". - intensity_thres : float, optional - Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 - max_dist_eye_km : float, optional - No wind speed calculation is done for centroids with a distance (in km) to the TC - center ("eye") larger than this parameter. Default: 300 - max_memory_gb : float, optional - To avoid memory issues, the computation is done for chunks of the track sequentially. - The chunk size is determined depending on the available memory (in GB). Default: 8 - - Raises - ------ - ValueError, KeyError - - Returns - ------- - haz : TropCyclone - """ - intensity_sparse, windfields_sparse = _compute_windfields_sparse( - track=track, - centroids=centroids, - idx_centr_filter=idx_centr_filter, - model=model, - model_kwargs=model_kwargs, - store_windfields=store_windfields, - metric=metric, - intensity_thres=intensity_thres, - max_dist_eye_km=max_dist_eye_km, - max_memory_gb=max_memory_gb, - ) - new_haz = cls(haz_type=HAZ_TYPE) - new_haz.intensity_thres = intensity_thres - new_haz.intensity = intensity_sparse - if store_windfields: - new_haz.windfields = [windfields_sparse] - new_haz.units = 'm/s' - new_haz.centroids = centroids - new_haz.event_id = np.array([1]) - new_haz.frequency = np.array([1]) - new_haz.event_name = [track.sid] - new_haz.fraction = sparse.csr_matrix(new_haz.intensity.shape) - # store first day of track as date - new_haz.date = np.array([ - dt.datetime(track.time.dt.year.values[0], - track.time.dt.month.values[0], - track.time.dt.day.values[0]).toordinal() - ]) - new_haz.orig = np.array([track.orig_event_flag]) - new_haz.category = np.array([track.category]) - # users that pickle TCTracks objects might still have data with the legacy basin attribute, - # so we have to deal with it here - new_haz.basin = [track.basin if isinstance(track.basin, str) - else str(track.basin.values[0])] - return new_haz - - def _apply_knutson_criterion( - self, - chg_int_freq: List, - scaling_rcp_year: float - ): - """ - Apply changes to intensities and cumulative frequencies. - - Parameters - ---------- - chg_int_freq : list(dict)) - list of criteria from climada.hazard.tc_clim_change - scaling_rcp_year : float - scale parameter because of chosen year and RCP - - Returns - ------- - tc_cc : climada.hazard.TropCyclone - Tropical cyclone with frequency and intensity scaled inspired by - the Knutson criterion. Returns a new instance of TropCyclone. - """ - - tc_cc = copy.deepcopy(self) - - # Criterion per basin - for basin in np.unique(tc_cc.basin): - - bas_sel = np.array(tc_cc.basin) == basin - - # Apply intensity change - inten_chg = [chg - for chg in chg_int_freq - if (chg['variable'] == 'intensity' and - chg['basin'] == basin) - ] - for chg in inten_chg: - sel_cat_chg = np.isin(tc_cc.category, chg['category']) & bas_sel - inten_scaling = 1 + (chg['change'] - 1) * scaling_rcp_year - tc_cc.intensity = sparse.diags( - np.where(sel_cat_chg, inten_scaling, 1) - ).dot(tc_cc.intensity) - - # Apply frequency change - freq_chg = [chg - for chg in chg_int_freq - if (chg['variable'] == 'frequency' and - chg['basin'] == basin) - ] - freq_chg.sort(reverse=False, key=lambda x: len(x['category'])) - - # Scale frequencies by category - cat_larger_list = [] - for chg in freq_chg: - cat_chg_list = [cat - for cat in chg['category'] - if cat not in cat_larger_list - ] - sel_cat_chg = np.isin(tc_cc.category, cat_chg_list) & bas_sel - if sel_cat_chg.any(): - freq_scaling = 1 + (chg['change'] - 1) * scaling_rcp_year - tc_cc.frequency[sel_cat_chg] *= freq_scaling - cat_larger_list += cat_chg_list - - if (tc_cc.frequency < 0).any(): - raise ValueError("The application of the given climate scenario" - "resulted in at least one negative frequency.") - - return tc_cc - -def _compute_windfields_sparse( - track: xr.Dataset, - centroids: Centroids, - idx_centr_filter: np.ndarray, - model: str = 'H08', +def compute_angular_windspeeds( + si_track: xr.Dataset, + d_centr: np.ndarray, + mask_centr_close: np.ndarray, + model: int, model_kwargs: Optional[dict] = None, - store_windfields: bool = False, - metric: str = "equirect", - intensity_thres: float = DEF_INTENSITY_THRES, - max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, - max_memory_gb: float = DEF_MAX_MEMORY_GB, -) -> Tuple[sparse.csr_matrix, Optional[sparse.csr_matrix]]: - """Version of ``compute_windfields`` that returns sparse matrices and limits memory usage + cyclostrophic: bool = False, +): + """Compute (absolute) angular wind speeds according to a parametric wind profile Parameters ---------- - track : xr.Dataset - Single tropical cyclone track. - centroids : Centroids - Centroids instance. - idx_centr_filter : np.ndarray - Indices of centroids to restrict to (e.g. sufficiently close to coast). - model : str, optional - Parametric wind field model, one of "H1980" (the prominent Holland 1980 model), - "H08" (Holland 1980 with b-value from Holland 2008), "H10" (Holland et al. 2010), or - "ER11" (Emanuel and Rotunno 2011). - Default: "H08". + si_track : xr.Dataset + Output of ``tctrack_to_si``. Which data variables are used in the computation of the wind + profile depends on the selected model. + d_centr : np.ndarray of shape (npositions, ncentroids) + Distance (in m) between centroids and track positions. + mask_centr_close : np.ndarray of shape (npositions, ncentroids) + For each track position one row indicating which centroids are within reach. + model : int + Wind profile model selection according to MODEL_VANG. model_kwargs: dict, optional If given, forward these kwargs to the selected model. Default: None - store_windfields : boolean, optional - If True, store windfields. Default: False. - metric : str, optional - Specify an approximation method to use for earth distances: "equirect" (faster) or - "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. - Default: "equirect". - intensity_thres : float, optional - Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 - max_dist_eye_km : float, optional - No wind speed calculation is done for centroids with a distance (in km) to the TC - center ("eye") larger than this parameter. Default: 300 - max_memory_gb : float, optional - To avoid memory issues, the computation is done for chunks of the track sequentially. - The chunk size is determined depending on the available memory (in GB). Default: 8 - - Raises - ------ - ValueError + cyclostrophic : bool, optional + If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0). + Default: False Returns ------- - intensity : csr_matrix - Maximum wind speed in each centroid over the whole storm life time. - windfields : csr_matrix or None - If store_windfields is True, the full velocity vectors at each centroid and track position - are stored in a sparse matrix of shape (npositions, ncentroids * 2) that can be reshaped - to a full ndarray of shape (npositions, ncentroids, 2). - If store_windfields is False, None is returned. + ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ - try: - mod_id = MODEL_VANG[model] - except KeyError as err: - raise ValueError(f'Model not implemented: {model}.') from err - - ncentroids = centroids.coord.shape[0] - npositions = track.sizes["time"] - windfields_shape = (npositions, ncentroids * 2) - intensity_shape = (1, ncentroids) - - # initialise arrays for the assumption that no centroids are within reach - windfields_sparse = ( - sparse.csr_matrix(([], ([], [])), shape=windfields_shape) - if store_windfields else None + model_kwargs = {} if model_kwargs is None else model_kwargs + compute_funs = { + MODEL_VANG['H1980']: _compute_angular_windspeeds_h1980, + MODEL_VANG['H08']: _compute_angular_windspeeds_h08, + MODEL_VANG['H10']: _compute_angular_windspeeds_h10, + MODEL_VANG['ER11']: _stat_er_2011, + } + if model not in compute_funs: + raise NotImplementedError(f"The specified wind model is not supported: {model}") + result = compute_funs[model]( + si_track, + d_centr, + mask_centr_close, + cyclostrophic=cyclostrophic, + **model_kwargs, ) - intensity_sparse = sparse.csr_matrix(([], ([], [])), shape=intensity_shape) + result[0, :] *= 0 + return result - # The wind field model requires at least two track positions because translational speed - # as well as the change in pressure (in case of H08) are required. - if npositions < 2: - return intensity_sparse, windfields_sparse - # convert track variables to SI units - si_track = tctrack_to_si(track, metric=metric) +def _compute_angular_windspeeds_h1980( + si_track: xr.Dataset, + d_centr: np.ndarray, + close_centr_msk: np.ndarray, + cyclostrophic: bool = False, + gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, + rho_air_const: float = DEF_RHO_AIR, +): + """Compute (absolute) angular wind speeds according to the Holland 1980 model - # When done properly, finding and storing the close centroids is not a memory bottle neck and - # can be done before chunking. Note that the longitudinal coordinates of `centroids_close` as - # returned by `get_close_centroids` are normalized to be consistent with the coordinates in - # `si_track`. - centroids_close, mask_centr, mask_centr_alongtrack = get_close_centroids( - si_track, centroids.coord[idx_centr_filter], max_dist_eye_km, metric=metric, - ) - idx_centr_filter = idx_centr_filter[mask_centr] - n_centr_close = centroids_close.shape[0] - if n_centr_close == 0: - return intensity_sparse, windfields_sparse + Parameters + ---------- + si_track : xr.Dataset + Output of `tctrack_to_si`. Which data variables are used in the computation of the wind + profile depends on the selected model. + d_centr : np.ndarray of shape (npositions, ncentroids) + Distance (in m) between centroids and track positions. + close_centr_msk : np.ndarray of shape (npositions, ncentroids) + For each track position one row indicating which centroids are within reach. + cyclostrophic : bool, optional + If True, don't apply the influence of the Coriolis force (set the Coriolis terms to 0). + Default: False + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use. The wind profile is computed on the + gradient level, and wind speeds are converted to the surface level using this factor. + Default: 0.9 + rho_air_const : float or None, optional + The constant value for air density (in kg/m³) to assume in the formulas from Holland 1980. + By default, the constant value suggested in Holland 1980 is used. If set to None, the air + density is computed from pressure following equation (9) in Holland et al. 2010. + Default: 1.15 - # the total memory requirement in GB if we compute everything without chunking: - # 8 Bytes per entry (float64), 10 arrays - total_memory_gb = npositions * n_centr_close * 8 * 10 / 1e9 - if total_memory_gb > max_memory_gb and npositions > 2: - # If the number of positions is down to 2 already, we cannot split any further. In that - # case, we just take the risk and try to do the computation anyway. It might still work - # since we have only computed an upper bound for the number of affected centroids. + Returns + ------- + ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location + """ + _vgrad(si_track, gradient_to_surface_winds) + _rho_air(si_track, rho_air_const) + _B_holland_1980(si_track) + result = _stat_holland_1980(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) + result *= gradient_to_surface_winds + return result - # Split the track into chunks, compute the result for each chunk, and combine: - return _compute_windfields_sparse_chunked( - mask_centr_alongtrack, - track, - centroids, - idx_centr_filter, - model=model, - model_kwargs=model_kwargs, - store_windfields=store_windfields, - metric=metric, - intensity_thres=intensity_thres, - max_dist_eye_km=max_dist_eye_km, - max_memory_gb=max_memory_gb, - ) - windfields, idx_centr_reachable = _compute_windfields( - si_track, - centroids_close, - mod_id, - model_kwargs=model_kwargs, - metric=metric, - max_dist_eye_km=max_dist_eye_km, - ) - idx_centr_filter = idx_centr_filter[idx_centr_reachable] - npositions = windfields.shape[0] - - intensity = np.linalg.norm(windfields, axis=-1).max(axis=0) - intensity[intensity < intensity_thres] = 0 - intensity_sparse = sparse.csr_matrix( - (intensity, idx_centr_filter, [0, intensity.size]), - shape=intensity_shape) - intensity_sparse.eliminate_zeros() - - windfields_sparse = None - if store_windfields: - n_centr_filter = idx_centr_filter.size - indices = np.zeros((npositions, n_centr_filter, 2), dtype=np.int64) - indices[:, :, 0] = 2 * idx_centr_filter[None] - indices[:, :, 1] = 2 * idx_centr_filter[None] + 1 - indices = indices.ravel() - indptr = np.arange(npositions + 1) * n_centr_filter * 2 - windfields_sparse = sparse.csr_matrix((windfields.ravel(), indices, indptr), - shape=windfields_shape) - windfields_sparse.eliminate_zeros() - - return intensity_sparse, windfields_sparse - -def _compute_windfields_sparse_chunked( - mask_centr_alongtrack: np.ndarray, - track: xr.Dataset, - *args, - max_memory_gb: float = DEF_MAX_MEMORY_GB, - **kwargs, -) -> Tuple[sparse.csr_matrix, Optional[sparse.csr_matrix]]: - """Call ``_compute_windfields_sparse`` for chunks of the track and re-assemble the results +def _compute_angular_windspeeds_h08( + si_track: xr.Dataset, + d_centr: np.ndarray, + close_centr_msk: np.ndarray, + cyclostrophic: bool = False, + gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, + rho_air_const: float = DEF_RHO_AIR, +): + """Compute (absolute) angular wind speeds according to the Holland 2008 model Parameters ---------- - mask_centr_alongtrack : np.ndarray of shape (npositions, ncentroids) - Each row is a mask that indicates the centroids within reach for one track position. - track : xr.Dataset - Single tropical cyclone track. - max_memory_gb : float, optional - Maximum memory requirements (in GB) for the computation of a single chunk of the track. - Default: 8 - args, kwargs : - The remaining arguments are passed on to ``_compute_windfields_sparse``. + si_track : xr.Dataset + Output of `tctrack_to_si`. Which data variables are used in the computation of the wind + profile depends on the selected model. + d_centr : np.ndarray of shape (npositions, ncentroids) + Distance (in m) between centroids and track positions. + close_centr_msk : np.ndarray of shape (npositions, ncentroids) + For each track position one row indicating which centroids are within reach. + cyclostrophic : bool, optional + If True, don't apply the influence of the Coriolis force (set the Coriolis terms to 0). + Default: False + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use. The wind profile is computed on the + surface level, but the clipping interval of the B-value depends on this factor. + Default: 0.9 + rho_air_const : float or None, optional + The constant value for air density (in kg/m³) to assume in the formula from Holland 1980. + By default, the constant value suggested in Holland 1980 is used. If set to None, the air + density is computed from pressure following equation (9) in Holland et al. 2010. + Default: 1.15 Returns ------- - intensity, windfields : - See ``_compute_windfields_sparse`` for a description of the return values. + ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ - npositions = track.sizes["time"] - # The memory requirements for each track position are estimated for the case of 10 arrays - # containing `nreachable` float64 (8 Byte) values each. The chunking is only relevant in - # extreme cases with a very high temporal and/or spatial resolution. - max_nreachable = max_memory_gb * 1e9 / (8 * 10 * npositions) - split_pos = [0] - chunk_size = 2 - while split_pos[-1] + chunk_size < npositions: - chunk_size += 1 - # create overlap between consecutive chunks - chunk_start = max(0, split_pos[-1] - 1) - chunk_end = chunk_start + chunk_size - nreachable = mask_centr_alongtrack[chunk_start:chunk_end].any(axis=0).sum() - if nreachable > max_nreachable: - split_pos.append(chunk_end - 1) - chunk_size = 2 - split_pos.append(npositions) - - intensity = [] - windfields = [] - for prev_chunk_end, chunk_end in zip(split_pos[:-1], split_pos[1:]): - chunk_start = max(0, prev_chunk_end - 1) - inten, win = _compute_windfields_sparse( - track.isel(time=slice(chunk_start, chunk_end)), *args, - max_memory_gb=max_memory_gb, **kwargs, - ) - intensity.append(inten) - windfields.append(win) + _rho_air(si_track, rho_air_const) + _bs_holland_2008(si_track, gradient_to_surface_winds=gradient_to_surface_winds) + return _stat_holland_1980(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) - intensity = sparse.csr_matrix(sparse.vstack(intensity).max(axis=0)) - if windfields[0] is not None: - # eliminate the overlap between consecutive chunks - windfields = [windfields[0]] + [win[1:, :] for win in windfields[1:]] - windfields = sparse.vstack(windfields, format="csr") - return intensity, windfields -def _compute_windfields( +def _compute_angular_windspeeds_h10( si_track: xr.Dataset, - centroids: np.ndarray, - model: int, - model_kwargs: Optional[dict] = None, - metric: str = "equirect", - max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, -) -> Tuple[np.ndarray, np.ndarray]: - """Compute 1-minute sustained winds (in m/s) at 10 meters above ground + d_centr: np.ndarray, + close_centr_msk: np.ndarray, + cyclostrophic: bool = True, + gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, + rho_air_const: float = DEF_RHO_AIR, + vmax_from_cen: bool = True, + vmax_in_brackets: bool = False, +): + """Compute (absolute) angular wind speeds according to the Holland et al. 2010 model - In a first step, centroids within reach of the track are determined so that wind fields will - only be computed and returned for those centroids. Still, since computing the distance of - the storm center to the centroids is computationally expensive, make sure to pre-filter the - centroids and call this function only for those centroids that are potentially affected. + Note that this model is always cyclostrophic, the parameter setting is ignored. Parameters ---------- si_track : xr.Dataset - Output of ``tctrack_to_si``. Which data variables are used in the computation of the wind - speeds depends on the selected model. - centroids : np.ndarray with two dimensions - Each row is a centroid [lat, lon]. Centroids that are not within reach of the track are - ignored. Longitudinal coordinates are assumed to be normalized consistently with the - longitudinal coordinates in ``si_track``. - model : int - Wind profile model selection according to MODEL_VANG. - model_kwargs: dict, optional - If given, forward these kwargs to the selected model. Default: None - metric : str, optional - Specify an approximation method to use for earth distances: "equirect" (faster) or - "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. - Default: "equirect". - max_dist_eye_km : float, optional - No wind speed calculation is done for centroids with a distance (in km) to the TC center - ("eye") larger than this parameter. Default: 300 + Output of `tctrack_to_si`. Which data variables are used in the computation of the wind + profile depends on the selected model. + d_centr : np.ndarray of shape (npositions, ncentroids) + Distance (in m) between centroids and track positions. + close_centr_msk : np.ndarray of shape (npositions, ncentroids) + For each track position one row indicating which centroids are within reach. + cyclostrophic : bool, optional + This parameter is ignored because this model is always cyclostrophic. Default: True + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use. the wind profile is computed on the + surface level, but the clipping interval of the B-value depends on this factor. + Default: 0.9 + rho_air_const : float or None, optional + The constant value for air density (in kg/m³) to assume in the formula for the B-value. By + default, the value suggested in Holland 1980 is used. If set to None, the air density is + computed from pressure following equation (9) in Holland et al. 2010. Default: 1.15 + vmax_from_cen : boolean, optional + If True, replace the recorded value of vmax along the track by an estimate from pressure, + following equation (8) in Holland et al. 2010. Default: True + vmax_in_brackets : bool, optional + Specifies which of the two formulas in equation (6) of Holland et al. 2010 to use. If + False, the formula with vmax outside of the brackets is used. Note that, a side-effect of + the formula with vmax inside of the brackets is that the wind speed maximum is attained a + bit farther away from the center than according to the recorded radius of maximum + winds (RMW). Default: False Returns ------- - windfields : np.ndarray of shape (npositions, nreachable, 2) - Directional wind fields for each track position on those centroids within reach - of the TC track. Note that the wind speeds at the first position are all zero because - the discrete time derivatives involved in the process are implemented using backward - differences. However, the first position is usually not relevant for impact calculations - since it is far off shore. - idx_centr_reachable : np.ndarray of shape (nreachable,) - List of indices of input centroids within reach of the TC track. + ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ - # start with the assumption that no centroids are within reach - npositions = si_track.sizes["time"] - idx_centr_reachable = np.zeros((0,), dtype=np.int64) - windfields = np.zeros((npositions, 0, 2), dtype=np.float64) + if not cyclostrophic: + LOGGER.debug( + 'The function _compute_angular_windspeeds_h10 was called with parameter ' + '"cyclostrophic" equal to false. Please be aware that this setting is ignored as the' + ' Holland et al. 2010 model is always cyclostrophic.') + _rho_air(si_track, rho_air_const) + if vmax_from_cen: + _bs_holland_2008(si_track, gradient_to_surface_winds=gradient_to_surface_winds) + _v_max_s_holland_2008(si_track) + else: + _B_holland_1980(si_track, gradient_to_surface_winds=gradient_to_surface_winds) + hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk, vmax_in_brackets=vmax_in_brackets) + return _stat_holland_2010( + si_track, d_centr, close_centr_msk, hol_x, vmax_in_brackets=vmax_in_brackets, + ) - # compute distances (in m) and vectors to all centroids - [d_centr], [v_centr_normed] = u_coord.dist_approx( - si_track["lat"].values[None], si_track["lon"].values[None], - centroids[None, :, 0], centroids[None, :, 1], - log=True, normalize=False, method=metric, units="m") - # exclude centroids that are too far from or too close to the eye - mask_centr_close = (d_centr <= max_dist_eye_km * KM_TO_M) & (d_centr > 1) - if not np.any(mask_centr_close): - return windfields, idx_centr_reachable +def _rho_air(si_track: xr.Dataset, const: Optional[float]): + """Eyewall density of air (in kg/m³) at each track node. - # restrict to the centroids that are within reach of any of the positions - mask_centr_close_any = mask_centr_close.any(axis=0) - mask_centr_close = mask_centr_close[:, mask_centr_close_any] - d_centr = d_centr[:, mask_centr_close_any] - v_centr_normed = v_centr_normed[:, mask_centr_close_any, :] + The track dataset is modified in place, with the "rho_air" data variable added. - # normalize the vectors pointing from the eye to the centroids - v_centr_normed[~mask_centr_close] = 0 - v_centr_normed[mask_centr_close] /= d_centr[mask_centr_close, None] + Parameters + ---------- + si_track : xr.Dataset + Track information as returned by `tctrack_to_si`. The data variables used by this function + are "lat", "cen", and "pdelta". The result is stored in place as new data + variable "rho_air". + const : float or None + A constant value for air density (in kg/m³) to assume. If None, the air density is + estimated from eyewall pressure following equation (9) in Holland et al. 2010. + """ + if const is not None: + si_track["rho_air"] = xr.full_like(si_track["time"], const, dtype=float) + return - # derive (absolute) angular velocity from parametric wind profile - v_ang_norm = compute_angular_windspeeds( - si_track, d_centr, mask_centr_close, model, model_kwargs=model_kwargs, cyclostrophic=False, - ) + # surface relative humidity (unitless), assumed constant following Holland et al. 2010 + surface_relative_humidity = 0.9 - # Influence of translational speed decreases with distance from eye. - # The "absorbing factor" is according to the following paper (see Fig. 7): - # - # Mouton, F. & Nordbeck, O. (2005). Cyclone Database Manager. A tool - # for converting point data from cyclone observations into tracks and - # wind speed profiles in a GIS. UNED/GRID-Geneva. - # https://unepgrid.ch/en/resource/19B7D302 - # - t_rad_bc = np.broadcast_to(si_track["rad"].values[:, None], d_centr.shape) - v_trans_corr = np.zeros_like(d_centr) - v_trans_corr[mask_centr_close] = np.fmin( - 1, t_rad_bc[mask_centr_close] / d_centr[mask_centr_close]) + # surface temperature (in °C), following equation (9) in Holland 2008 + temp_s = 28.0 - 3.0 * (si_track["lat"] - 10.0) / 20.0 - if model in [MODEL_VANG['H08'], MODEL_VANG['H10']]: - # In these models, v_ang_norm already contains vtrans_norm, so subtract it first, before - # converting to vectors and then adding (vectorial) vtrans again. Make sure to apply the - # "absorbing factor" in both steps: - vtrans_norm_bc = np.broadcast_to(si_track["vtrans_norm"].values[:, None], d_centr.shape) - v_ang_norm[mask_centr_close] -= ( - vtrans_norm_bc[mask_centr_close] * v_trans_corr[mask_centr_close] - ) + # eyewall surface pressure (in Pa), following equation (6) in Holland 2008 + pres_eyewall = si_track["cen"] + si_track["pdelta"] / np.exp(1) - # vectorial angular velocity - windfields = ( - si_track.attrs["latsign"] * np.array([1.0, -1.0])[..., :] * v_centr_normed[:, :, ::-1] - ) - windfields[mask_centr_close] *= v_ang_norm[mask_centr_close, None] + # mixing ratio (in kg/kg), estimated from temperature, using formula for saturation vapor + # pressure in Bolton 1980 (multiplied by the ratio of molar masses of water vapor and dry air) + # We multiply by 100, since the formula by Bolton is in hPa (mbar), and we use Pa. + r_mix = 100 * 3.802 / pres_eyewall * np.exp(17.67 * temp_s / (243.5 + temp_s)) - # add angular and corrected translational velocity vectors - windfields[1:] += si_track["vtrans"].values[1:, None, :] * v_trans_corr[1:, :, None] - windfields[np.isnan(windfields)] = 0 - windfields[0, :, :] = 0 - [idx_centr_reachable] = mask_centr_close_any.nonzero() - return windfields, idx_centr_reachable + # virtual surface temperature (in K) + temp_vs = (T_ICE_K + temp_s) * (1 + 0.81 * surface_relative_humidity * r_mix) -def tctrack_to_si( - track: xr.Dataset, - metric: str = "equirect", -) -> xr.Dataset: - """Convert track variables to SI units and prepare for wind field computation + # specific gas constant of dry air (in J/kgK) + r_dry_air = 286.9 - In addition to unit conversion, the variable names are shortened, the longitudinal coordinates - are normalized and additional variables are defined: + # density of air (in kg/m³); when checking the units, note that J/Pa = m³ + si_track["rho_air"] = pres_eyewall / (r_dry_air * temp_vs) - * cp (coriolis parameter) - * pdelta (difference between environmental and central pressure, always strictly positive) - * vtrans (translational velocity vectors) - * vtrans_norm (absolute value of translational speed) - Furthermore, some scalar variables are stored as attributes: +def _bs_holland_2008( + si_track: xr.Dataset, + gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, +): + """Holland's 2008 b-value estimate for sustained surface winds. + (This is also one option of how to estimate the b-value in Holland 2010, + for the other option consult '_bs_holland_2010_v2'.) - * latsign (1.0 if the track is located on the northern and -1.0 if on southern hemisphere) - * mid_lon (the central longitude that was used to normalize the longitudinal coordinates) + The result is stored in place as a new data variable "hol_b". - Finally, some corrections are applied to variables: + Unlike the original 1980 formula (see ``_B_holland_1980``), this approach does not require any + wind speed measurements, but is based on the more reliable pressure information. - * clip central pressure values so that environmental pressure values are never exceeded - * extrapolate radius of max wind from pressure if missing + The parameter applies to 1-minute sustained winds at 10 meters above ground. + It is taken from equation (11) in the following paper: - Parameters - ---------- - track : xr.Dataset - Track information. - metric : str, optional - Specify an approximation method to use for earth distances: "equirect" (faster) or - "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. - Default: "equirect". + Holland, G. (2008). A revised hurricane pressure-wind model. Monthly + Weather Review, 136(9), 3432–3445. https://doi.org/10.1175/2008MWR2395.1 - Returns - ------- - xr.Dataset - """ - si_track = track[["lat", "lon", "time"]].copy() - si_track["tstep"] = track["time_step"] * H_TO_S - si_track["env"] = track["environmental_pressure"] * MBAR_TO_PA + For reference, it reads - # we support some non-standard unit names - unit_replace = {"mb": "mbar", "kn": "knots"} - configs = [ - ("central_pressure", "cen", "Pa"), - ("max_sustained_wind", "vmax", "m/s"), - ] - for long_name, var_name, si_unit in configs: - unit = track.attrs[f"{long_name}_unit"] - unit = unit_replace.get(unit, unit) - try: - conv_factor = ureg(unit).to(si_unit).magnitude - except Exception as ex: - raise ValueError( - f"The {long_name}_unit '{unit}' in the provided track is not supported." - ) from ex - si_track[var_name] = track[long_name] * conv_factor + b_s = -4.4 * 1e-5 * (penv - pcen)^2 + 0.01 * (penv - pcen) + + 0.03 * (dp/dt) - 0.014 * |lat| + 0.15 * (v_trans)^hol_xx + 1.0 - # normalize longitudinal coordinates - si_track.attrs["mid_lon"] = 0.5 * sum(u_coord.lon_bounds(si_track["lon"].values)) - u_coord.lon_normalize(si_track["lon"].values, center=si_track.attrs["mid_lon"]) + where ``dp/dt`` is the time derivative of central pressure and ``hol_xx`` is Holland's x + parameter: hol_xx = 0.6 * (1 - (penv - pcen) / 215) - # make sure that central pressure never exceeds environmental pressure - pres_exceed_msk = (si_track["cen"] > si_track["env"]).values - si_track["cen"].values[pres_exceed_msk] = si_track["env"].values[pres_exceed_msk] + The equation for b_s has been fitted statistically using hurricane best track records for + central pressure and maximum wind. It therefore performs best in the North Atlantic. - # extrapolate radius of max wind from pressure if not given - si_track["rad"] = track["radius_max_wind"].copy() - si_track["rad"].values[:] = estimate_rmw( - si_track["rad"].values, si_track["cen"].values / MBAR_TO_PA, + Furthermore, b_s has been fitted under the assumption of a "cyclostrophic" wind field which + means that the influence from Coriolis forces is assumed to be small. This is reasonable close + to the radius of maximum wind where the Coriolis term (r*f/2) is small compared to the rest + (see ``_stat_holland_1980``). More precisely: At the radius of maximum wind speeds, the typical + order of the Coriolis term is 1 while wind speed is 50 (which changes away from the + radius of maximum winds and as the TC moves away from the equator). + + Parameters + ---------- + si_track : xr.Dataset + Output of ``tctrack_to_si``. The data variables used by this function are "lat", "tstep", + "vtrans_norm", "cen", and "pdelta". The result is stored in place as a new data + variable "hol_b". + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use when determining the clipping + interval. Default: 0.9 + """ + # adjust pressure at previous track point + prev_cen = np.zeros_like(si_track["cen"].values) + prev_cen[1:] = si_track["cen"].values[:-1].copy() + msk = prev_cen < 850 * MBAR_TO_PA + prev_cen[msk] = si_track["cen"].values[msk] + + # The formula assumes that pressure values are in millibar (hPa) instead of SI units (Pa), + # and time steps are in hours instead of seconds, but translational wind speed is still + # expected to be in m/s. + pdelta = si_track["pdelta"] / MBAR_TO_PA + hol_xx = 0.6 * (1. - pdelta / 215) + si_track["hol_b"] = ( + -4.4e-5 * pdelta**2 + 0.01 * pdelta + + 0.03 * (si_track["cen"] - prev_cen) / si_track["tstep"] * (H_TO_S / MBAR_TO_PA) + - 0.014 * abs(si_track["lat"]) + + 0.15 * si_track["vtrans_norm"]**hol_xx + 1.0 ) - si_track["rad"] *= NM_TO_KM * KM_TO_M + clip_interval = _b_holland_clip_interval(gradient_to_surface_winds) + si_track["hol_b"] = np.clip(si_track["hol_b"], *clip_interval) - hemisphere = 'N' - if np.count_nonzero(si_track["lat"] < 0) > np.count_nonzero(si_track["lat"] > 0): - hemisphere = 'S' - si_track.attrs["latsign"] = 1.0 if hemisphere == 'N' else -1.0 - # add translational speed of track at every node (in m/s) - _vtrans(si_track, metric=metric) +def _v_max_s_holland_2008(si_track: xr.Dataset): + """Compute maximum surface winds from pressure according to Holland 2008. - # add Coriolis parameter - si_track["cp"] = ("time", _coriolis_parameter(si_track["lat"].values)) + The result is stored in place as a data variable "vmax". If a variable of that name already + exists, its values are overwritten. - # add pressure drop - si_track["pdelta"] = np.fmax(np.spacing(1), si_track["env"] - si_track["cen"]) + This function implements equation (11) in the following paper: - return si_track + Holland, G. (2008). A revised hurricane pressure-wind model. Monthly + Weather Review, 136(9), 3432–3445. https://doi.org/10.1175/2008MWR2395.1 -def _vgrad(si_track, gradient_to_surface_winds): - """Gradient wind speeds (in m/s) without translational influence at each track node + For reference, it reads - The track dataset is modified in place, with the "vgrad" data variable added. + v_ms = [b_s / (rho * e) * (penv - pcen)]^0.5 + + where ``b_s`` is Holland b-value (see ``_bs_holland_2008``), e is Euler's number, rho is the + density of air, ``penv`` is environmental, and ``pcen`` is central pressure. Parameters ---------- si_track : xr.Dataset - Track information as returned by `tctrack_to_si`. The data variables used by this function - are "vmax" and "vtrans_norm". The result is stored in place as new data variable "vgrad". - gradient_to_surface_winds : float - The gradient-to-surface wind reduction factor to use. + Output of ``tctrack_to_si`` with "hol_b" (see _bs_holland_2008) and "rho_air" (see + _rho_air) variables. The data variables used by this function are "pdelta", "hol_b", + and "rho_air". The results are stored in place as a new data variable "vmax". If a variable + of that name already exists, its values are overwritten. """ - si_track["vgrad"] = ( - np.fmax(0, si_track["vmax"] - si_track["vtrans_norm"]) / gradient_to_surface_winds + si_track["vmax"] = np.sqrt( + si_track["hol_b"] / (si_track["rho_air"] * np.exp(1)) * si_track["pdelta"] ) -def compute_angular_windspeeds( + +def _B_holland_1980( # pylint: disable=invalid-name si_track: xr.Dataset, - d_centr: np.ndarray, - mask_centr_close: np.ndarray, - model: int, - model_kwargs: Optional[dict] = None, - cyclostrophic: bool = False, + gradient_to_surface_winds: Optional[float] = None, ): - """Compute (absolute) angular wind speeds according to a parametric wind profile + """Holland's 1980 B-value computation for gradient-level winds. + + The result is stored in place as a new data variable "hol_b". + + The parameter applies to gradient-level winds (about 1000 metres above the earth's surface). + The formula for B is derived from equations (5) and (6) in the following paper: + + Holland, G.J. (1980): An Analytic Model of the Wind and Pressure Profiles + in Hurricanes. Monthly Weather Review 108(8): 1212–1218. + https://doi.org/10.1175/1520-0493(1980)108<1212:AAMOTW>2.0.CO;2 + + For reference, inserting (6) into (5) and solving for B at r = RMW yields: + + B = v^2 * e * rho / (penv - pcen) + + where v are maximum gradient-level winds ``gradient_winds``, e is Euler's number, rho is the + density of air, ``penv`` is environmental, and ``pcen`` is central pressure. Parameters ---------- si_track : xr.Dataset - Output of ``tctrack_to_si``. Which data variables are used in the computation of the wind - profile depends on the selected model. - d_centr : np.ndarray of shape (npositions, ncentroids) - Distance (in m) between centroids and track positions. - mask_centr_close : np.ndarray of shape (npositions, ncentroids) - For each track position one row indicating which centroids are within reach. - model : int - Wind profile model selection according to MODEL_VANG. - model_kwargs: dict, optional - If given, forward these kwargs to the selected model. Default: None - cyclostrophic : bool, optional - If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0). - Default: False - - Returns - ------- - ndarray of shape (npositions, ncentroids) - containing the magnitude of the angular windspeed per track position per centroid location + Output of ``tctrack_to_si`` with "rho_air" variable (see _rho_air). The data variables + used by this function are "vgrad" (or "vmax" if gradient_to_surface_winds is different from + 1.0), "pdelta", and "rho_air". The results are stored in place as a new data + variable "hol_b". + gradient_to_surface_winds : float, optional + The gradient-to-surface wind reduction factor to use when determining the clipping + interval. By default, the gradient level values are assumed. Default: None """ - model_kwargs = {} if model_kwargs is None else model_kwargs - compute_funs = { - MODEL_VANG['H1980']: _compute_angular_windspeeds_h1980, - MODEL_VANG['H08']: _compute_angular_windspeeds_h08, - MODEL_VANG['H10']: _compute_angular_windspeeds_h10, - MODEL_VANG['ER11']: _stat_er_2011, - } - if model not in compute_funs: - raise NotImplementedError(f"The specified wind model is not supported: {model}") - result = compute_funs[model]( - si_track, - d_centr, - mask_centr_close, - cyclostrophic=cyclostrophic, - **model_kwargs, + windvar = "vgrad" if gradient_to_surface_winds is None else "vmax" + + si_track["hol_b"] = ( + si_track[windvar]**2 * np.exp(1) * si_track["rho_air"] / si_track["pdelta"] ) - result[0, :] *= 0 - return result -def _compute_angular_windspeeds_h1980( - si_track: xr.Dataset, - d_centr: np.ndarray, - close_centr_msk: np.ndarray, - cyclostrophic: bool = False, - gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, - rho_air_const: float = DEF_RHO_AIR, -): - """Compute (absolute) angular wind speeds according to the Holland 1980 model + clip_interval = _b_holland_clip_interval(gradient_to_surface_winds) + si_track["hol_b"] = np.clip(si_track["hol_b"], *clip_interval) + + +def _b_holland_clip_interval(gradient_to_surface_winds): + """The clip interval to use for the Holland B-value + + The default clip interval for gradient level B-values is taken to be (1.0, 2.5), following + Holland 1980. Parameters ---------- - si_track : xr.Dataset - Output of `tctrack_to_si`. Which data variables are used in the computation of the wind - profile depends on the selected model. - d_centr : np.ndarray of shape (npositions, ncentroids) - Distance (in m) between centroids and track positions. - close_centr_msk : np.ndarray of shape (npositions, ncentroids) - For each track position one row indicating which centroids are within reach. - cyclostrophic : bool, optional - If True, don't apply the influence of the Coriolis force (set the Coriolis terms to 0). - Default: False - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use. The wind profile is computed on the - gradient level, and wind speeds are converted to the surface level using this factor. - Default: 0.9 - rho_air_const : float or None, optional - The constant value for air density (in kg/m³) to assume in the formulas from Holland 1980. - By default, the constant value suggested in Holland 1980 is used. If set to None, the air - density is computed from pressure following equation (9) in Holland et al. 2010. - Default: 1.15 + gradient_to_surface_winds : float or None + The gradient-to-surface wind reduction factor to use when rescaling the gradient-level + clip interval (1.0, 2.5) proposed in Holland 1980. If None, no rescaling is applied. Returns ------- - ndarray of shape (npositions, ncentroids) - containing the magnitude of the angular windspeed per track position per centroid location + b_min, b_max : float + Minimum and maximum value of the clip interval. """ - _vgrad(si_track, gradient_to_surface_winds) - _rho_air(si_track, rho_air_const) - _B_holland_1980(si_track) - result = _stat_holland_1980(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) - result *= gradient_to_surface_winds - return result + clip_interval = (1.0, 2.5) + if gradient_to_surface_winds is not None: + fact = gradient_to_surface_winds**2 + clip_interval = tuple(c * fact for c in clip_interval) + return clip_interval -def _compute_angular_windspeeds_h08( + +def _x_holland_2010( si_track: xr.Dataset, d_centr: np.ndarray, - close_centr_msk: np.ndarray, - cyclostrophic: bool = False, - gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, - rho_air_const: float = DEF_RHO_AIR, -): - """Compute (absolute) angular wind speeds according to the Holland 2008 model + mask_centr_close: np.ndarray, + v_n: Union[float, np.ndarray] = 17.0, + r_n_km: Union[float, np.ndarray] = 300.0, + vmax_in_brackets: bool = False, +) -> np.ndarray: + """Compute exponent for wind model according to Holland et al. 2010. + + This function implements equation (10) from the following paper: + + Holland et al. (2010): A Revised Model for Radial Profiles of Hurricane Winds. Monthly + Weather Review 138(12): 4393–4401. https://doi.org/10.1175/2010MWR3317.1 + + For reference, it reads + + x = 0.5 [for r < r_max] + x = 0.5 + (r - r_max) * (x_n - 0.5) / (r_n - r_max) [for r >= r_max] + + The peripheral exponent x_n is adjusted to fit the peripheral observation of wind speeds ``v_n`` + at radius ``r_n``. Parameters ---------- si_track : xr.Dataset - Output of `tctrack_to_si`. Which data variables are used in the computation of the wind - profile depends on the selected model. - d_centr : np.ndarray of shape (npositions, ncentroids) - Distance (in m) between centroids and track positions. - close_centr_msk : np.ndarray of shape (npositions, ncentroids) - For each track position one row indicating which centroids are within reach. - cyclostrophic : bool, optional - If True, don't apply the influence of the Coriolis force (set the Coriolis terms to 0). - Default: False - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use. The wind profile is computed on the - surface level, but the clipping interval of the B-value depends on this factor. - Default: 0.9 - rho_air_const : float or None, optional - The constant value for air density (in kg/m³) to assume in the formula from Holland 1980. - By default, the constant value suggested in Holland 1980 is used. If set to None, the air - density is computed from pressure following equation (9) in Holland et al. 2010. - Default: 1.15 + Output of ``tctrack_to_si`` with "hol_b" variable (see _bs_holland_2008). The data variables + used by this function are "rad", "vmax", and "hol_b". + d_centr : np.ndarray of shape (nnodes, ncentroids) + Distance (in m) between centroids and track nodes. + mask_centr_close : np.ndarray of shape (nnodes, ncentroids) + Mask indicating for each track node which centroids are within reach of the windfield. + v_n : np.ndarray of shape (nnodes,) or float, optional + Peripheral wind speeds (in m/s) at radius ``r_n`` outside of radius of maximum winds + ``r_max``. In absence of a second wind speed measurement, this value defaults to 17 m/s + following Holland et al. 2010 (at a radius of 300 km). + r_n_km : np.ndarray of shape (nnodes,) or float, optional + Radius (in km) where the peripheral wind speed ``v_n`` is measured (or assumed). + In absence of a second wind speed measurement, this value defaults to 300 km following + Holland et al. 2010. + vmax_in_brackets : bool, optional + If True, use the alternative formula in equation (6) to solve for the peripheral exponent + x_n from the second measurement. Note that, a side-effect of the formula with vmax inside + of the brackets is that the wind speed maximum is attained a bit farther away from the + center than according to the recorded radius of maximum winds (RMW). Default: False Returns ------- - ndarray of shape (npositions, ncentroids) - containing the magnitude of the angular windspeed per track position per centroid location + hol_x : np.ndarray of shape (nnodes, ncentroids) + Exponents according to Holland et al. 2010. """ - _rho_air(si_track, rho_air_const) - _bs_holland_2008(si_track, gradient_to_surface_winds=gradient_to_surface_winds) - return _stat_holland_1980(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) + hol_x = np.zeros_like(d_centr) + r_max, v_max_s, hol_b, d_centr, v_n, r_n = [ + np.broadcast_to(ar, d_centr.shape)[mask_centr_close] + for ar in [ + si_track["rad"].values[:, None], + si_track["vmax"].values[:, None], + si_track["hol_b"].values[:, None], + d_centr, + np.atleast_1d(v_n)[:, None], + np.atleast_1d(r_n_km)[:, None], + ] + ] -def _compute_angular_windspeeds_h10( + # convert to SI units + r_n *= KM_TO_M + + # compute peripheral exponent from second measurement + # (equation (6) from Holland et al. 2010 solved for x) + r_max_norm = (r_max / r_n)**hol_b + if vmax_in_brackets: + x_n = np.log(v_n) / np.log(v_max_s**2 * r_max_norm * np.exp(1 - r_max_norm)) + + # With `vmax_in_brackets`, the maximum is shifted away from the recorded RMW. We truncate + # here to avoid an exaggerated shift. The value 1.0 has been found to be reasonable by + # manual testing of thresholds. Note that the truncation means that the peripheral wind + # speed v_n is not exactly attained in some cases. + x_n = np.fmin(x_n, 1.0) + else: + x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) + + # linearly interpolate between max exponent and peripheral exponent + x_max = 0.5 + hol_x[mask_centr_close] = x_max + np.fmax(0, d_centr - r_max) * (x_n - x_max) / (r_n - r_max) + + # Truncate to prevent wind speed from increasing again towards the peripheral radius (which is + # unphysical). A value of 0.4 has been found to be reasonable by manual testing of thresholds. + # Note that this means that the peripheral wind speed v_n is not exactly attained sometimes. + hol_x[mask_centr_close] = np.fmax(hol_x[mask_centr_close], 0.4) + + return hol_x + + +def _stat_holland_2010( si_track: xr.Dataset, d_centr: np.ndarray, - close_centr_msk: np.ndarray, - cyclostrophic: bool = True, - gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, - rho_air_const: float = DEF_RHO_AIR, - vmax_from_cen: bool = True, + mask_centr_close: np.ndarray, + hol_x: Union[float, np.ndarray], vmax_in_brackets: bool = False, -): - """Compute (absolute) angular wind speeds according to the Holland et al. 2010 model +) -> np.ndarray: + """Symmetric and static surface wind fields (in m/s) according to Holland et al. 2010 - Note that this model is always cyclostrophic, the parameter setting is ignored. + This function applies the cyclostrophic surface wind model expressed in equation (6) from + + Holland et al. (2010): A Revised Model for Radial Profiles of Hurricane Winds. Monthly + Weather Review 138(12): 4393–4401. https://doi.org/10.1175/2010MWR3317.1 + + More precisely, this function implements the following equation: + + V(r) = v_max_s * [(r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x + + In terms of this function's arguments, b_s is ``hol_b`` and r is ``d_centr``. + + If ``vmax_in_brackets`` is True, the alternative formula in (6) is used: + + .. math:: + + V(r) = [v_max_s^2 * (r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x Parameters ---------- si_track : xr.Dataset - Output of `tctrack_to_si`. Which data variables are used in the computation of the wind - profile depends on the selected model. - d_centr : np.ndarray of shape (npositions, ncentroids) - Distance (in m) between centroids and track positions. - close_centr_msk : np.ndarray of shape (npositions, ncentroids) - For each track position one row indicating which centroids are within reach. - cyclostrophic : bool, optional - This parameter is ignored because this model is always cyclostrophic. Default: True - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use. the wind profile is computed on the - surface level, but the clipping interval of the B-value depends on this factor. - Default: 0.9 - rho_air_const : float or None, optional - The constant value for air density (in kg/m³) to assume in the formula for the B-value. By - default, the value suggested in Holland 1980 is used. If set to None, the air density is - computed from pressure following equation (9) in Holland et al. 2010. Default: 1.15 - vmax_from_cen : boolean, optional - If True, replace the recorded value of vmax along the track by an estimate from pressure, - following equation (8) in Holland et al. 2010. Default: True + Output of ``tctrack_to_si`` with "hol_b" (see _bs_holland_2008) data variables. The data + variables used by this function are "vmax", "rad", and "hol_b". + d_centr : np.ndarray of shape (nnodes, ncentroids) + Distance (in m) between centroids and track nodes. + mask_centr_close : np.ndarray of shape (nnodes, ncentroids) + Mask indicating for each track node which centroids are within reach of the windfield. + hol_x : np.ndarray of shape (nnodes, ncentroids) or float + The exponent according to ``_x_holland_2010``. vmax_in_brackets : bool, optional - Specifies which of the two formulas in equation (6) of Holland et al. 2010 to use. If - False, the formula with vmax outside of the brackets is used. Note that, a side-effect of - the formula with vmax inside of the brackets is that the wind speed maximum is attained a - bit farther away from the center than according to the recorded radius of maximum + If True, use the alternative formula in equation (6). Note that, a side-effect of the + formula with vmax inside of the brackets is that the wind speed maximum is attained a bit + farther away from the center than according to the recorded radius of maximum winds (RMW). Default: False Returns ------- - ndarray of shape (npositions, ncentroids) - containing the magnitude of the angular windspeed per track position per centroid location + v_ang : np.ndarray (nnodes, ncentroids) + Absolute values of wind speeds (in m/s) in angular direction. """ - if not cyclostrophic: - LOGGER.debug( - 'The function _compute_angular_windspeeds_h10 was called with parameter ' - '"cyclostrophic" equal to false. Please be aware that this setting is ignored as the' - ' Holland et al. 2010 model is always cyclostrophic.') - _rho_air(si_track, rho_air_const) - if vmax_from_cen: - _bs_holland_2008(si_track, gradient_to_surface_winds=gradient_to_surface_winds) - _v_max_s_holland_2008(si_track) + v_ang = np.zeros_like(d_centr) + v_max_s, r_max, hol_b, d_centr, hol_x = [ + np.broadcast_to(ar, d_centr.shape)[mask_centr_close] + for ar in [ + si_track["vmax"].values[:, None], + si_track["rad"].values[:, None], + si_track["hol_b"].values[:, None], + d_centr, + hol_x, + ] + ] + + r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b + if vmax_in_brackets: + v_ang[mask_centr_close] = (v_max_s**2 * r_max_norm * np.exp(1 - r_max_norm))**hol_x else: - _B_holland_1980(si_track, gradient_to_surface_winds=gradient_to_surface_winds) - hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk, vmax_in_brackets=vmax_in_brackets) - return _stat_holland_2010( - si_track, d_centr, close_centr_msk, hol_x, vmax_in_brackets=vmax_in_brackets, - ) + v_ang[mask_centr_close] = v_max_s * (r_max_norm * np.exp(1 - r_max_norm))**hol_x + return v_ang -def get_close_centroids( - si_track: xr.Dataset, - centroids: np.ndarray, - buffer_km: float, - metric: str = "equirect", -) -> np.ndarray: - """Check whether centroids lay within a buffer around track positions - Note that, hypothetically, a problem occurs when a TC track is travelling so far in longitude - that adding a buffer exceeds 360 degrees (i.e. crosses the antimeridian), which is physically - impossible, but might happen with synthetical or test data. +def _stat_holland_1980( + si_track: xr.Dataset, + d_centr: np.ndarray, + mask_centr_close: np.ndarray, + cyclostrophic: bool = False +) -> np.ndarray: + """Symmetric and static wind fields (in m/s) according to Holland 1980. + + This function applies the gradient wind model expressed in equation (4) (combined with + equation (6)) from + + Holland, G.J. (1980): An Analytic Model of the Wind and Pressure Profiles in Hurricanes. + Monthly Weather Review 108(8): 1212–1218. + + More precisely, this function implements the following equation: + + V(r) = [(B/rho) * (r_max/r)^B * (penv - pcen) * e^(-(r_max/r)^B) + (r*f/2)^2]^0.5 - (r*f/2) + + In terms of this function's arguments, B is ``hol_b`` and r is ``d_centr``. The air density + rho and the Coriolis parameter f are taken from ``si_track``. + + Even though the equation has been derived originally for gradient winds (when combined with the + output of ``_B_holland_1980``), it can be used for surface winds by adjusting the parameter + ``hol_b`` (see function ``_bs_holland_2008``). Parameters ---------- - si_track : xr.Dataset with dimension "time" - Track information as returned by ``tctrack_to_si``. Hence, longitudinal coordinates are - normalized around the central longitude stored in the "mid_lon" attribute. This makes sure - that the buffered bounding box around the track does not cross the antimeridian. The data - variables used by this function are "lat", and "lon". - centroids : np.ndarray of shape (ncentroids, 2) - Coordinates of centroids, each row is a pair [lat, lon]. The longitudinal coordinates are - normalized within this function to be consistent with the track coordinates. - buffer_km : float - Size of the buffer (in km). The buffer is converted to a lat/lon buffer, rescaled in - longitudinal direction according to the t_lat coordinates. - metric : str, optional - Specify an approximation method to use for earth distances: "equirect" (faster) or - "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. - Default: "equirect". + si_track : xr.Dataset + Output of ``tctrack_to_si`` with "hol_b" (see, e.g., _B_holland_1980) and + "rho_air" (see _rho_air) data variable. The data variables used by this function + are "lat", "cp", "rad", "pdelta", "hol_b", and "rho_air". + d_centr : np.ndarray of shape (nnodes, ncentroids) + Distance (in m) between centroids and track nodes. + mask_centr_close : np.ndarray of shape (nnodes, ncentroids) + Mask indicating for each track node which centroids are within reach of the windfield. + cyclostrophic : bool, optional + If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0). + Default: False Returns ------- - centroids_close_normalized : np.ndarray of shape (nclose, 2) - Coordinates of close centroids, each row is a pair [lat, lon]. The normalization of - longitudinal coordinates is consistent with the track coordinates. - mask_centr : np.ndarray of shape (ncentroids,) - Mask that is True for close centroids and False for other centroids. - mask_centr_alongtrack : np.ndarray of shape (npositions, nclose) - Each row is a mask that indicates the centroids within reach for one track position. Note - that these masks refer only to the "close centroids" to reduce memory requirements. The - number of positions ``npositions`` corresponds to the size of the "time" dimension of - ``si_track``. + v_ang : np.ndarray (nnodes, ncentroids) + Absolute values of wind speeds (m/s) in angular direction. """ - npositions = si_track.sizes["time"] - ncentroids = centroids.shape[0] - t_lat, t_lon = si_track["lat"].values, si_track["lon"].values - centr_lat, centr_lon = centroids[:, 0].copy(), centroids[:, 1].copy() + v_ang = np.zeros_like(d_centr) + r_max, hol_b, pdelta, coriolis_p, rho_air, d_centr = [ + np.broadcast_to(ar, d_centr.shape)[mask_centr_close] + for ar in [ + si_track["rad"].values[:, None], + si_track["hol_b"].values[:, None], + si_track["pdelta"].values[:, None], + si_track["cp"].values[:, None], + si_track["rho_air"].values[:, None], + d_centr, + ] + ] - # Normalize longitudinal coordinates of centroids. - u_coord.lon_normalize(centr_lon, center=si_track.attrs["mid_lon"]) + r_coriolis = 0 + if not cyclostrophic: + r_coriolis = 0.5 * d_centr * coriolis_p - # Restrict to the bounding box of the whole track first (this can already reduce the number of - # centroids that are considered by a factor larger than 30). - buffer_lat = buffer_km / u_const.ONE_LAT_KM - buffer_lon = buffer_km / ( - u_const.ONE_LAT_KM * np.cos(np.radians( - np.fmin(89.999, np.abs(centr_lat) + buffer_lat) - )) - ) - [idx_close] = ( - (t_lat.min() - centr_lat <= buffer_lat) - & (centr_lat - t_lat.max() <= buffer_lat) - & (t_lon.min() - centr_lon <= buffer_lon) - & (centr_lon - t_lon.max() <= buffer_lon) - ).nonzero() - centr_lat = centr_lat[idx_close] - centr_lon = centr_lon[idx_close] + r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b + sqrt_term = hol_b / rho_air * r_max_norm * pdelta * np.exp(-r_max_norm) + r_coriolis**2 + v_ang[mask_centr_close] = np.sqrt(np.fmax(0, sqrt_term)) - r_coriolis + return v_ang - # Restrict to bounding boxes of each track position. - buffer_lat = buffer_km / u_const.ONE_LAT_KM - buffer_lon = buffer_km / (u_const.ONE_LAT_KM * np.cos(np.radians( - np.fmin(89.999, np.abs(t_lat[:, None]) + buffer_lat) - ))) - [idx_close_sub] = ( - (t_lat[:, None] - buffer_lat <= centr_lat[None]) - & (t_lat[:, None] + buffer_lat >= centr_lat[None]) - & (t_lon[:, None] - buffer_lon <= centr_lon[None]) - & (t_lon[:, None] + buffer_lon >= centr_lon[None]) - ).any(axis=0).nonzero() - idx_close = idx_close[idx_close_sub] - centr_lat = centr_lat[idx_close_sub] - centr_lon = centr_lon[idx_close_sub] - # Restrict to metric distance radius around each track position. - # - # We do the distance computation for chunks of the track since computing the distance requires - # npositions*ncentroids*8*3 Bytes of memory. For example, Hurricane FAITH's life time was more - # than 500 hours. At 0.5-hourly resolution and 1,000,000 centroids, that's 24 GB of memory for - # FAITH. With a chunk size of 10, this figure is down to 240 MB. The final along-track mask - # will require 1.0 GB of memory. - chunk_size = 10 - chunks = np.split(np.arange(npositions), np.arange(chunk_size, npositions, chunk_size)) - mask_centr_alongtrack = np.concatenate([ - ( - u_coord.dist_approx( - t_lat[None, chunk], t_lon[None, chunk], - centr_lat[None], centr_lon[None], - normalize=False, method=metric, units="km", - )[0] <= buffer_km - ) for chunk in chunks - ], axis=0) - [idx_close_sub] = mask_centr_alongtrack.any(axis=0).nonzero() - idx_close = idx_close[idx_close_sub] - centr_lat = centr_lat[idx_close_sub] - centr_lon = centr_lon[idx_close_sub] - mask_centr_alongtrack = mask_centr_alongtrack[:, idx_close_sub] +def _stat_er_2011( + si_track: xr.Dataset, + d_centr: np.ndarray, + mask_centr_close: np.ndarray, + cyclostrophic: bool = False, +) -> np.ndarray: + """Symmetric and static wind fields (in m/s) according to Emanuel and Rotunno 2011 - # Derive mask from index. - mask_centr = np.zeros((ncentroids,), dtype=bool) - mask_centr[idx_close] = True + Emanuel, K., Rotunno, R. (2011): Self-Stratification of Tropical Cyclone Outflow. Part I: + Implications for Storm Structure. Journal of the Atmospheric Sciences 68(10): 2236–2249. + https://dx.doi.org/10.1175/JAS-D-10-05024.1 + + The wind speeds ``v_ang`` are extracted from the momentum via the relationship M = v_ang * r, + where r corresponds to ``d_centr``. On the other hand, the momentum is derived from the momentum + at the peak wind position using equation (36) from Emanuel and Rotunno 2011 with Ck == Cd: + + M = M_max * [2 * (r / r_max)^2 / (1 + (r / r_max)^2)]. + + The momentum at the peak wind position is + + M_max = r_max * v_max + 0.5 * f * r_max**2, + + where the Coriolis parameter f is computed from the latitude ``lat`` using the constant rotation + rate of the earth. + + Parameters + ---------- + si_track : xr.Dataset + Output of ``tctrack_to_si``. The data variables used by this function are "lat", "cp", + "rad", and "vmax". + d_centr : np.ndarray of shape (nnodes, ncentroids) + Distance (in m) between centroids and track nodes. + mask_centr_close : np.ndarray of shape (nnodes, ncentroids) + Mask indicating for each track node which centroids are within reach of the windfield. + cyclostrophic : bool, optional + If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0) in + the computation of M_max. Default: False + + Returns + ------- + v_ang : np.ndarray (nnodes, ncentroids) + Absolute values of wind speeds (m/s) in angular direction. + """ + v_ang = np.zeros_like(d_centr) + r_max, v_max, coriolis_p, d_centr = [ + np.broadcast_to(ar, d_centr.shape)[mask_centr_close] + for ar in [ + si_track["rad"].values[:, None], + si_track["vmax"].values[:, None], + si_track["cp"].values[:, None], + d_centr, + ] + ] + + # compute the momentum at the maximum + momentum_max = r_max * v_max + + if not cyclostrophic: + # add the influence of the Coriolis force + momentum_max += 0.5 * coriolis_p * r_max**2 + + # rescale the momentum using formula (36) in Emanuel and Rotunno 2011 with Ck == Cd + r_max_norm = (d_centr / r_max)**2 + momentum = momentum_max * 2 * r_max_norm / (1 + r_max_norm) + + # extract the velocity from the rescaled momentum through division by r + v_ang[mask_centr_close] = np.fmax(0, momentum / (d_centr + 1e-11)) + return v_ang + + +KN_TO_MS = (1.0 * ureg.knot).to(ureg.meter / ureg.second).magnitude +V_ANG_EARTH = 7.29e-5 - centroids_close_normalized = np.stack([centr_lat, centr_lon], axis=1) - return centroids_close_normalized, mask_centr, mask_centr_alongtrack def _vtrans(si_track: xr.Dataset, metric: str = "equirect"): """Translational vector and velocity (in m/s) at each track node. @@ -1517,6 +844,7 @@ def _vtrans(si_track: xr.Dataset, metric: str = "equirect"): si_track["vtrans"].values[msk, :] *= fact[:, None] si_track["vtrans_norm"].values[msk] *= fact + def _coriolis_parameter(lat: np.ndarray) -> np.ndarray: """Compute the Coriolis parameter from latitude. @@ -1532,500 +860,534 @@ def _coriolis_parameter(lat: np.ndarray) -> np.ndarray: """ return 2 * V_ANG_EARTH * np.sin(np.radians(np.abs(lat))) -def _rho_air(si_track: xr.Dataset, const: Optional[float]): - """Eyewall density of air (in kg/m³) at each track node. - - The track dataset is modified in place, with the "rho_air" data variable added. +def compute_windfields_sparse( + track: xr.Dataset, + centroids: Centroids, + idx_centr_filter: np.ndarray, + model: str = 'H08', + model_kwargs: Optional[dict] = None, + store_windfields: bool = False, + metric: str = "equirect", + intensity_thres: float = DEF_INTENSITY_THRES, + max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, + max_memory_gb: float = DEF_MAX_MEMORY_GB, +) -> Tuple[sparse.csr_matrix, Optional[sparse.csr_matrix]]: + """Version of ``compute_windfields`` that returns sparse matrices and limits memory usage Parameters ---------- - si_track : xr.Dataset - Track information as returned by `tctrack_to_si`. The data variables used by this function - are "lat", "cen", and "pdelta". The result is stored in place as new data - variable "rho_air". - const : float or None - A constant value for air density (in kg/m³) to assume. If None, the air density is - estimated from eyewall pressure following equation (9) in Holland et al. 2010. - """ - if const is not None: - si_track["rho_air"] = xr.full_like(si_track["time"], const, dtype=float) - return - - # surface relative humidity (unitless), assumed constant following Holland et al. 2010 - surface_relative_humidity = 0.9 - - # surface temperature (in °C), following equation (9) in Holland 2008 - temp_s = 28.0 - 3.0 * (si_track["lat"] - 10.0) / 20.0 - - # eyewall surface pressure (in Pa), following equation (6) in Holland 2008 - pres_eyewall = si_track["cen"] + si_track["pdelta"] / np.exp(1) - - # mixing ratio (in kg/kg), estimated from temperature, using formula for saturation vapor - # pressure in Bolton 1980 (multiplied by the ratio of molar masses of water vapor and dry air) - # We multiply by 100, since the formula by Bolton is in hPa (mbar), and we use Pa. - r_mix = 100 * 3.802 / pres_eyewall * np.exp(17.67 * temp_s / (243.5 + temp_s)) - - # virtual surface temperature (in K) - temp_vs = (T_ICE_K + temp_s) * (1 + 0.81 * surface_relative_humidity * r_mix) - - # specific gas constant of dry air (in J/kgK) - r_dry_air = 286.9 - - # density of air (in kg/m³); when checking the units, note that J/Pa = m³ - si_track["rho_air"] = pres_eyewall / (r_dry_air * temp_vs) - -def _bs_holland_2008( - si_track: xr.Dataset, - gradient_to_surface_winds: float = DEF_GRADIENT_TO_SURFACE_WINDS, -): - """Holland's 2008 b-value estimate for sustained surface winds. - (This is also one option of how to estimate the b-value in Holland 2010, - for the other option consult '_bs_holland_2010_v2'.) - - The result is stored in place as a new data variable "hol_b". - - Unlike the original 1980 formula (see ``_B_holland_1980``), this approach does not require any - wind speed measurements, but is based on the more reliable pressure information. - - The parameter applies to 1-minute sustained winds at 10 meters above ground. - It is taken from equation (11) in the following paper: - - Holland, G. (2008). A revised hurricane pressure-wind model. Monthly - Weather Review, 136(9), 3432–3445. https://doi.org/10.1175/2008MWR2395.1 - - For reference, it reads - - b_s = -4.4 * 1e-5 * (penv - pcen)^2 + 0.01 * (penv - pcen) - + 0.03 * (dp/dt) - 0.014 * |lat| + 0.15 * (v_trans)^hol_xx + 1.0 - - where ``dp/dt`` is the time derivative of central pressure and ``hol_xx`` is Holland's x - parameter: hol_xx = 0.6 * (1 - (penv - pcen) / 215) - - The equation for b_s has been fitted statistically using hurricane best track records for - central pressure and maximum wind. It therefore performs best in the North Atlantic. + track : xr.Dataset + Single tropical cyclone track. + centroids : Centroids + Centroids instance. + idx_centr_filter : np.ndarray + Indices of centroids to restrict to (e.g. sufficiently close to coast). + model : str, optional + Parametric wind field model, one of "H1980" (the prominent Holland 1980 model), + "H08" (Holland 1980 with b-value from Holland 2008), "H10" (Holland et al. 2010), or + "ER11" (Emanuel and Rotunno 2011). + Default: "H08". + model_kwargs: dict, optional + If given, forward these kwargs to the selected model. Default: None + store_windfields : boolean, optional + If True, store windfields. Default: False. + metric : str, optional + Specify an approximation method to use for earth distances: "equirect" (faster) or + "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. + Default: "equirect". + intensity_thres : float, optional + Wind speeds (in m/s) below this threshold are stored as 0. Default: 17.5 + max_dist_eye_km : float, optional + No wind speed calculation is done for centroids with a distance (in km) to the TC + center ("eye") larger than this parameter. Default: 300 + max_memory_gb : float, optional + To avoid memory issues, the computation is done for chunks of the track sequentially. + The chunk size is determined depending on the available memory (in GB). Default: 8 - Furthermore, b_s has been fitted under the assumption of a "cyclostrophic" wind field which - means that the influence from Coriolis forces is assumed to be small. This is reasonable close - to the radius of maximum wind where the Coriolis term (r*f/2) is small compared to the rest - (see ``_stat_holland_1980``). More precisely: At the radius of maximum wind speeds, the typical - order of the Coriolis term is 1 while wind speed is 50 (which changes away from the - radius of maximum winds and as the TC moves away from the equator). + Raises + ------ + ValueError - Parameters - ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si``. The data variables used by this function are "lat", "tstep", - "vtrans_norm", "cen", and "pdelta". The result is stored in place as a new data - variable "hol_b". - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use when determining the clipping - interval. Default: 0.9 + Returns + ------- + intensity : csr_matrix + Maximum wind speed in each centroid over the whole storm life time. + windfields : csr_matrix or None + If store_windfields is True, the full velocity vectors at each centroid and track position + are stored in a sparse matrix of shape (npositions, ncentroids * 2) that can be reshaped + to a full ndarray of shape (npositions, ncentroids, 2). + If store_windfields is False, None is returned. """ - # adjust pressure at previous track point - prev_cen = np.zeros_like(si_track["cen"].values) - prev_cen[1:] = si_track["cen"].values[:-1].copy() - msk = prev_cen < 850 * MBAR_TO_PA - prev_cen[msk] = si_track["cen"].values[msk] - - # The formula assumes that pressure values are in millibar (hPa) instead of SI units (Pa), - # and time steps are in hours instead of seconds, but translational wind speed is still - # expected to be in m/s. - pdelta = si_track["pdelta"] / MBAR_TO_PA - hol_xx = 0.6 * (1. - pdelta / 215) - si_track["hol_b"] = ( - -4.4e-5 * pdelta**2 + 0.01 * pdelta - + 0.03 * (si_track["cen"] - prev_cen) / si_track["tstep"] * (H_TO_S / MBAR_TO_PA) - - 0.014 * abs(si_track["lat"]) - + 0.15 * si_track["vtrans_norm"]**hol_xx + 1.0 - ) - clip_interval = _b_holland_clip_interval(gradient_to_surface_winds) - si_track["hol_b"] = np.clip(si_track["hol_b"], *clip_interval) - -def _v_max_s_holland_2008(si_track: xr.Dataset): - """Compute maximum surface winds from pressure according to Holland 2008. - - The result is stored in place as a data variable "vmax". If a variable of that name already - exists, its values are overwritten. - - This function implements equation (11) in the following paper: - - Holland, G. (2008). A revised hurricane pressure-wind model. Monthly - Weather Review, 136(9), 3432–3445. https://doi.org/10.1175/2008MWR2395.1 - - For reference, it reads - - v_ms = [b_s / (rho * e) * (penv - pcen)]^0.5 + try: + mod_id = MODEL_VANG[model] + except KeyError as err: + raise ValueError(f'Model not implemented: {model}.') from err - where ``b_s`` is Holland b-value (see ``_bs_holland_2008``), e is Euler's number, rho is the - density of air, ``penv`` is environmental, and ``pcen`` is central pressure. + ncentroids = centroids.coord.shape[0] + npositions = track.sizes["time"] + windfields_shape = (npositions, ncentroids * 2) + intensity_shape = (1, ncentroids) - Parameters - ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si`` with "hol_b" (see _bs_holland_2008) and "rho_air" (see - _rho_air) variables. The data variables used by this function are "pdelta", "hol_b", - and "rho_air". The results are stored in place as a new data variable "vmax". If a variable - of that name already exists, its values are overwritten. - """ - si_track["vmax"] = np.sqrt( - si_track["hol_b"] / (si_track["rho_air"] * np.exp(1)) * si_track["pdelta"] + # initialise arrays for the assumption that no centroids are within reach + windfields_sparse = ( + sparse.csr_matrix(([], ([], [])), shape=windfields_shape) + if store_windfields else None ) + intensity_sparse = sparse.csr_matrix(([], ([], [])), shape=intensity_shape) -def _B_holland_1980( # pylint: disable=invalid-name - si_track: xr.Dataset, - gradient_to_surface_winds: Optional[float] = None, -): - """Holland's 1980 B-value computation for gradient-level winds. - - The result is stored in place as a new data variable "hol_b". + # The wind field model requires at least two track positions because translational speed + # as well as the change in pressure (in case of H08) are required. + if npositions < 2: + return intensity_sparse, windfields_sparse - The parameter applies to gradient-level winds (about 1000 metres above the earth's surface). - The formula for B is derived from equations (5) and (6) in the following paper: + # convert track variables to SI units + si_track = tctrack_to_si(track, metric=metric) - Holland, G.J. (1980): An Analytic Model of the Wind and Pressure Profiles - in Hurricanes. Monthly Weather Review 108(8): 1212–1218. - https://doi.org/10.1175/1520-0493(1980)108<1212:AAMOTW>2.0.CO;2 + # When done properly, finding and storing the close centroids is not a memory bottle neck and + # can be done before chunking. Note that the longitudinal coordinates of `centroids_close` as + # returned by `get_close_centroids` are normalized to be consistent with the coordinates in + # `si_track`. + centroids_close, mask_centr, mask_centr_alongtrack = get_close_centroids( + si_track, centroids.coord[idx_centr_filter], max_dist_eye_km, metric=metric, + ) + idx_centr_filter = idx_centr_filter[mask_centr] + n_centr_close = centroids_close.shape[0] + if n_centr_close == 0: + return intensity_sparse, windfields_sparse - For reference, inserting (6) into (5) and solving for B at r = RMW yields: + # the total memory requirement in GB if we compute everything without chunking: + # 8 Bytes per entry (float64), 10 arrays + total_memory_gb = npositions * n_centr_close * 8 * 10 / 1e9 + if total_memory_gb > max_memory_gb and npositions > 2: + # If the number of positions is down to 2 already, we cannot split any further. In that + # case, we just take the risk and try to do the computation anyway. It might still work + # since we have only computed an upper bound for the number of affected centroids. - B = v^2 * e * rho / (penv - pcen) + # Split the track into chunks, compute the result for each chunk, and combine: + return _compute_windfields_sparse_chunked( + mask_centr_alongtrack, + track, + centroids, + idx_centr_filter, + model=model, + model_kwargs=model_kwargs, + store_windfields=store_windfields, + metric=metric, + intensity_thres=intensity_thres, + max_dist_eye_km=max_dist_eye_km, + max_memory_gb=max_memory_gb, + ) - where v are maximum gradient-level winds ``gradient_winds``, e is Euler's number, rho is the - density of air, ``penv`` is environmental, and ``pcen`` is central pressure. + windfields, idx_centr_reachable = _compute_windfields( + si_track, + centroids_close, + mod_id, + model_kwargs=model_kwargs, + metric=metric, + max_dist_eye_km=max_dist_eye_km, + ) + idx_centr_filter = idx_centr_filter[idx_centr_reachable] + npositions = windfields.shape[0] - Parameters - ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si`` with "rho_air" variable (see _rho_air). The data variables - used by this function are "vgrad" (or "vmax" if gradient_to_surface_winds is different from - 1.0), "pdelta", and "rho_air". The results are stored in place as a new data - variable "hol_b". - gradient_to_surface_winds : float, optional - The gradient-to-surface wind reduction factor to use when determining the clipping - interval. By default, the gradient level values are assumed. Default: None - """ - windvar = "vgrad" if gradient_to_surface_winds is None else "vmax" + intensity = np.linalg.norm(windfields, axis=-1).max(axis=0) + intensity[intensity < intensity_thres] = 0 + intensity_sparse = sparse.csr_matrix( + (intensity, idx_centr_filter, [0, intensity.size]), + shape=intensity_shape) + intensity_sparse.eliminate_zeros() - si_track["hol_b"] = ( - si_track[windvar]**2 * np.exp(1) * si_track["rho_air"] / si_track["pdelta"] - ) + windfields_sparse = None + if store_windfields: + n_centr_filter = idx_centr_filter.size + indices = np.zeros((npositions, n_centr_filter, 2), dtype=np.int64) + indices[:, :, 0] = 2 * idx_centr_filter[None] + indices[:, :, 1] = 2 * idx_centr_filter[None] + 1 + indices = indices.ravel() + indptr = np.arange(npositions + 1) * n_centr_filter * 2 + windfields_sparse = sparse.csr_matrix((windfields.ravel(), indices, indptr), + shape=windfields_shape) + windfields_sparse.eliminate_zeros() - clip_interval = _b_holland_clip_interval(gradient_to_surface_winds) - si_track["hol_b"] = np.clip(si_track["hol_b"], *clip_interval) + return intensity_sparse, windfields_sparse -def _b_holland_clip_interval(gradient_to_surface_winds): - """The clip interval to use for the Holland B-value - The default clip interval for gradient level B-values is taken to be (1.0, 2.5), following - Holland 1980. +def _compute_windfields_sparse_chunked( + mask_centr_alongtrack: np.ndarray, + track: xr.Dataset, + *args, + max_memory_gb: float = DEF_MAX_MEMORY_GB, + **kwargs, +) -> Tuple[sparse.csr_matrix, Optional[sparse.csr_matrix]]: + """Call ``compute_windfields_sparse`` for chunks of the track and re-assemble the results Parameters ---------- - gradient_to_surface_winds : float or None - The gradient-to-surface wind reduction factor to use when rescaling the gradient-level - clip interval (1.0, 2.5) proposed in Holland 1980. If None, no rescaling is applied. + mask_centr_alongtrack : np.ndarray of shape (npositions, ncentroids) + Each row is a mask that indicates the centroids within reach for one track position. + track : xr.Dataset + Single tropical cyclone track. + max_memory_gb : float, optional + Maximum memory requirements (in GB) for the computation of a single chunk of the track. + Default: 8 + args, kwargs : + The remaining arguments are passed on to ``compute_windfields_sparse``. Returns ------- - b_min, b_max : float - Minimum and maximum value of the clip interval. + intensity, windfields : + See ``compute_windfields_sparse`` for a description of the return values. """ - clip_interval = (1.0, 2.5) - if gradient_to_surface_winds is not None: - fact = gradient_to_surface_winds**2 - clip_interval = tuple(c * fact for c in clip_interval) - return clip_interval - -def _x_holland_2010( - si_track: xr.Dataset, - d_centr: np.ndarray, - mask_centr_close: np.ndarray, - v_n: Union[float, np.ndarray] = 17.0, - r_n_km: Union[float, np.ndarray] = 300.0, - vmax_in_brackets: bool = False, -) -> np.ndarray: - """Compute exponent for wind model according to Holland et al. 2010. + npositions = track.sizes["time"] + # The memory requirements for each track position are estimated for the case of 10 arrays + # containing `nreachable` float64 (8 Byte) values each. The chunking is only relevant in + # extreme cases with a very high temporal and/or spatial resolution. + max_nreachable = max_memory_gb * 1e9 / (8 * 10 * npositions) + split_pos = [0] + chunk_size = 2 + while split_pos[-1] + chunk_size < npositions: + chunk_size += 1 + # create overlap between consecutive chunks + chunk_start = max(0, split_pos[-1] - 1) + chunk_end = chunk_start + chunk_size + nreachable = mask_centr_alongtrack[chunk_start:chunk_end].any(axis=0).sum() + if nreachable > max_nreachable: + split_pos.append(chunk_end - 1) + chunk_size = 2 + split_pos.append(npositions) - This function implements equation (10) from the following paper: + intensity = [] + windfields = [] + for prev_chunk_end, chunk_end in zip(split_pos[:-1], split_pos[1:]): + chunk_start = max(0, prev_chunk_end - 1) + inten, win = compute_windfields_sparse( + track.isel(time=slice(chunk_start, chunk_end)), *args, + max_memory_gb=max_memory_gb, **kwargs, + ) + intensity.append(inten) + windfields.append(win) - Holland et al. (2010): A Revised Model for Radial Profiles of Hurricane Winds. Monthly - Weather Review 138(12): 4393–4401. https://doi.org/10.1175/2010MWR3317.1 + intensity = sparse.csr_matrix(sparse.vstack(intensity).max(axis=0)) + if windfields[0] is not None: + # eliminate the overlap between consecutive chunks + windfields = [windfields[0]] + [win[1:, :] for win in windfields[1:]] + windfields = sparse.vstack(windfields, format="csr") + return intensity, windfields - For reference, it reads - x = 0.5 [for r < r_max] - x = 0.5 + (r - r_max) * (x_n - 0.5) / (r_n - r_max) [for r >= r_max] +def _compute_windfields( + si_track: xr.Dataset, + centroids: np.ndarray, + model: int, + model_kwargs: Optional[dict] = None, + metric: str = "equirect", + max_dist_eye_km: float = DEF_MAX_DIST_EYE_KM, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute 1-minute sustained winds (in m/s) at 10 meters above ground - The peripheral exponent x_n is adjusted to fit the peripheral observation of wind speeds ``v_n`` - at radius ``r_n``. + In a first step, centroids within reach of the track are determined so that wind fields will + only be computed and returned for those centroids. Still, since computing the distance of + the storm center to the centroids is computationally expensive, make sure to pre-filter the + centroids and call this function only for those centroids that are potentially affected. Parameters ---------- si_track : xr.Dataset - Output of ``tctrack_to_si`` with "hol_b" variable (see _bs_holland_2008). The data variables - used by this function are "rad", "vmax", and "hol_b". - d_centr : np.ndarray of shape (nnodes, ncentroids) - Distance (in m) between centroids and track nodes. - mask_centr_close : np.ndarray of shape (nnodes, ncentroids) - Mask indicating for each track node which centroids are within reach of the windfield. - v_n : np.ndarray of shape (nnodes,) or float, optional - Peripheral wind speeds (in m/s) at radius ``r_n`` outside of radius of maximum winds - ``r_max``. In absence of a second wind speed measurement, this value defaults to 17 m/s - following Holland et al. 2010 (at a radius of 300 km). - r_n_km : np.ndarray of shape (nnodes,) or float, optional - Radius (in km) where the peripheral wind speed ``v_n`` is measured (or assumed). - In absence of a second wind speed measurement, this value defaults to 300 km following - Holland et al. 2010. - vmax_in_brackets : bool, optional - If True, use the alternative formula in equation (6) to solve for the peripheral exponent - x_n from the second measurement. Note that, a side-effect of the formula with vmax inside - of the brackets is that the wind speed maximum is attained a bit farther away from the - center than according to the recorded radius of maximum winds (RMW). Default: False + Output of ``tctrack_to_si``. Which data variables are used in the computation of the wind + speeds depends on the selected model. + centroids : np.ndarray with two dimensions + Each row is a centroid [lat, lon]. Centroids that are not within reach of the track are + ignored. Longitudinal coordinates are assumed to be normalized consistently with the + longitudinal coordinates in ``si_track``. + model : int + Wind profile model selection according to MODEL_VANG. + model_kwargs: dict, optional + If given, forward these kwargs to the selected model. Default: None + metric : str, optional + Specify an approximation method to use for earth distances: "equirect" (faster) or + "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. + Default: "equirect". + max_dist_eye_km : float, optional + No wind speed calculation is done for centroids with a distance (in km) to the TC center + ("eye") larger than this parameter. Default: 300 Returns ------- - hol_x : np.ndarray of shape (nnodes, ncentroids) - Exponents according to Holland et al. 2010. + windfields : np.ndarray of shape (npositions, nreachable, 2) + Directional wind fields for each track position on those centroids within reach + of the TC track. Note that the wind speeds at the first position are all zero because + the discrete time derivatives involved in the process are implemented using backward + differences. However, the first position is usually not relevant for impact calculations + since it is far off shore. + idx_centr_reachable : np.ndarray of shape (nreachable,) + List of indices of input centroids within reach of the TC track. """ - hol_x = np.zeros_like(d_centr) - r_max, v_max_s, hol_b, d_centr, v_n, r_n = [ - np.broadcast_to(ar, d_centr.shape)[mask_centr_close] - for ar in [ - si_track["rad"].values[:, None], - si_track["vmax"].values[:, None], - si_track["hol_b"].values[:, None], - d_centr, - np.atleast_1d(v_n)[:, None], - np.atleast_1d(r_n_km)[:, None], - ] - ] - - # convert to SI units - r_n *= KM_TO_M - - # compute peripheral exponent from second measurement - # (equation (6) from Holland et al. 2010 solved for x) - r_max_norm = (r_max / r_n)**hol_b - if vmax_in_brackets: - x_n = np.log(v_n) / np.log(v_max_s**2 * r_max_norm * np.exp(1 - r_max_norm)) - - # With `vmax_in_brackets`, the maximum is shifted away from the recorded RMW. We truncate - # here to avoid an exaggerated shift. The value 1.0 has been found to be reasonable by - # manual testing of thresholds. Note that the truncation means that the peripheral wind - # speed v_n is not exactly attained in some cases. - x_n = np.fmin(x_n, 1.0) - else: - x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) - - # linearly interpolate between max exponent and peripheral exponent - x_max = 0.5 - hol_x[mask_centr_close] = x_max + np.fmax(0, d_centr - r_max) * (x_n - x_max) / (r_n - r_max) - - # Truncate to prevent wind speed from increasing again towards the peripheral radius (which is - # unphysical). A value of 0.4 has been found to be reasonable by manual testing of thresholds. - # Note that this means that the peripheral wind speed v_n is not exactly attained sometimes. - hol_x[mask_centr_close] = np.fmax(hol_x[mask_centr_close], 0.4) - - return hol_x - -def _stat_holland_2010( - si_track: xr.Dataset, - d_centr: np.ndarray, - mask_centr_close: np.ndarray, - hol_x: Union[float, np.ndarray], - vmax_in_brackets: bool = False, -) -> np.ndarray: - """Symmetric and static surface wind fields (in m/s) according to Holland et al. 2010 - - This function applies the cyclostrophic surface wind model expressed in equation (6) from - - Holland et al. (2010): A Revised Model for Radial Profiles of Hurricane Winds. Monthly - Weather Review 138(12): 4393–4401. https://doi.org/10.1175/2010MWR3317.1 + # start with the assumption that no centroids are within reach + npositions = si_track.sizes["time"] + idx_centr_reachable = np.zeros((0,), dtype=np.int64) + windfields = np.zeros((npositions, 0, 2), dtype=np.float64) - More precisely, this function implements the following equation: + # compute distances (in m) and vectors to all centroids + [d_centr], [v_centr_normed] = u_coord.dist_approx( + si_track["lat"].values[None], si_track["lon"].values[None], + centroids[None, :, 0], centroids[None, :, 1], + log=True, normalize=False, method=metric, units="m") - V(r) = v_max_s * [(r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x + # exclude centroids that are too far from or too close to the eye + mask_centr_close = (d_centr <= max_dist_eye_km * KM_TO_M) & (d_centr > 1) + if not np.any(mask_centr_close): + return windfields, idx_centr_reachable - In terms of this function's arguments, b_s is ``hol_b`` and r is ``d_centr``. + # restrict to the centroids that are within reach of any of the positions + mask_centr_close_any = mask_centr_close.any(axis=0) + mask_centr_close = mask_centr_close[:, mask_centr_close_any] + d_centr = d_centr[:, mask_centr_close_any] + v_centr_normed = v_centr_normed[:, mask_centr_close_any, :] - If ``vmax_in_brackets`` is True, the alternative formula in (6) is used: + # normalize the vectors pointing from the eye to the centroids + v_centr_normed[~mask_centr_close] = 0 + v_centr_normed[mask_centr_close] /= d_centr[mask_centr_close, None] - .. math:: + # derive (absolute) angular velocity from parametric wind profile + v_ang_norm = compute_angular_windspeeds( + si_track, d_centr, mask_centr_close, model, model_kwargs=model_kwargs, cyclostrophic=False, + ) - V(r) = [v_max_s^2 * (r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x + # Influence of translational speed decreases with distance from eye. + # The "absorbing factor" is according to the following paper (see Fig. 7): + # + # Mouton, F. & Nordbeck, O. (2005). Cyclone Database Manager. A tool + # for converting point data from cyclone observations into tracks and + # wind speed profiles in a GIS. UNED/GRID-Geneva. + # https://unepgrid.ch/en/resource/19B7D302 + # + t_rad_bc = np.broadcast_to(si_track["rad"].values[:, None], d_centr.shape) + v_trans_corr = np.zeros_like(d_centr) + v_trans_corr[mask_centr_close] = np.fmin( + 1, t_rad_bc[mask_centr_close] / d_centr[mask_centr_close]) - Parameters - ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si`` with "hol_b" (see _bs_holland_2008) data variables. The data - variables used by this function are "vmax", "rad", and "hol_b". - d_centr : np.ndarray of shape (nnodes, ncentroids) - Distance (in m) between centroids and track nodes. - mask_centr_close : np.ndarray of shape (nnodes, ncentroids) - Mask indicating for each track node which centroids are within reach of the windfield. - hol_x : np.ndarray of shape (nnodes, ncentroids) or float - The exponent according to ``_x_holland_2010``. - vmax_in_brackets : bool, optional - If True, use the alternative formula in equation (6). Note that, a side-effect of the - formula with vmax inside of the brackets is that the wind speed maximum is attained a bit - farther away from the center than according to the recorded radius of maximum - winds (RMW). Default: False + if model in [MODEL_VANG['H08'], MODEL_VANG['H10']]: + # In these models, v_ang_norm already contains vtrans_norm, so subtract it first, before + # converting to vectors and then adding (vectorial) vtrans again. Make sure to apply the + # "absorbing factor" in both steps: + vtrans_norm_bc = np.broadcast_to(si_track["vtrans_norm"].values[:, None], d_centr.shape) + v_ang_norm[mask_centr_close] -= ( + vtrans_norm_bc[mask_centr_close] * v_trans_corr[mask_centr_close] + ) - Returns - ------- - v_ang : np.ndarray (nnodes, ncentroids) - Absolute values of wind speeds (in m/s) in angular direction. - """ - v_ang = np.zeros_like(d_centr) - v_max_s, r_max, hol_b, d_centr, hol_x = [ - np.broadcast_to(ar, d_centr.shape)[mask_centr_close] - for ar in [ - si_track["vmax"].values[:, None], - si_track["rad"].values[:, None], - si_track["hol_b"].values[:, None], - d_centr, - hol_x, - ] - ] + # vectorial angular velocity + windfields = ( + si_track.attrs["latsign"] * np.array([1.0, -1.0])[..., :] * v_centr_normed[:, :, ::-1] + ) + windfields[mask_centr_close] *= v_ang_norm[mask_centr_close, None] - r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b - if vmax_in_brackets: - v_ang[mask_centr_close] = (v_max_s**2 * r_max_norm * np.exp(1 - r_max_norm))**hol_x - else: - v_ang[mask_centr_close] = v_max_s * (r_max_norm * np.exp(1 - r_max_norm))**hol_x - return v_ang + # add angular and corrected translational velocity vectors + windfields[1:] += si_track["vtrans"].values[1:, None, :] * v_trans_corr[1:, :, None] + windfields[np.isnan(windfields)] = 0 + windfields[0, :, :] = 0 + [idx_centr_reachable] = mask_centr_close_any.nonzero() + return windfields, idx_centr_reachable -def _stat_holland_1980( - si_track: xr.Dataset, - d_centr: np.ndarray, - mask_centr_close: np.ndarray, - cyclostrophic: bool = False -) -> np.ndarray: - """Symmetric and static wind fields (in m/s) according to Holland 1980. - This function applies the gradient wind model expressed in equation (4) (combined with - equation (6)) from +def tctrack_to_si( + track: xr.Dataset, + metric: str = "equirect", +) -> xr.Dataset: + """Convert track variables to SI units and prepare for wind field computation - Holland, G.J. (1980): An Analytic Model of the Wind and Pressure Profiles in Hurricanes. - Monthly Weather Review 108(8): 1212–1218. + In addition to unit conversion, the variable names are shortened, the longitudinal coordinates + are normalized and additional variables are defined: - More precisely, this function implements the following equation: + * cp (coriolis parameter) + * pdelta (difference between environmental and central pressure, always strictly positive) + * vtrans (translational velocity vectors) + * vtrans_norm (absolute value of translational speed) - V(r) = [(B/rho) * (r_max/r)^B * (penv - pcen) * e^(-(r_max/r)^B) + (r*f/2)^2]^0.5 - (r*f/2) + Furthermore, some scalar variables are stored as attributes: - In terms of this function's arguments, B is ``hol_b`` and r is ``d_centr``. The air density - rho and the Coriolis parameter f are taken from ``si_track``. + * latsign (1.0 if the track is located on the northern and -1.0 if on southern hemisphere) + * mid_lon (the central longitude that was used to normalize the longitudinal coordinates) - Even though the equation has been derived originally for gradient winds (when combined with the - output of ``_B_holland_1980``), it can be used for surface winds by adjusting the parameter - ``hol_b`` (see function ``_bs_holland_2008``). + Finally, some corrections are applied to variables: + + * clip central pressure values so that environmental pressure values are never exceeded + * extrapolate radius of max wind from pressure if missing Parameters ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si`` with "hol_b" (see, e.g., _B_holland_1980) and - "rho_air" (see _rho_air) data variable. The data variables used by this function - are "lat", "cp", "rad", "pdelta", "hol_b", and "rho_air". - d_centr : np.ndarray of shape (nnodes, ncentroids) - Distance (in m) between centroids and track nodes. - mask_centr_close : np.ndarray of shape (nnodes, ncentroids) - Mask indicating for each track node which centroids are within reach of the windfield. - cyclostrophic : bool, optional - If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0). - Default: False + track : xr.Dataset + Track information. + metric : str, optional + Specify an approximation method to use for earth distances: "equirect" (faster) or + "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. + Default: "equirect". Returns ------- - v_ang : np.ndarray (nnodes, ncentroids) - Absolute values of wind speeds (m/s) in angular direction. + xr.Dataset """ - v_ang = np.zeros_like(d_centr) - r_max, hol_b, pdelta, coriolis_p, rho_air, d_centr = [ - np.broadcast_to(ar, d_centr.shape)[mask_centr_close] - for ar in [ - si_track["rad"].values[:, None], - si_track["hol_b"].values[:, None], - si_track["pdelta"].values[:, None], - si_track["cp"].values[:, None], - si_track["rho_air"].values[:, None], - d_centr, - ] + si_track = track[["lat", "lon", "time"]].copy() + si_track["tstep"] = track["time_step"] * H_TO_S + si_track["env"] = track["environmental_pressure"] * MBAR_TO_PA + + # we support some non-standard unit names + unit_replace = {"mb": "mbar", "kn": "knots"} + configs = [ + ("central_pressure", "cen", "Pa"), + ("max_sustained_wind", "vmax", "m/s"), ] + for long_name, var_name, si_unit in configs: + unit = track.attrs[f"{long_name}_unit"] + unit = unit_replace.get(unit, unit) + try: + conv_factor = ureg(unit).to(si_unit).magnitude + except Exception as ex: + raise ValueError( + f"The {long_name}_unit '{unit}' in the provided track is not supported." + ) from ex + si_track[var_name] = track[long_name] * conv_factor - r_coriolis = 0 - if not cyclostrophic: - r_coriolis = 0.5 * d_centr * coriolis_p + # normalize longitudinal coordinates + si_track.attrs["mid_lon"] = 0.5 * sum(u_coord.lon_bounds(si_track["lon"].values)) + u_coord.lon_normalize(si_track["lon"].values, center=si_track.attrs["mid_lon"]) - r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b - sqrt_term = hol_b / rho_air * r_max_norm * pdelta * np.exp(-r_max_norm) + r_coriolis**2 - v_ang[mask_centr_close] = np.sqrt(np.fmax(0, sqrt_term)) - r_coriolis - return v_ang + # make sure that central pressure never exceeds environmental pressure + pres_exceed_msk = (si_track["cen"] > si_track["env"]).values + si_track["cen"].values[pres_exceed_msk] = si_track["env"].values[pres_exceed_msk] -def _stat_er_2011( - si_track: xr.Dataset, - d_centr: np.ndarray, - mask_centr_close: np.ndarray, - cyclostrophic: bool = False, -) -> np.ndarray: - """Symmetric and static wind fields (in m/s) according to Emanuel and Rotunno 2011 + # extrapolate radius of max wind from pressure if not given + si_track["rad"] = track["radius_max_wind"].copy() + si_track["rad"].values[:] = estimate_rmw( + si_track["rad"].values, si_track["cen"].values / MBAR_TO_PA, + ) + si_track["rad"] *= NM_TO_KM * KM_TO_M - Emanuel, K., Rotunno, R. (2011): Self-Stratification of Tropical Cyclone Outflow. Part I: - Implications for Storm Structure. Journal of the Atmospheric Sciences 68(10): 2236–2249. - https://dx.doi.org/10.1175/JAS-D-10-05024.1 + hemisphere = 'N' + if np.count_nonzero(si_track["lat"] < 0) > np.count_nonzero(si_track["lat"] > 0): + hemisphere = 'S' + si_track.attrs["latsign"] = 1.0 if hemisphere == 'N' else -1.0 - The wind speeds ``v_ang`` are extracted from the momentum via the relationship M = v_ang * r, - where r corresponds to ``d_centr``. On the other hand, the momentum is derived from the momentum - at the peak wind position using equation (36) from Emanuel and Rotunno 2011 with Ck == Cd: + # add translational speed of track at every node (in m/s) + _vtrans(si_track, metric=metric) - M = M_max * [2 * (r / r_max)^2 / (1 + (r / r_max)^2)]. + # add Coriolis parameter + si_track["cp"] = ("time", _coriolis_parameter(si_track["lat"].values)) - The momentum at the peak wind position is + # add pressure drop + si_track["pdelta"] = np.fmax(np.spacing(1), si_track["env"] - si_track["cen"]) - M_max = r_max * v_max + 0.5 * f * r_max**2, + return si_track - where the Coriolis parameter f is computed from the latitude ``lat`` using the constant rotation - rate of the earth. + +def get_close_centroids( + si_track: xr.Dataset, + centroids: np.ndarray, + buffer_km: float, + metric: str = "equirect", +) -> np.ndarray: + """Check whether centroids lay within a buffer around track positions + + Note that, hypothetically, a problem occurs when a TC track is travelling so far in longitude + that adding a buffer exceeds 360 degrees (i.e. crosses the antimeridian), which is physically + impossible, but might happen with synthetical or test data. Parameters ---------- - si_track : xr.Dataset - Output of ``tctrack_to_si``. The data variables used by this function are "lat", "cp", - "rad", and "vmax". - d_centr : np.ndarray of shape (nnodes, ncentroids) - Distance (in m) between centroids and track nodes. - mask_centr_close : np.ndarray of shape (nnodes, ncentroids) - Mask indicating for each track node which centroids are within reach of the windfield. - cyclostrophic : bool, optional - If True, do not apply the influence of the Coriolis force (set the Coriolis terms to 0) in - the computation of M_max. Default: False + si_track : xr.Dataset with dimension "time" + Track information as returned by ``tctrack_to_si``. Hence, longitudinal coordinates are + normalized around the central longitude stored in the "mid_lon" attribute. This makes sure + that the buffered bounding box around the track does not cross the antimeridian. The data + variables used by this function are "lat", and "lon". + centroids : np.ndarray of shape (ncentroids, 2) + Coordinates of centroids, each row is a pair [lat, lon]. The longitudinal coordinates are + normalized within this function to be consistent with the track coordinates. + buffer_km : float + Size of the buffer (in km). The buffer is converted to a lat/lon buffer, rescaled in + longitudinal direction according to the t_lat coordinates. + metric : str, optional + Specify an approximation method to use for earth distances: "equirect" (faster) or + "geosphere" (more accurate). See ``dist_approx`` function in ``climada.util.coordinates``. + Default: "equirect". Returns ------- - v_ang : np.ndarray (nnodes, ncentroids) - Absolute values of wind speeds (m/s) in angular direction. + centroids_close_normalized : np.ndarray of shape (nclose, 2) + Coordinates of close centroids, each row is a pair [lat, lon]. The normalization of + longitudinal coordinates is consistent with the track coordinates. + mask_centr : np.ndarray of shape (ncentroids,) + Mask that is True for close centroids and False for other centroids. + mask_centr_alongtrack : np.ndarray of shape (npositions, nclose) + Each row is a mask that indicates the centroids within reach for one track position. Note + that these masks refer only to the "close centroids" to reduce memory requirements. The + number of positions ``npositions`` corresponds to the size of the "time" dimension of + ``si_track``. """ - v_ang = np.zeros_like(d_centr) - r_max, v_max, coriolis_p, d_centr = [ - np.broadcast_to(ar, d_centr.shape)[mask_centr_close] - for ar in [ - si_track["rad"].values[:, None], - si_track["vmax"].values[:, None], - si_track["cp"].values[:, None], - d_centr, - ] - ] + npositions = si_track.sizes["time"] + ncentroids = centroids.shape[0] + t_lat, t_lon = si_track["lat"].values, si_track["lon"].values + centr_lat, centr_lon = centroids[:, 0].copy(), centroids[:, 1].copy() - # compute the momentum at the maximum - momentum_max = r_max * v_max + # Normalize longitudinal coordinates of centroids. + u_coord.lon_normalize(centr_lon, center=si_track.attrs["mid_lon"]) - if not cyclostrophic: - # add the influence of the Coriolis force - momentum_max += 0.5 * coriolis_p * r_max**2 + # Restrict to the bounding box of the whole track first (this can already reduce the number of + # centroids that are considered by a factor larger than 30). + buffer_lat = buffer_km / u_const.ONE_LAT_KM + buffer_lon = buffer_km / ( + u_const.ONE_LAT_KM * np.cos(np.radians( + np.fmin(89.999, np.abs(centr_lat) + buffer_lat) + )) + ) + [idx_close] = ( + (t_lat.min() - centr_lat <= buffer_lat) + & (centr_lat - t_lat.max() <= buffer_lat) + & (t_lon.min() - centr_lon <= buffer_lon) + & (centr_lon - t_lon.max() <= buffer_lon) + ).nonzero() + centr_lat = centr_lat[idx_close] + centr_lon = centr_lon[idx_close] - # rescale the momentum using formula (36) in Emanuel and Rotunno 2011 with Ck == Cd - r_max_norm = (d_centr / r_max)**2 - momentum = momentum_max * 2 * r_max_norm / (1 + r_max_norm) + # Restrict to bounding boxes of each track position. + buffer_lat = buffer_km / u_const.ONE_LAT_KM + buffer_lon = buffer_km / (u_const.ONE_LAT_KM * np.cos(np.radians( + np.fmin(89.999, np.abs(t_lat[:, None]) + buffer_lat) + ))) + [idx_close_sub] = ( + (t_lat[:, None] - buffer_lat <= centr_lat[None]) + & (t_lat[:, None] + buffer_lat >= centr_lat[None]) + & (t_lon[:, None] - buffer_lon <= centr_lon[None]) + & (t_lon[:, None] + buffer_lon >= centr_lon[None]) + ).any(axis=0).nonzero() + idx_close = idx_close[idx_close_sub] + centr_lat = centr_lat[idx_close_sub] + centr_lon = centr_lon[idx_close_sub] - # extract the velocity from the rescaled momentum through division by r - v_ang[mask_centr_close] = np.fmax(0, momentum / (d_centr + 1e-11)) - return v_ang + # Restrict to metric distance radius around each track position. + # + # We do the distance computation for chunks of the track since computing the distance requires + # npositions*ncentroids*8*3 Bytes of memory. For example, Hurricane FAITH's life time was more + # than 500 hours. At 0.5-hourly resolution and 1,000,000 centroids, that's 24 GB of memory for + # FAITH. With a chunk size of 10, this figure is down to 240 MB. The final along-track mask + # will require 1.0 GB of memory. + chunk_size = 10 + chunks = np.split(np.arange(npositions), np.arange(chunk_size, npositions, chunk_size)) + mask_centr_alongtrack = np.concatenate([ + ( + u_coord.dist_approx( + t_lat[None, chunk], t_lon[None, chunk], + centr_lat[None], centr_lon[None], + normalize=False, method=metric, units="km", + )[0] <= buffer_km + ) for chunk in chunks + ], axis=0) + [idx_close_sub] = mask_centr_alongtrack.any(axis=0).nonzero() + idx_close = idx_close[idx_close_sub] + centr_lat = centr_lat[idx_close_sub] + centr_lon = centr_lon[idx_close_sub] + mask_centr_alongtrack = mask_centr_alongtrack[:, idx_close_sub] + + # Derive mask from index. + mask_centr = np.zeros((ncentroids,), dtype=bool) + mask_centr[idx_close] = True + + centroids_close_normalized = np.stack([centr_lat, centr_lon], axis=1) + return centroids_close_normalized, mask_centr, mask_centr_alongtrack diff --git a/doc/climada/climada.hazard.rst b/doc/climada/climada.hazard.rst index 520c9aae3c..a555ad0083 100644 --- a/doc/climada/climada.hazard.rst +++ b/doc/climada/climada.hazard.rst @@ -64,8 +64,14 @@ climada\.hazard\.tc\_tracks\_synth module climada\.hazard\.trop\_cyclone module ------------------------------------- -.. automodule:: climada.hazard.trop_cyclone +Module for representing the hazard of a tropical cyclone (wind field). + +.. automodule:: climada.hazard.trop_cyclone.trop_cyclone :members: :undoc-members: :show-inheritance: +.. automodule:: climada.hazard.trop_cyclone.trop_cyclone_windfields + :members: + :undoc-members: + :show-inheritance: