From 038959eadca4b1d9520029ac56356dd9d51fb0eb Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Thu, 25 Jan 2024 18:36:00 +0100 Subject: [PATCH 01/13] Introduce a second version of the Holland 2010 model --- climada/hazard/test/test_trop_cyclone.py | 25 ++++++++--- climada/hazard/trop_cyclone.py | 54 ++++++++++++++++++++++-- 2 files changed, 70 insertions(+), 9 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index 992fdeb5bb..067193ede9 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -33,8 +33,8 @@ 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, + _bs_holland_2010_v2, _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.centroids.centr import Centroids import climada.hazard.test as hazard_test @@ -139,19 +139,23 @@ def test_windfield_models(self): 24.745521, 25.596484, 26.475329, 24.690914, 28.650107, 31.584395, 21.723546, 26.140293, 28.94964, 28.051915, 18.49378, 35.312152, ], - # Holland 1980 and Emanuel & Rotunno 2011 use recorded wind speeds, while the above use - # pressure values only. That's why the results are so different: + # Holland 1980, version 2 of Holland 2010 and Emanuel & Rotunno 2011 use recorded wind speeds, + # while the above use pressure values only. That's why the results are so different: "H1980": [21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207], "ER11": [23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795, 18.996389, 26.102109, 30.780737, 29.498453, 0, 38.368805], + "H10_v2": [ + 28.067253, 28.544574, 28.975862, 27.907048, 30.450187, 32.091651, + 25.79227 , 28.140679, 29.615011, 28.693995, 22.191709, 32.245298 + ], } tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK) tc_track.equal_timestep() tc_track.data = tc_track.data[:1] - for model in ["H08", "H10", "H1980", "ER11"]: + for model in ["H08", "H10", "H1980", "ER11", "H10_v2"]: tc_haz = TropCyclone.from_tracks(tc_track, centroids=CENTR_TEST_BRB, model=model) np.testing.assert_array_almost_equal( tc_haz.intensity[0, intensity_idx].toarray()[0], intensity_values[model]) @@ -299,6 +303,17 @@ def test_bs_holland_2008_pass(self): np.testing.assert_array_almost_equal( si_track["hol_b"], [np.nan, 1.27085617, 1.26555341]) + def test_bs_holland_2010_v2_pass(self): + """Test _bs_holland_2010_v2 function. Compare to MATLAB reference.""" + si_track = xr.Dataset({ + "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010, 1010])), + "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 1005.2633, 1005.2682, 1005.2731])), + "vmax": ("time", [np.nan, 22.5, 25.4, 42.5]), + }) + _bs_holland_2010_v2(si_track) + np.testing.assert_array_almost_equal( + si_track["hol_b"], [np.nan, 1.75439, 2.238092, 2.5]) + def test_v_max_s_holland_2008_pass(self): """Test _v_max_s_holland_2008 function.""" # Numbers analogous to test_B_holland_1980_pass diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 125063012b..d7de6fcda9 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -59,11 +59,13 @@ DEF_MAX_MEMORY_GB = 8 """Default value of the memory limit (in GB) for windfield computations (in each thread).""" -MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3} +MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3, 'H10_v2': 4} """Enumerate different symmetric wind field models.""" RHO_AIR = 1.15 """Air density. Assumed constant, following Holland 1980.""" +RHO_AIR_SURFACE = 1.2 +"""Air density at surface level. Assumed constant, following Holland 2010""" GRADIENT_LEVEL_TO_SURFACE_WINDS = 0.9 """Gradient-to-surface wind reduction factor according to the 90%-rule: @@ -969,7 +971,7 @@ def _compute_windfields( v_trans_corr[close_centr_msk] = np.fmin( 1, t_rad_bc[close_centr_msk] / d_centr[close_centr_msk]) - if model in [MODEL_VANG['H08'], MODEL_VANG['H10']]: + if model in [MODEL_VANG['H08'], MODEL_VANG['H10'], MODEL_VANG['H10_v2']]: # 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: @@ -1106,7 +1108,8 @@ def compute_angular_windspeeds(si_track, d_centr, close_centr_msk, model, cyclos _B_holland_1980(si_track) elif model in [MODEL_VANG['H08'], MODEL_VANG['H10']]: _bs_holland_2008(si_track) - + elif model in [MODEL_VANG['H10_v2']]: + _bs_holland_2010_v2(si_track) if model in [MODEL_VANG['H1980'], MODEL_VANG['H08']]: result = _stat_holland_1980( si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic, @@ -1120,6 +1123,9 @@ def compute_angular_windspeeds(si_track, d_centr, close_centr_msk, model, cyclos result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) elif model == MODEL_VANG['ER11']: result = _stat_er_2011(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) + elif model in [MODEL_VANG['H10_v2'],]: + hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk) + result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) else: raise NotImplementedError @@ -1255,6 +1261,8 @@ def _coriolis_parameter(lat: np.ndarray) -> np.ndarray: def _bs_holland_2008(si_track: xr.Dataset): """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". @@ -1311,6 +1319,44 @@ def _bs_holland_2008(si_track: xr.Dataset): ) si_track["hol_b"] = np.clip(si_track["hol_b"], 1, 2.5) +def _bs_holland_2010_v2(si_track: xr.Dataset): + """Holland 2010's second version of how to estimate b-value for sustained surface winds. + For version 1 see "_bs_holland_2008" + + The result is stored in place as a new data variable "hol_b". + + Like the original 1980 formula (see `_B_holland_1980`), this approach does also require + wind speed measurements, if the wind speed measurements are not available or not reliable, consider the second + option proposed in Holland 2010 (and Holland 2008) implemented as "_bs_holland_2008" + + The parameter applies to 1-minute sustained winds at 10 meters above ground. + It is taken from equation (7) in 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 + b_s = vmax^2 * rho^e / ( 100 * (penv - pcen) ) + where: + rho is the air density ρ (in kg m−3) at the surface level + !penv and pcen is assumed to be in hPa in this formula - not Pa, as in our tracks + + Parameters + ---------- + si_track : xr.Dataset + Output of `tctrack_to_si`. The data variables used by this function are + "cen", "env", and "vmax". The result is stored in place as a new data + variable "hol_b". + """ + + si_track["hol_b"] = ( + si_track["vmax"]**2 * RHO_AIR_SURFACE**np.exp(1) + / ( si_track["env"] - si_track["cen"] ) # we do not need the factor 100 as in the original + # formula, because environmental pressure and central pressure are already saved in Pa, not hPa + ) + si_track["hol_b"] = np.clip(si_track["hol_b"], 0.4, 2.5) # reconsider the clip as b_s is b_g/1.6, but does this also relate to limits? + + def _v_max_s_holland_2008(si_track: xr.Dataset): """Compute maximum surface winds from pressure according to Holland 2008. @@ -1433,7 +1479,7 @@ def _x_holland_2010( # compute peripheral exponent from second measurement r_max_norm = (r_max / r_n)**hol_b - x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) + x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) # holland 2010 equation 6 solved for x # linearly interpolate between max exponent and peripheral exponent x_max = 0.5 From a40bb21a3df613c733ec0a1c284f29bb5da55c21 Mon Sep 17 00:00:00 2001 From: Thomas Vogt Date: Fri, 26 Jan 2024 14:05:45 +0100 Subject: [PATCH 02/13] TropCyclone: implement model_kwargs --- climada/hazard/test/test_trop_cyclone.py | 98 +++-- climada/hazard/trop_cyclone.py | 477 +++++++++++++++++------ 2 files changed, 416 insertions(+), 159 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index 067193ede9..26a85cde25 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -33,8 +33,8 @@ from climada.hazard.tc_tracks import TCTracks from climada.hazard.trop_cyclone import ( TropCyclone, get_close_centroids, _vtrans, _B_holland_1980, _bs_holland_2008, - _bs_holland_2010_v2, _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, + _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.centroids.centr import Centroids import climada.hazard.test as hazard_test @@ -130,36 +130,47 @@ def test_set_one_pass(self): def test_windfield_models(self): """Test _tc_from_track function with different wind field models.""" intensity_idx = [0, 1, 2, 3, 80, 100, 120, 200, 220, 250, 260, 295] - intensity_values = { - "H08": [ - 22.74903, 23.784691, 24.82255, 22.67403, 27.218706, 30.593959, - 18.980878, 24.540069, 27.826407, 26.846293, 0., 34.568898, - ], - "H10": [ - 24.745521, 25.596484, 26.475329, 24.690914, 28.650107, 31.584395, - 21.723546, 26.140293, 28.94964, 28.051915, 18.49378, 35.312152, - ], - # Holland 1980, version 2 of Holland 2010 and Emanuel & Rotunno 2011 use recorded wind speeds, - # while the above use pressure values only. That's why the results are so different: - "H1980": [21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, - 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207], - "ER11": [23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795, - 18.996389, 26.102109, 30.780737, 29.498453, 0, 38.368805], - "H10_v2": [ - 28.067253, 28.544574, 28.975862, 27.907048, 30.450187, 32.091651, - 25.79227 , 28.140679, 29.615011, 28.693995, 22.191709, 32.245298 - ], - } + intensity_values = [ + ("H08", None, [ + 22.74903, 23.784691, 24.82255, 22.67403, 27.218706, 30.593959, + 18.980878, 24.540069, 27.826407, 26.846293, 0., 34.568898, + ]), + ("H10", None, [ + 24.745521, 25.596484, 26.475329, 24.690914, 28.650107, 31.584395, + 21.723546, 26.140293, 28.94964, 28.051915, 18.49378, 35.312152, + ]), + # The following model configurations use recorded wind speeds, while the above use + # pressure values only. That's why some of the values are so different. + ("H10", dict(vmax_from_cen=False, rho_air_const=1.2), [ + 23.702232, 24.327615, 24.947161, 23.589233, 26.616085, 29.389295, + 21.338178, 24.257067, 26.472543, 25.662313, 18.535842, 31.886041, + ]), + ("H10", dict(vmax_from_cen=False, rho_air_const=None), [ + 24.244162, 24.835561, 25.432454, 24.139294, 27.127457, 29.719196, + 21.910658, 24.692637, 26.783575, 25.971516, 19.005555, 31.904048, + ]), + ("H1980", None, [ + 21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, + 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207, + ]), + ("ER11", None, [ + 23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795, + 18.996389, 26.102109, 30.780737, 29.498453, 0, 38.368805, + ]), + ] tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK) tc_track.equal_timestep() tc_track.data = tc_track.data[:1] - for model in ["H08", "H10", "H1980", "ER11", "H10_v2"]: - tc_haz = TropCyclone.from_tracks(tc_track, centroids=CENTR_TEST_BRB, model=model) + for model, model_kwargs, inten_ref in intensity_values: + tc_haz = TropCyclone.from_tracks( + tc_track, centroids=CENTR_TEST_BRB, model=model, model_kwargs=model_kwargs, + ) np.testing.assert_array_almost_equal( - tc_haz.intensity[0, intensity_idx].toarray()[0], intensity_values[model]) - for idx, val in zip(intensity_idx, intensity_values[model]): + tc_haz.intensity[0, intensity_idx].toarray()[0], inten_ref, + ) + for idx, val in zip(intensity_idx, inten_ref): if val == 0: self.assertEqual(tc_haz.intensity[0, idx], 0) @@ -172,10 +183,14 @@ def test_windfield_models_different_windunits(self): intensity_values = { # Holland 1980 and Emanuel & Rotunno 2011 use recorded wind speeds, that is why checking them for different # windspeed units is so important: - "H1980": [21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, - 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207], - "ER11": [23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795, - 18.996389, 26.102109, 30.780737, 29.498453, 0, 38.368805], + "H1980": [ + 21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, + 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207, + ], + "ER11": [ + 23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795, + 18.996389, 26.102109, 30.780737, 29.498453, 0, 38.368805, + ], } tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK) @@ -286,10 +301,20 @@ def test_B_holland_1980_pass(self): "env": ("time", MBAR_TO_PA * np.array([1010, 1010])), "cen": ("time", MBAR_TO_PA * np.array([995, 980])), "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({ + "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010, 1010])), + "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 995, 980, 970])), + "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({ @@ -303,17 +328,6 @@ def test_bs_holland_2008_pass(self): np.testing.assert_array_almost_equal( si_track["hol_b"], [np.nan, 1.27085617, 1.26555341]) - def test_bs_holland_2010_v2_pass(self): - """Test _bs_holland_2010_v2 function. Compare to MATLAB reference.""" - si_track = xr.Dataset({ - "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010, 1010])), - "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 1005.2633, 1005.2682, 1005.2731])), - "vmax": ("time", [np.nan, 22.5, 25.4, 42.5]), - }) - _bs_holland_2010_v2(si_track) - np.testing.assert_array_almost_equal( - si_track["hol_b"], [np.nan, 1.75439, 2.238092, 2.5]) - def test_v_max_s_holland_2008_pass(self): """Test _v_max_s_holland_2008 function.""" # Numbers analogous to test_B_holland_1980_pass @@ -321,6 +335,7 @@ def test_v_max_s_holland_2008_pass(self): "env": ("time", MBAR_TO_PA * np.array([1010, 1010])), "cen": ("time", MBAR_TO_PA * np.array([995, 980])), "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]) @@ -396,6 +411,7 @@ def test_stat_holland_1980(self): "cen": ("time", MBAR_TO_PA * np.array([970.8727666672957, 1005.268166666671])), "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) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index d7de6fcda9..ba667ea9f8 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -59,20 +59,21 @@ DEF_MAX_MEMORY_GB = 8 """Default value of the memory limit (in GB) for windfield computations (in each thread).""" -MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3, 'H10_v2': 4} +MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3} """Enumerate different symmetric wind field models.""" -RHO_AIR = 1.15 -"""Air density. Assumed constant, following Holland 1980.""" -RHO_AIR_SURFACE = 1.2 -"""Air density at surface level. Assumed constant, following Holland 2010""" +DEF_RHO_AIR = 1.15 +"""Default value for air density (in kg/m³), following Holland 1980.""" -GRADIENT_LEVEL_TO_SURFACE_WINDS = 0.9 -"""Gradient-to-surface wind reduction factor according to the 90%-rule: +DEF_GRADIENT_TO_SURFACE_WINDS = 0.9 +"""Default gradient-to-surface wind reduction factor, following the 90%-rule mentioned in: Franklin, J.L., Black, M.L., Valde, K. (2003): GPS Dropwindsonde Wind Profiles in Hurricanes and Their Operational Implications. Weather and Forecasting 18(1): 32–44. https://doi.org/10.1175/1520-0434(2003)018<0032:GDWPIH>2.0.CO;2 + +According to Table 2, this is a reasonable factor for the 750 hPa level in the eyewall region. For +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 @@ -83,6 +84,9 @@ 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)""" @@ -183,6 +187,7 @@ def from_tracks( centroids: Optional[Centroids] = None, 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", @@ -224,6 +229,28 @@ def from_tracks( "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 wind model. None of the parameters is + currently supported by the ER11 model. 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 + + Default: None ignore_distance_to_coast : boolean, optional If True, centroids far from coast are not ignored. Default: False. store_windfields : boolean, optional @@ -296,21 +323,32 @@ def from_tracks( & (t_lat_min <= coastal_centroids[:, 0]) & (coastal_centroids[:, 0] <= t_lat_max))] - LOGGER.info('Mapping %s tracks to %s coastal centroids.', str(tracks.size), - str(coastal_idx.size)) + # prepare keyword arguments to pass to `from_single_track` + kwargs_from_single_track = dict( + centroids=centroids, + coastal_idx=coastal_idx, + 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, coastal_idx.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, - itertools.repeat(centroids, num_tracks), - itertools.repeat(coastal_idx, num_tracks), - itertools.repeat(model, num_tracks), - itertools.repeat(store_windfields, num_tracks), - itertools.repeat(metric, num_tracks), - itertools.repeat(intensity_thres, num_tracks), - itertools.repeat(max_dist_eye_km, num_tracks), - itertools.repeat(max_memory_gb, num_tracks), - chunksize=chunksize) + cls.from_single_track, + tracks.data, + *kwargs_repeated, + chunksize=chunksize, + ) else: last_perc = 0 tc_haz_list = [] @@ -320,11 +358,8 @@ def from_tracks( LOGGER.info("Progress: %d%%", perc) last_perc = perc tc_haz_list.append( - cls.from_single_track(track, centroids, coastal_idx, - model=model, store_windfields=store_windfields, - metric=metric, intensity_thres=intensity_thres, - max_dist_eye_km=max_dist_eye_km, - max_memory_gb=max_memory_gb)) + cls.from_single_track(track, **kwargs_from_single_track) + ) if last_perc < 100: LOGGER.info("Progress: 100%") @@ -511,6 +546,7 @@ def from_single_track( centroids: Centroids, coastal_idx: np.ndarray, model: str = 'H08', + model_kwargs: Optional[dict] = None, store_windfields: bool = False, metric: str = "equirect", intensity_thres: float = DEF_INTENSITY_THRES, @@ -533,6 +569,8 @@ def from_single_track( "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 @@ -561,6 +599,7 @@ def from_single_track( centroids=centroids, coastal_idx=coastal_idx, model=model, + model_kwargs=model_kwargs, store_windfields=store_windfields, metric=metric, intensity_thres=intensity_thres, @@ -667,6 +706,7 @@ def _compute_windfields_sparse( centroids: Centroids, coastal_idx: np.ndarray, model: str = 'H08', + model_kwargs: Optional[dict] = None, store_windfields: bool = False, metric: str = "equirect", intensity_thres: float = DEF_INTENSITY_THRES, @@ -688,6 +728,8 @@ def _compute_windfields_sparse( "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 @@ -787,6 +829,7 @@ def _compute_windfields_sparse( centroids, coastal_idx, model=model, + model_kwargs=model_kwargs, store_windfields=store_windfields, metric=metric, intensity_thres=intensity_thres, @@ -795,7 +838,12 @@ def _compute_windfields_sparse( ) windfields, reachable_centr_idx = _compute_windfields( - si_track, coastal_centr, mod_id, metric=metric, max_dist_eye_km=max_dist_eye_km, + si_track, + coastal_centr, + mod_id, + model_kwargs=model_kwargs, + metric=metric, + max_dist_eye_km=max_dist_eye_km, ) reachable_coastal_centr_idx = coastal_idx[reachable_centr_idx] npositions = windfields.shape[0] @@ -887,6 +935,7 @@ 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]: @@ -907,6 +956,8 @@ def _compute_windfields( Centroids that are not within reach of the track are ignored. 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`. @@ -954,10 +1005,9 @@ def _compute_windfields( # derive (absolute) angular velocity from parametric wind profile v_ang_norm = compute_angular_windspeeds( - si_track, d_centr, close_centr_msk, model, cyclostrophic=False, + si_track, d_centr, close_centr_msk, model, model_kwargs=model_kwargs, cyclostrophic=False, ) - # Influence of translational speed decreases with distance from eye. # The "absorbing factor" is according to the following paper (see Fig. 7): # @@ -971,7 +1021,7 @@ def _compute_windfields( v_trans_corr[close_centr_msk] = np.fmin( 1, t_rad_bc[close_centr_msk] / d_centr[close_centr_msk]) - if model in [MODEL_VANG['H08'], MODEL_VANG['H10'], MODEL_VANG['H10_v2']]: + 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: @@ -1073,16 +1123,35 @@ def tctrack_to_si( # add translational speed of track at every node (in m/s) _vtrans(si_track, metric=metric) - # convert surface winds to gradient winds without translational influence - si_track["vgrad"] = ( - np.fmax(0, si_track["vmax"] - si_track["vtrans_norm"]) / GRADIENT_LEVEL_TO_SURFACE_WINDS - ) - si_track["cp"] = ("time", _coriolis_parameter(si_track["lat"].values)) return si_track -def compute_angular_windspeeds(si_track, d_centr, close_centr_msk, model, cyclostrophic=False): +def _vgrad(si_track, gradient_to_surface_winds): + """Gradient wind speeds (in m/s) without translational influence at each track node + + The track dataset is modified in place, with the "vgrad" data variable added. + + 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. + """ + si_track["vgrad"] = ( + np.fmax(0, si_track["vmax"] - si_track["vtrans_norm"]) / gradient_to_surface_winds + ) + +def compute_angular_windspeeds( + si_track: xr.Dataset, + d_centr: np.ndarray, + close_centr_msk: np.ndarray, + model: int, + model_kwargs: Optional[dict] = None, + cyclostrophic: bool = False, +): """Compute (absolute) angular wind speeds according to a parametric wind profile Parameters @@ -1096,6 +1165,8 @@ def compute_angular_windspeeds(si_track, d_centr, close_centr_msk, model, cyclos 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, don't apply the influence of the Coriolis force (set the Coriolis terms to 0). Default: False @@ -1104,35 +1175,157 @@ def compute_angular_windspeeds(si_track, d_centr, close_centr_msk, model, cyclos ------- ndarray of shape (npositions, ncentroids) """ - if model == MODEL_VANG['H1980']: - _B_holland_1980(si_track) - elif model in [MODEL_VANG['H08'], MODEL_VANG['H10']]: - _bs_holland_2008(si_track) - elif model in [MODEL_VANG['H10_v2']]: - _bs_holland_2010_v2(si_track) - if model in [MODEL_VANG['H1980'], MODEL_VANG['H08']]: - result = _stat_holland_1980( - si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic, - ) - if model == MODEL_VANG['H1980']: - result *= GRADIENT_LEVEL_TO_SURFACE_WINDS - elif model == MODEL_VANG['H10']: - # this model is always cyclostrophic - _v_max_s_holland_2008(si_track) - hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk) - result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) - elif model == MODEL_VANG['ER11']: - result = _stat_er_2011(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic) - elif model in [MODEL_VANG['H10_v2'],]: - hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk) - result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) - else: - raise NotImplementedError - + 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, + close_centr_msk, + cyclostrophic=cyclostrophic, + **model_kwargs, + ) 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 + + 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 + Returns + ------- + ndarray of shape (npositions, ncentroids) + """ + _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 +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 + ---------- + 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 + ------- + ndarray of shape (npositions, ncentroids) + """ + _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) + +def _compute_angular_windspeeds_h10( + 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, + vmax_from_cen: bool = True, +): + """Compute (absolute) angular wind speeds according to the Holland et al. 2010 model + + 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 + 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: 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 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 + + Returns + ------- + ndarray of shape (npositions, ncentroids) + """ + _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) + return _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) + def get_close_centroids( t_lat: np.ndarray, t_lon: np.ndarray, @@ -1259,7 +1452,51 @@ def _coriolis_parameter(lat: np.ndarray) -> np.ndarray: """ return 2 * V_ANG_EARTH * np.sin(np.radians(np.abs(lat))) -def _bs_holland_2008(si_track: xr.Dataset): +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. + + Parameters + ---------- + si_track : xr.Dataset + Track information as returned by `tctrack_to_si`. The data variables used by this function + are "lat" and "cen". 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["env"] - si_track["cen"]) / 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'.) @@ -1299,6 +1536,9 @@ def _bs_holland_2008(si_track: xr.Dataset): Output of `tctrack_to_si`. The data variables used by this function are "lat", "tstep", "vtrans_norm", "cen", and "env". 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) @@ -1317,45 +1557,8 @@ def _bs_holland_2008(si_track: xr.Dataset): - 0.014 * abs(si_track["lat"]) + 0.15 * si_track["vtrans_norm"]**hol_xx + 1.0 ) - si_track["hol_b"] = np.clip(si_track["hol_b"], 1, 2.5) - -def _bs_holland_2010_v2(si_track: xr.Dataset): - """Holland 2010's second version of how to estimate b-value for sustained surface winds. - For version 1 see "_bs_holland_2008" - - The result is stored in place as a new data variable "hol_b". - - Like the original 1980 formula (see `_B_holland_1980`), this approach does also require - wind speed measurements, if the wind speed measurements are not available or not reliable, consider the second - option proposed in Holland 2010 (and Holland 2008) implemented as "_bs_holland_2008" - - The parameter applies to 1-minute sustained winds at 10 meters above ground. - It is taken from equation (7) in 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 - b_s = vmax^2 * rho^e / ( 100 * (penv - pcen) ) - where: - rho is the air density ρ (in kg m−3) at the surface level - !penv and pcen is assumed to be in hPa in this formula - not Pa, as in our tracks - - Parameters - ---------- - si_track : xr.Dataset - Output of `tctrack_to_si`. The data variables used by this function are - "cen", "env", and "vmax". The result is stored in place as a new data - variable "hol_b". - """ - - si_track["hol_b"] = ( - si_track["vmax"]**2 * RHO_AIR_SURFACE**np.exp(1) - / ( si_track["env"] - si_track["cen"] ) # we do not need the factor 100 as in the original - # formula, because environmental pressure and central pressure are already saved in Pa, not hPa - ) - si_track["hol_b"] = np.clip(si_track["hol_b"], 0.4, 2.5) # reconsider the clip as b_s is b_g/1.6, but does this also relate to limits? - + 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. @@ -1378,15 +1581,18 @@ def _v_max_s_holland_2008(si_track: xr.Dataset): 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 "env", "cen", and "hol_b". The results are stored in place as - a new data variable "vmax". If a variable of that name already exists, its values are - overwritten. + 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 "env", "cen", "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. """ pdelta = si_track["env"] - si_track["cen"] - si_track["vmax"] = np.sqrt(si_track["hol_b"] / (RHO_AIR * np.exp(1)) * pdelta) + si_track["vmax"] = np.sqrt(si_track["hol_b"] / (si_track["rho_air"] * np.exp(1)) * pdelta) -def _B_holland_1980(si_track: xr.Dataset): # pylint: disable=invalid-name +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". @@ -1408,13 +1614,46 @@ def _B_holland_1980(si_track: xr.Dataset): # pylint: disable=invalid-name Parameters ---------- si_track : xr.Dataset - Output of `tctrack_to_si` with "vgrad" variable (see _vgrad). The data variables - used by this function are "vgrad", "env", and "cen". The results are stored in place as - a new data variable "hol_b". + 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), "env", "cen", 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" + pdelta = si_track["env"] - si_track["cen"] - si_track["hol_b"] = si_track["vgrad"]**2 * np.exp(1) * RHO_AIR / np.fmax(np.spacing(1), pdelta) - si_track["hol_b"] = np.clip(si_track["hol_b"], 1, 2.5) + si_track["hol_b"] = ( + si_track[windvar]**2 * np.exp(1) * si_track["rho_air"] / np.fmax(np.spacing(1), pdelta) + ) + + 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 + ---------- + 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 + ------- + b_min, b_max : float + Minimum and maximum value of the clip interval. + """ + 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, @@ -1478,8 +1717,9 @@ def _x_holland_2010( 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 - x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) # holland 2010 equation 6 solved for x + 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 @@ -1563,9 +1803,8 @@ def _stat_holland_1980( 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 is assumed to be constant while the Coriolis parameter f is computed - from the latitude `lat` using the constant rotation rate of the earth. + 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 @@ -1574,8 +1813,9 @@ def _stat_holland_1980( Parameters ---------- si_track : xr.Dataset - Output of `tctrack_to_si` with "hol_b" (see, e.g., _B_holland_1980) data variable. The data - variables used by this function are "lat", "cp", "rad", "cen", "env", and "hol_b". + 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", "cen", "env", "hol_b", and "rho_air". d_centr : np.ndarray of shape (nnodes, ncentroids) Distance (in m) between centroids and track nodes. close_centr : np.ndarray of shape (nnodes, ncentroids) @@ -1590,7 +1830,7 @@ def _stat_holland_1980( Absolute values of wind speeds (m/s) in angular direction. """ v_ang = np.zeros_like(d_centr) - r_max, hol_b, penv, pcen, coriolis_p, d_centr = [ + r_max, hol_b, penv, pcen, coriolis_p, rho_air, d_centr = [ np.broadcast_to(ar, d_centr.shape)[close_centr] for ar in [ si_track["rad"].values[:, None], @@ -1598,6 +1838,7 @@ def _stat_holland_1980( si_track["env"].values[:, None], si_track["cen"].values[:, None], si_track["cp"].values[:, None], + si_track["rho_air"].values[:, None], d_centr, ] ] @@ -1607,7 +1848,7 @@ def _stat_holland_1980( r_coriolis = 0.5 * d_centr * coriolis_p r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b - sqrt_term = hol_b / RHO_AIR * r_max_norm * (penv - pcen) * np.exp(-r_max_norm) + r_coriolis**2 + sqrt_term = hol_b / rho_air * r_max_norm * (penv - pcen) * np.exp(-r_max_norm) + r_coriolis**2 v_ang[close_centr] = np.sqrt(np.fmax(0, sqrt_term)) - r_coriolis return v_ang From 5b833adf7a0e5004f3c1bea56e9a30168768634e Mon Sep 17 00:00:00 2001 From: Thomas Vogt Date: Wed, 7 Feb 2024 18:10:28 +0100 Subject: [PATCH 03/13] trop_cyclone: add pdelta to si_track --- climada/hazard/test/test_trop_cyclone.py | 31 +++++++++----------- climada/hazard/trop_cyclone.py | 37 ++++++++++++++---------- 2 files changed, 35 insertions(+), 33 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index 26a85cde25..c70a850aaf 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -298,8 +298,7 @@ def test_get_close_centroids_pass(self): def test_B_holland_1980_pass(self): """Test _B_holland_1980 function.""" si_track = xr.Dataset({ - "env": ("time", MBAR_TO_PA * np.array([1010, 1010])), - "cen": ("time", MBAR_TO_PA * np.array([995, 980])), + "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), "vgrad": ("time", [35, 40]), "rho_air": ("time", [1.15, 1.15]) }) @@ -307,8 +306,7 @@ def test_B_holland_1980_pass(self): np.testing.assert_array_almost_equal(si_track["hol_b"], [2.5, 1.667213]) si_track = xr.Dataset({ - "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010, 1010])), - "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 995, 980, 970])), + "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]) }) @@ -320,20 +318,20 @@ def test_bs_holland_2008_pass(self): 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]), - "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010])), + "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_array_almost_equal( - si_track["hol_b"], [np.nan, 1.27085617, 1.26555341]) + 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({ - "env": ("time", MBAR_TO_PA * np.array([1010, 1010])), - "cen": ("time", MBAR_TO_PA * np.array([995, 980])), + "pdelta": ("time", MBAR_TO_PA * np.array([15, 30])), "hol_b": ("time", [2.5, 1.67]), "rho_air": ("time", [1.15, 1.15]), }) @@ -407,23 +405,22 @@ def test_stat_holland_1980(self): si_track = xr.Dataset({ "rad": ("time", KM_TO_M * np.array([40.665454622610511, 75.547902916671745])), "hol_b": ("time", [1.486076257880692, 1.265551666104679]), - "env": ("time", MBAR_TO_PA * np.array([1010.0, 1010.0])), - "cen": ("time", MBAR_TO_PA * np.array([970.8727666672957, 1005.268166666671])), + "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_array_almost_equal(v_ang_norm, - [[11.279764005440288, 11.682978583939310, 11.610940769149384, 42.412845], - [5.384115724400597, 0, 5.281356766052531, 12.763087]]) + 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_array_almost_equal(v_ang_norm, - [[15.719924, 16.037052, 15.980323, 43.128461], - [8.836768, 0, 8.764678, 13.807452]]) + 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) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index ba667ea9f8..3f4ebf51e2 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -1053,6 +1053,7 @@ def tctrack_to_si( are normalized and additional variables are defined: * cp (coriolis parameter) + * pdelta (difference between environmental and central pressure, always strictly positive) * vtrans (translational velocity vectors) * vtrans_norm (absolute value of translational speed) @@ -1123,8 +1124,12 @@ def tctrack_to_si( # add translational speed of track at every node (in m/s) _vtrans(si_track, metric=metric) + # add Coriolis parameter si_track["cp"] = ("time", _coriolis_parameter(si_track["lat"].values)) + # add pressure drop + si_track["pdelta"] = np.fmax(np.spacing(1), si_track["env"] - si_track["cen"]) + return si_track def _vgrad(si_track, gradient_to_surface_winds): @@ -1461,7 +1466,8 @@ def _rho_air(si_track: xr.Dataset, const: Optional[float]): ---------- si_track : xr.Dataset Track information as returned by `tctrack_to_si`. The data variables used by this function - are "lat" and "cen". The result is stored in place as new data variable "rho_air". + 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. @@ -1477,7 +1483,7 @@ def _rho_air(si_track: xr.Dataset, const: Optional[float]): 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["env"] - si_track["cen"]) / np.exp(1) + 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) @@ -1534,7 +1540,7 @@ def _bs_holland_2008( ---------- si_track : xr.Dataset Output of `tctrack_to_si`. The data variables used by this function are "lat", "tstep", - "vtrans_norm", "cen", and "env". The result is stored in place as a new data + "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 @@ -1549,7 +1555,7 @@ def _bs_holland_2008( # 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["env"] - si_track["cen"]) / MBAR_TO_PA + 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 @@ -1582,12 +1588,13 @@ def _v_max_s_holland_2008(si_track: xr.Dataset): ---------- 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 "env", "cen", "hol_b" + _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. """ - pdelta = si_track["env"] - si_track["cen"] - si_track["vmax"] = np.sqrt(si_track["hol_b"] / (si_track["rho_air"] * np.exp(1)) * pdelta) + si_track["vmax"] = np.sqrt( + si_track["hol_b"] / (si_track["rho_air"] * np.exp(1)) * si_track["pdelta"] + ) def _B_holland_1980( # pylint: disable=invalid-name si_track: xr.Dataset, @@ -1616,17 +1623,16 @@ def _B_holland_1980( # pylint: disable=invalid-name 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), "env", "cen", and "rho_air". The results are stored in place as a new data variable - "hol_b". + 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" - pdelta = si_track["env"] - si_track["cen"] si_track["hol_b"] = ( - si_track[windvar]**2 * np.exp(1) * si_track["rho_air"] / np.fmax(np.spacing(1), pdelta) + si_track[windvar]**2 * np.exp(1) * si_track["rho_air"] / si_track["pdelta"] ) clip_interval = _b_holland_clip_interval(gradient_to_surface_winds) @@ -1815,7 +1821,7 @@ def _stat_holland_1980( 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", "cen", "env", "hol_b", and "rho_air". + "cp", "rad", "pdelta", "hol_b", and "rho_air". d_centr : np.ndarray of shape (nnodes, ncentroids) Distance (in m) between centroids and track nodes. close_centr : np.ndarray of shape (nnodes, ncentroids) @@ -1830,13 +1836,12 @@ def _stat_holland_1980( Absolute values of wind speeds (m/s) in angular direction. """ v_ang = np.zeros_like(d_centr) - r_max, hol_b, penv, pcen, coriolis_p, rho_air, d_centr = [ + r_max, hol_b, pdelta, coriolis_p, rho_air, d_centr = [ np.broadcast_to(ar, d_centr.shape)[close_centr] for ar in [ si_track["rad"].values[:, None], si_track["hol_b"].values[:, None], - si_track["env"].values[:, None], - si_track["cen"].values[:, None], + si_track["pdelta"].values[:, None], si_track["cp"].values[:, None], si_track["rho_air"].values[:, None], d_centr, @@ -1848,7 +1853,7 @@ def _stat_holland_1980( r_coriolis = 0.5 * d_centr * coriolis_p r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b - sqrt_term = hol_b / rho_air * r_max_norm * (penv - pcen) * np.exp(-r_max_norm) + r_coriolis**2 + sqrt_term = hol_b / rho_air * r_max_norm * pdelta * np.exp(-r_max_norm) + r_coriolis**2 v_ang[close_centr] = np.sqrt(np.fmax(0, sqrt_term)) - r_coriolis return v_ang From 027d3ee288734e5615d5e3beb0767795fbc1fec7 Mon Sep 17 00:00:00 2001 From: Thomas Vogt Date: Tue, 20 Feb 2024 16:30:21 +0100 Subject: [PATCH 04/13] trop_cyclone: truncate hol_x parameter --- climada/hazard/test/test_trop_cyclone.py | 20 ++++++++++---------- climada/hazard/trop_cyclone.py | 9 ++++----- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index c70a850aaf..633a94abf8 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -142,12 +142,12 @@ def test_windfield_models(self): # The following model configurations use recorded wind speeds, while the above use # pressure values only. That's why some of the values are so different. ("H10", dict(vmax_from_cen=False, rho_air_const=1.2), [ - 23.702232, 24.327615, 24.947161, 23.589233, 26.616085, 29.389295, - 21.338178, 24.257067, 26.472543, 25.662313, 18.535842, 31.886041, + 23.702232, 24.327615, 24.947161, 23.589233, 26.616085, 29.389295, + 21.338178, 24.257067, 26.472543, 25.662313, 18.535842, 31.886041, ]), ("H10", dict(vmax_from_cen=False, rho_air_const=None), [ - 24.244162, 24.835561, 25.432454, 24.139294, 27.127457, 29.719196, - 21.910658, 24.692637, 26.783575, 25.971516, 19.005555, 31.904048, + 24.244162, 24.835561, 25.432454, 24.139294, 27.127457, 29.719196, + 21.910658, 24.692637, 26.783575, 25.971516, 19.005555, 31.904048, ]), ("H1980", None, [ 21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, @@ -380,21 +380,21 @@ def test_holland_2010_pass(self): 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.289203, 0.000000], - [0.5, 0.497620, 0.131072, 0.000000, 0.000000], + [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_array_almost_equal(v_ang_norm, [ + 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.0000000, 35.000000, 21.181497, 17.00000, 12.103461], - [1.2964800, 34.990037, 21.593755, 17.00000, 0.0000000], - [0.3219518, 15.997500, 13.585498, 16.00000, 16.000000], + [1.2964800, 34.990037, 21.593755, 12.89131, 0.0000000], + [0.3219518, 15.997500, 9.7120060, 8.087240, 6.2289690], [24.823469, 89.992938, 24.381965, 17.00000, 1.9292020], - ]) + ], atol=1e-5) def test_stat_holland_1980(self): """Test _stat_holland_1980 function. Compare to MATLAB reference.""" diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 3f4ebf51e2..49e5c7b675 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -1731,11 +1731,10 @@ def _x_holland_2010( x_max = 0.5 hol_x[close_centr] = x_max + np.fmax(0, d_centr - r_max) * (x_n - x_max) / (r_n - r_max) - # Negative hol_x values appear when v_max_s is very close to or even lower than v_n (which - # should never happen in theory). In those cases, wind speeds might decrease outside of the eye - # wall and increase again towards the peripheral radius (which is actually unphysical). - # We clip hol_x to 0, otherwise wind speeds keep increasing indefinitely away from the eye: - hol_x[close_centr] = np.fmax(hol_x[close_centr], 0.0) + # 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[close_centr] = np.fmax(hol_x[close_centr], 0.4) return hol_x From 11f02d44591a8e1f73d37c47479671d8c8702352 Mon Sep 17 00:00:00 2001 From: Thomas Vogt Date: Tue, 20 Feb 2024 17:01:39 +0100 Subject: [PATCH 05/13] trop_cyclone: implement vmax_in_brackets model kwarg --- climada/hazard/test/test_trop_cyclone.py | 4 ++ climada/hazard/trop_cyclone.py | 51 ++++++++++++++++++++++-- 2 files changed, 51 insertions(+), 4 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index 85c30656c6..0b1fe859cf 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -168,6 +168,10 @@ def test_windfield_models(self): 24.244162, 24.835561, 25.432454, 24.139294, 27.127457, 29.719196, 21.910658, 24.692637, 26.783575, 25.971516, 19.005555, 31.904048, ]), + ("H10", dict(vmax_from_cen=False, rho_air_const=None, vmax_in_brackets=True), [ + 23.592924, 24.208169, 24.817104, 23.483053, 26.468975, 29.221715, + 21.260867, 24.150879, 26.34288 , 25.543635, 18.487385, 31.904048 + ]), ("H1980", None, [ 21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303, 19.220149, 21.984516, 24.196388, 23.449116, 0, 31.550207, diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index fd3a9e4a60..4f2490ed72 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -249,6 +249,12 @@ def from_tracks( 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 Default: None ignore_distance_to_coast : boolean, optional @@ -1279,6 +1285,7 @@ def _compute_angular_windspeeds_h10( 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 @@ -1306,6 +1313,12 @@ def _compute_angular_windspeeds_h10( 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 ------- @@ -1317,8 +1330,10 @@ def _compute_angular_windspeeds_h10( _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) - return _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x) + 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, + ) def get_close_centroids( si_track: xr.Dataset, @@ -1698,6 +1713,7 @@ def _x_holland_2010( 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. @@ -1731,6 +1747,11 @@ def _x_holland_2010( 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 ------- @@ -1756,7 +1777,16 @@ def _x_holland_2010( # 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 - x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm)) + 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 @@ -1774,6 +1804,7 @@ def _stat_holland_2010( 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 @@ -1788,6 +1819,10 @@ def _stat_holland_2010( 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: + + V(r) = [v_max_s^2 * (r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x + Parameters ---------- si_track : xr.Dataset @@ -1799,6 +1834,11 @@ def _stat_holland_2010( 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 Returns ------- @@ -1818,7 +1858,10 @@ def _stat_holland_2010( ] r_max_norm = (r_max / np.fmax(1, d_centr))**hol_b - v_ang[mask_centr_close] = v_max_s * (r_max_norm * np.exp(1 - r_max_norm))**hol_x + 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 def _stat_holland_1980( From 9c16bc8a16b1b46ad25e21a6202679ed5f0a4941 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20R=C3=B6=C3=B6sli?= Date: Tue, 21 May 2024 18:23:55 +0200 Subject: [PATCH 06/13] Use math in doc string Co-authored-by: Lukas Riedel <34276446+peanutfun@users.noreply.github.com> --- climada/hazard/trop_cyclone.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 07e69e8c58..5f8ed14088 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -1821,9 +1821,11 @@ def _stat_holland_2010( 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: + If ``vmax_in_brackets`` is True, the alternative formula in (6) is used: - V(r) = [v_max_s^2 * (r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x + .. math:: + + V(r) = [v_max_s^2 * (r_max / r)^b_s * e^(1 - (r_max / r)^b_s)]^x Parameters ---------- From 147ad92f553983c8f31c7d7b9fb4300fc74ce672 Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Tue, 21 May 2024 18:44:18 +0200 Subject: [PATCH 07/13] Update docstring with DOI --- climada/hazard/trop_cyclone.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 5f8ed14088..9a6307769d 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -225,9 +225,23 @@ def from_tracks( description : str, optional Description of the event set. Default: "". model : str, optional - Parametric wind field model to use: 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). + Parametric wind field model to use: + "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, + or + "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. Default: "H08". model_kwargs : dict, optional If given, forward these kwargs to the selected wind model. None of the parameters is From 8e409f84f546d8c924e53783df189bf5c4efb905 Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Tue, 21 May 2024 18:47:18 +0200 Subject: [PATCH 08/13] remove docstring of inextisting parameter --- climada/hazard/trop_cyclone.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 9a6307769d..3288d6274e 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -222,8 +222,6 @@ def from_tracks( 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 - description : str, optional - Description of the event set. Default: "". model : str, optional Parametric wind field model to use: "H1980" (the prominent Holland 1980 model) from the paper: From 91adf6529d08d34cbd05d5325c3977ef9dfe7acb Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Tue, 21 May 2024 19:12:09 +0200 Subject: [PATCH 09/13] add description of returned array for angular windspeed functions --- climada/hazard/trop_cyclone.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index 3288d6274e..a5155ff066 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -1188,6 +1188,7 @@ def compute_angular_windspeeds( Returns ------- ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ model_kwargs = {} if model_kwargs is None else model_kwargs compute_funs = { @@ -1243,6 +1244,7 @@ def _compute_angular_windspeeds_h1980( 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) @@ -1286,6 +1288,7 @@ def _compute_angular_windspeeds_h08( Returns ------- ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ _rho_air(si_track, rho_air_const) _bs_holland_2008(si_track, gradient_to_surface_winds=gradient_to_surface_winds) @@ -1337,6 +1340,7 @@ def _compute_angular_windspeeds_h10( Returns ------- ndarray of shape (npositions, ncentroids) + containing the magnitude of the angular windspeed per track position per centroid location """ _rho_air(si_track, rho_air_const) if vmax_from_cen: From 4f83c006ba24c3257b099dfc83f6cfcab8f61313 Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Tue, 21 May 2024 19:57:15 +0200 Subject: [PATCH 10/13] correct values and tolerance for test holland 2010 pass --- climada/hazard/test/test_trop_cyclone.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index ad0f710622..b58bf268eb 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -437,11 +437,11 @@ def test_holland_2010_pass(self): # 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.0000000, 35.000000, 21.181497, 17.00000, 12.103461], - [1.2964800, 34.990037, 21.593755, 12.89131, 0.0000000], - [0.3219518, 15.997500, 9.7120060, 8.087240, 6.2289690], - [24.823469, 89.992938, 24.381965, 17.00000, 1.9292020], - ], atol=1e-5) + [ 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.""" From ca8338dfd5e54bca67444c4cf905073af1e882ab Mon Sep 17 00:00:00 2001 From: Lukas Riedel <34276446+peanutfun@users.noreply.github.com> Date: Thu, 13 Jun 2024 16:07:46 +0200 Subject: [PATCH 11/13] Format docstring of TropCyclone.from_tracks --- climada/hazard/trop_cyclone.py | 132 +++++++++++++++++---------------- 1 file changed, 69 insertions(+), 63 deletions(-) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index a5155ff066..acc8c6f67a 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -201,9 +201,9 @@ def from_tracks( 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 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. @@ -219,90 +219,96 @@ def from_tracks( tracks : climada.hazard.TCTracks Tracks of storm events. centroids : Centroids, optional - Centroids where to model TC. Default: global centroids at 360 arc-seconds resolution. + 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 + Pool that will be used for parallel computation of wind fields. Default: + None model : str, optional - Parametric wind field model to use: - "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, - or - "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. - Default: "H08". + 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. The Holland models support the following - parameters, in alphabetical order: + 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 + 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 + 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 + 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 + 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 - Default: None 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. + `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. + 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. + * "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 + 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 + 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 + 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 + 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 ------ @@ -1840,7 +1846,7 @@ def _stat_holland_2010( 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 From 48c44cca2a241533abe7407966386faf78eab464 Mon Sep 17 00:00:00 2001 From: Thomas Roosli Date: Fri, 21 Jun 2024 15:08:08 +0200 Subject: [PATCH 12/13] Issue debug message if Holland 2010 model is called with cyclostrophic=False --- climada/hazard/trop_cyclone.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index acc8c6f67a..b1bf7a72c6 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -1304,7 +1304,7 @@ def _compute_angular_windspeeds_h10( si_track: xr.Dataset, d_centr: np.ndarray, close_centr_msk: np.ndarray, - cyclostrophic: bool = False, + 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, @@ -1324,7 +1324,7 @@ def _compute_angular_windspeeds_h10( 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: False + 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. @@ -1348,6 +1348,11 @@ def _compute_angular_windspeeds_h10( ndarray of shape (npositions, ncentroids) containing the magnitude of the angular windspeed per track position per centroid location """ + 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) From ff8700d463a69b0e317de96d0e1800bb762b9e6c Mon Sep 17 00:00:00 2001 From: Lukas Riedel <34276446+peanutfun@users.noreply.github.com> Date: Fri, 21 Jun 2024 15:28:31 +0200 Subject: [PATCH 13/13] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f2814c97b9..b382ce2400 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ Code freeze date: YYYY-MM-DD CLIMADA tutorials. [#872](https://github.com/CLIMADA-project/climada_python/pull/872) - Centroids complete overhaul. Most function should be backward compatible. Internal data is stored in a geodataframe attribute. Raster are now stored as points, and the meta attribute is removed. Several methds were deprecated or removed. [#787](https://github.com/CLIMADA-project/climada_python/pull/787) - Improved error messages produced by `ImpactCalc.impact()` in case impact function in the exposures is not found in impf_set [#863](https://github.com/CLIMADA-project/climada_python/pull/863) +- Update the Holland et al. 2010 TC windfield model and introduce `model_kwargs` parameter to adjust model parameters [#846](https://github.com/CLIMADA-project/climada_python/pull/846) ### Fixed