From 104dfdb64ed00fce6c49a5d67be371f4fca016dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20R=C3=B6=C3=B6sli?= Date: Thu, 11 Jan 2024 10:48:36 +0100 Subject: [PATCH] Feature/tc holland vtrans (#833) * Change how vtrans is considered in the implementation of Holland 2008 and Holland 2010 model * trop_cyclone: broadcast_arrays -> broadcast_to * Update CHANGELOG.md --------- Co-authored-by: Thomas Roosli Co-authored-by: Thomas Vogt <57705593+tovogt@users.noreply.github.com> Co-authored-by: Lukas Riedel <34276446+peanutfun@users.noreply.github.com> --- CHANGELOG.md | 1 + climada/hazard/test/test_trop_cyclone.py | 33 ++++++---- climada/hazard/trop_cyclone.py | 76 ++++++++++++++++-------- 3 files changed, 71 insertions(+), 39 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c89d56c8ea..6183d3eae5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ Code freeze date: YYYY-MM-DD - `Hazard.from_xarray_raster` now stores strings as default values for `Hazard.event_name` [#795](https://github.com/CLIMADA-project/climada_python/pull/795) - Fix the dist_approx util function when used with method="geosphere" and log=True and points that are very close. [#792](https://github.com/CLIMADA-project/climada_python/pull/792) - `climada.util.yearsets.sample_from_poisson`: fix a bug ([#819](https://github.com/CLIMADA-project/climada_python/issues/819)) and inconsistency that occurs when lambda events per year (`lam`) are set to 1. [[#823](https://github.com/CLIMADA-project/climada_python/pull/823)] +- In the TropCyclone class in the Holland model 2008 and 2010 implementation, a doublecounting of translational velocity is removed [#833](https://github.com/CLIMADA-project/climada_python/pull/833) ### Deprecated diff --git a/climada/hazard/test/test_trop_cyclone.py b/climada/hazard/test/test_trop_cyclone.py index a52fedc5e5..ffa1ff9e36 100644 --- a/climada/hazard/test/test_trop_cyclone.py +++ b/climada/hazard/test/test_trop_cyclone.py @@ -60,9 +60,10 @@ def test_memory_limit(self): tc_haz = TropCyclone.from_tracks(tc_track, centroids=CENTR_TEST_BRB, max_memory_gb=0.001) intensity_idx = [0, 1, 2, 3, 80, 100, 120, 200, 220, 250, 260, 295] intensity_values = [ - 25.60778909, 26.90887264, 28.26624642, 25.54092386, 31.21941738, 36.16596567, - 21.11399856, 28.01452136, 32.65076804, 31.33884098, 0, 40.27002104, + 22.74903, 23.784691, 24.82255, 22.67403, 27.218706, 30.593959, + 18.980878, 24.540069, 27.826407, 26.846293, 0., 34.568898, ] + np.testing.assert_array_almost_equal( tc_haz.intensity[0, intensity_idx].toarray()[0], intensity_values, @@ -72,12 +73,14 @@ def test_set_one_pass(self): """Test _tc_from_track function.""" intensity_idx = [0, 1, 2, 3, 80, 100, 120, 200, 220, 250, 260, 295] intensity_values = { - "geosphere": [25.60794285, 26.90906280, 28.26649026, 25.54076797, 31.21986961, - 36.17171808, 21.11408573, 28.01457948, 32.65349378, 31.34027741, 0, - 40.27362679], - "equirect": [25.60778909, 26.90887264, 28.26624642, 25.54092386, 31.21941738, - 36.16596567, 21.11399856, 28.01452136, 32.65076804, 31.33884098, 0, - 40.27002104] + "geosphere": [ + 22.74927, 23.78498, 24.822908, 22.674202, 27.220042, 30.602122, + 18.981022, 24.540138, 27.830925, 26.8489, 0., 34.572391, + ], + "equirect": [ + 22.74903, 23.784691, 24.82255, 22.67403, 27.218706, 30.593959, + 18.980878, 24.540069, 27.826407, 26.846293, 0., 34.568898, + ] } # the values for the two metrics should agree up to first digit at least for i, val in enumerate(intensity_values["geosphere"]): @@ -108,7 +111,7 @@ def test_set_one_pass(self): self.assertTrue(isinstance(tc_haz.intensity, sparse.csr_matrix)) self.assertEqual(tc_haz.intensity.shape, (1, 296)) - self.assertEqual(np.nonzero(tc_haz.intensity)[0].size, 280) + self.assertEqual(np.nonzero(tc_haz.intensity)[0].size, 255) np.testing.assert_array_almost_equal( tc_haz.intensity[0, intensity_idx].toarray()[0], intensity_values[metric]) @@ -127,10 +130,14 @@ 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": [25.60778909, 26.90887264, 28.26624642, 25.54092386, 31.21941738, 36.16596567, - 21.11399856, 28.01452136, 32.65076804, 31.33884098, 0, 40.27002104], - "H10": [27.604317, 28.720708, 29.894993, 27.52234 , 32.512395, 37.114355, - 23.848917, 29.614752, 33.775593, 32.545347, 19.957627, 41.014578], + "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 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, diff --git a/climada/hazard/trop_cyclone.py b/climada/hazard/trop_cyclone.py index bdd1b98977..125063012b 100644 --- a/climada/hazard/trop_cyclone.py +++ b/climada/hazard/trop_cyclone.py @@ -955,11 +955,6 @@ def _compute_windfields( si_track, d_centr, close_centr_msk, model, cyclostrophic=False, ) - # vectorial angular velocity - windfields = ( - si_track.attrs["latsign"] * np.array([1.0, -1.0])[..., :] * v_centr_normed[:, :, ::-1] - ) - windfields[close_centr_msk] *= v_ang_norm[close_centr_msk, None] # Influence of translational speed decreases with distance from eye. # The "absorbing factor" is according to the following paper (see Fig. 7): @@ -969,11 +964,26 @@ def _compute_windfields( # wind speed profiles in a GIS. UNED/GRID-Geneva. # https://unepgrid.ch/en/resource/19B7D302 # - t_rad_bc = np.broadcast_arrays(si_track["rad"].values[:, None], d_centr)[0] + 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[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']]: + # 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[close_centr_msk] -= ( + vtrans_norm_bc[close_centr_msk] * v_trans_corr[close_centr_msk] + ) + + # vectorial angular velocity + windfields = ( + si_track.attrs["latsign"] * np.array([1.0, -1.0])[..., :] * v_centr_normed[:, :, ::-1] + ) + windfields[close_centr_msk] *= v_ang_norm[close_centr_msk, None] + # 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 @@ -1407,11 +1417,15 @@ def _x_holland_2010( """ hol_x = np.zeros_like(d_centr) r_max, v_max_s, hol_b, d_centr, v_n, r_n = [ - ar[close_centr] for ar in np.broadcast_arrays( - 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], - ) + np.broadcast_to(ar, d_centr.shape)[close_centr] + 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 @@ -1470,11 +1484,15 @@ def _stat_holland_2010( Absolute values of wind speeds (in m/s) in angular direction. """ v_ang = np.zeros_like(d_centr) - d_centr, v_max_s, r_max, hol_b, hol_x = [ - ar[close_centr] for ar in np.broadcast_arrays( - d_centr, si_track["vmax"].values[:, None], si_track["rad"].values[:, None], - si_track["hol_b"].values[:, None], hol_x, - ) + v_max_s, r_max, hol_b, d_centr, hol_x = [ + np.broadcast_to(ar, d_centr.shape)[close_centr] + 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 @@ -1526,12 +1544,16 @@ def _stat_holland_1980( Absolute values of wind speeds (m/s) in angular direction. """ v_ang = np.zeros_like(d_centr) - d_centr, r_max, hol_b, penv, pcen, coriolis_p = [ - ar[close_centr] for ar in np.broadcast_arrays( - d_centr, si_track["rad"].values[:, None], si_track["hol_b"].values[:, None], - si_track["env"].values[:, None], si_track["cen"].values[:, None], - si_track["cp"].values[:, None] - ) + r_max, hol_b, penv, pcen, coriolis_p, 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["cp"].values[:, None], + d_centr, + ] ] r_coriolis = 0 @@ -1587,12 +1609,14 @@ def _stat_er_2011( Absolute values of wind speeds (m/s) in angular direction. """ v_ang = np.zeros_like(d_centr) - d_centr, r_max, v_max, coriolis_p = [ - ar[close_centr] for ar in np.broadcast_arrays( - d_centr, si_track["rad"].values[:, None], + r_max, v_max, coriolis_p, d_centr = [ + np.broadcast_to(ar, d_centr.shape)[close_centr] + 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