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 <[email protected]> Co-authored-by: Thomas Vogt <[email protected]> Co-authored-by: Lukas Riedel <[email protected]>
CLIMADA-project · Jan 11, 2024 · 104dfdb · 104dfdb
1 parent e81249f
commit 104dfdb
Show file tree

Hide file tree

Showing 3 changed files with 71 additions and 39 deletions.
diff --git 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
@@ -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
@@ -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