Merge pull request #174 from nismod/feature/taper_advective_winds

Improve robustness of wind field estimation

config/config.yaml

@@ -112,7 +112,7 @@ storm_sets:
     IBTrACS_black-marble-validation: 'config/storm_sets/20230120_black_marble.json'
 
 # consider countries at risk of a storm if within this many degrees of any storm track point
-max_track_search_radius_deg: 2
+max_track_search_radius_deg: 3
 
 # wind speed is constant within a square area of sides wind_grid_resolution_deg
 # note that for large domains, e.g. USA or CHN, 0.05 deg resolution requires

src/open_gira/tests/test_wind.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from open_gira.wind import holland_wind_model, advective_vector, rotational_field
+from open_gira.wind import holland_wind_model, advective_vector, sigmoid_decay, estimate_wind_field
 
 
 class TestHollandWindModel:
@@ -84,35 +84,109 @@ def test_default_argument_regression(self):
         )
 
 
-class TestRotationalField:
+class TestSigmoidDecay:
+    """Test sigmoid/logistic decay function."""
+
+    def test_sigmoid_decay(self):
+        expected = [
+            0.98201379, 0.97340301, 0.96083428, 0.94267582, 0.9168273 ,
+            0.88079708, 0.83201839, 0.76852478, 0.68997448, 0.59868766,
+            0.5       , 0.40131234, 0.31002552, 0.23147522, 0.16798161,
+            0.11920292, 0.0831727 , 0.05732418, 0.03916572, 0.02659699,
+            0.01798621
+        ]
+
+        np.testing.assert_allclose(
+            expected,
+            sigmoid_decay(np.arange(0, 1001, 50), 500, 0.004),
+            rtol=1E-6
+        )
+
+
+class TestEstimateWindField:
     """
-    Test function which computes field of rotational vectors, that is intensity
-    from the Holland model, and rotation direction from hemisphere.
+    Test function which computes fields of advective and rotational vectors and combines them.
+
+    Testing a transect in longitude either side of storm eye, with single valued latitude.
     """
 
-    def test_rotational_field(self):
-        eye_lon = -30
-        lon_width = 5
-        max_lon = eye_lon + (lon_width - 1) / 2
-        min_lon = eye_lon - (lon_width - 1) / 2
-        lon_arr = np.linspace(min_lon, max_lon, lon_width)
-        eye_lat = 15
+    def test_estimate_wind_field_rotation_only(self):
+        # stationary, rotating with maximum velocity of 80ms-1 at 10,000m radius
+        # should rotate anticlockwise in northern hemisphere
+        result = estimate_wind_field(np.linspace(-30.5, -29.5, 9), 15, -30, 15, 10000, 80, 92000, 100000, 0, 0)
+        np.testing.assert_allclose(
+            np.abs(result),
+            np.array(
+                [[16.58826006, 23.81669229, 38.22875421, 72.35172612, np.nan,
+                72.35172612, 38.22875421, 23.81669229, 16.58826006]]
+            )
+        )
+        np.testing.assert_allclose(
+            np.angle(result, deg=True),
+            np.array(
+                # degrees CCW from positive real axis (mathematical convention)
+                [[-89.93529486, -89.95147127, -89.96764757, -89.9838238 ,
+                np.nan,  89.9838238 ,  89.96764757,  89.95147127, 89.93529486]]
+            )
+        )
 
-        expected = np.array(
-            [[0.00063 -0.139398j, 0.029922-13.247382j, np.nan + np.nan * 1j, 0.029922+13.247382j, 0.00063  +0.139398j]]
+    def test_estimate_wind_field_rotation_only_southern_hemisphere(self):
+        # stationary, rotating with maximum velocity of 80ms-1 at 10,000m radius
+        # should rotate clockwise in southern hemisphere
+        result = estimate_wind_field(np.linspace(-30.5, -29.5, 9), -15, -30, -15, 10000, 80, 92000, 100000, 0, 0)
+        np.testing.assert_allclose(
+            np.abs(result),
+            np.array(
+                [[16.58826006, 23.81669229, 38.22875421, 72.35172612, np.nan,
+                72.35172612, 38.22875421, 23.81669229, 16.58826006]]
+            )
         )
-        # test for regression
-        north = rotational_field(lon_arr, np.array([eye_lat]), eye_lon, eye_lat, 50_000, 100, 95_000, 100_000)
         np.testing.assert_allclose(
-            north,
-            expected,
-            rtol=1E-5
+            np.angle(result, deg=True),
+            np.array(
+                # flipped sign w.r.t. test_estimate_wind_field_rotation_only
+                # degrees anticlockwise from positive real axis (mathematical convention)
+                [[ 89.93529486,  89.95147127,  89.96764757,  89.9838238 , np.nan,
+                -89.9838238 , -89.96764757, -89.95147127, -89.93529486]]
+            )
         )
 
-        # check that the direction of rotation is reversed in the southern hemisphere
-        south = rotational_field(lon_arr, np.array([-eye_lat]), eye_lon, -eye_lat, 50_000, 100, 95_000, 100_000)
+    def test_estimate_wind_field_advection_only(self):
+        # storm moving N (heading 0 degrees) at 10ms-1, no pressure drop and no rotation (what a perverse storm)
+        result = estimate_wind_field(np.linspace(-30.5, -29.5, 9), 15, -30, 15, 10000, 1E-6, 99999, 100000, 0, 10)
+        # recovering alpha parameter, ~0.56 factor reduction from eye speed
+        np.testing.assert_allclose(
+            np.abs(result) / 10,
+            np.array(
+                [[0.54467011, 0.5461928 , 0.5475677 , 0.54880849, np.nan,
+                0.54880849, 0.5475677 , 0.5461928 , 0.54467011]]
+            )
+        )
+        # recovering beta parameter, 19.2 degrees over-rotation
+        np.testing.assert_allclose(
+            np.angle(result, deg=True) - 90,
+            np.array(
+                # degrees CCW from positive real axis (mathematical convention)
+                [[19.2, 19.2, 19.2, 19.2, np.nan, 19.2, 19.2, 19.2, 19.2]]
+            )
+        )
 
+    def test_estimate_wind_field_rotation_and_advection(self):
+        # storm moving N at 10ms-1, rotating at 80ms-1 with RMW=10km, 100mbar pressure drop
+        result = estimate_wind_field(np.linspace(-30.5, -29.5, 9), 15, -30, 15, 10000, 80, 90000, 100000, 0, 10)
+        # greater speeds on east of eye, where advection and rotation combine
         np.testing.assert_allclose(
-            np.angle(north),
-            -1 * np.angle(south)
+            np.abs(result),
+            np.array(
+                [[24.45427248, 31.92847355, 44.33622784, 65.62543488, np.nan,
+                75.9877587 , 54.67172047, 42.23278268, 34.72300576]]
+            )
         )
+        np.testing.assert_allclose(
+            np.angle(result, deg=True),
+            np.array(
+                # degrees CCW from positive real axis (mathematical convention)
+                [[-94.12223634, -93.16869142, -92.29164582, -91.55850813,
+                np.nan,  91.34593479,  91.8582487 ,  92.39504393, 92.90189219]]
+            )
+        )

src/open_gira/wind.py

@@ -167,11 +167,11 @@ def advective_vector(
     Geophys. Res., 117, D09120, doi:10.1029/2011JD017126
 
     Arguments:
-        eye_heading_deg (float): Heading of eye in degrees clockwise from north
-        eye_speed_ms (float): Speed of eye in metres per second
-        hemisphere (int): +1 for northern, -1 for southern
-        alpha (float): Fractional reduction of advective wind speed from eye speed
-        beta (float): Degrees advective winds tend to rotate past storm track (in direction of rotation)
+        eye_heading_deg: Heading of eye in degrees clockwise from north
+        eye_speed_ms: Speed of eye in metres per second
+        hemisphere: +1 for northern, -1 for southern
+        alpha: Fractional reduction of advective wind speed from eye speed
+        beta: Degrees advective winds tend to rotate past storm track (in direction of rotation)
 
     Returns:
         np.complex128: Advective wind vector
@@ -187,7 +187,20 @@ def advective_vector(
     return mag_v_a * np.sin(phi_a) + mag_v_a * np.cos(phi_a) * 1j
 
 
-def rotational_field(
+@numba.njit
+def sigmoid_decay(x: np.ndarray, midpoint: float, slope: float) -> np.ndarray:
+    """
+    Transform input by a sigmoid shape decay to return values in range [0, 1].
+
+    Args:
+        x: Input array to transform
+        midpoint: Decay midpoint in x
+        slope: Larger values decay faster
+    """
+    return 0.5 * (1 + np.tanh(slope * (midpoint - x)))
+
+
+def estimate_wind_field(
     longitude: np.ndarray,
     latitude: np.ndarray,
     eye_long: float,
@@ -196,24 +209,52 @@ def rotational_field(
     max_wind_speed_ms: float,
     min_pressure_pa: float,
     env_pressure_pa: float,
+    advection_azimuth_deg: float,
+    eye_speed_ms: float,
 ) -> np.ndarray:
     """
-    Calculate the rotational component of a storm's vector wind field
+    Given a spatial domain and tropical cyclone attributes, estimate a vector wind field.
+
+    The rotational component uses a modified Holland model. The advective
+    component (from the storm's translational motion) is modelled to fall off
+    from approximately 10 maximum wind radii. These two components are vector
+    added and returned.
 
     Args:
-        longitude (np.ndarray[float]): Grid values to evaluate on
-        latitude (np.ndarray[float]): Grid values to evaluate on
-        eye_long (float): Location of eye in degrees
-        eye_lat (float): Location of eye in degrees
-        radius_to_max_winds_m (float): Distance from eye centre to maximum wind speed in metres
-        max_wind_speed_ms (float): Maximum linear wind speed (relative to storm eye)
-        min_pressure_pa (float): Minimum pressure in storm eye in Pascals
-        env_pressure_pa (float): Environmental pressure, typical for this locale, in Pascals
+        longitude: Grid values to evaluate on
+        latitude: Grid values to evaluate on
+        eye_long: Location of eye in degrees
+        eye_lat: Location of eye in degrees
+        radius_to_max_winds_m: Distance from eye centre to maximum wind speed in metres
+        max_wind_speed_ms: Maximum wind speed (relative to ground)
+        min_pressure_pa: Minimum pressure in storm eye in Pascals
+        env_pressure_pa: Environmental pressure, typical for this locale, in Pascals
+        eye_heading_deg: Heading of eye in degrees clockwise from north
+        eye_speed_ms: Speed of eye in metres per second
 
     Returns:
-        np.ndarray[complex]: Grid of wind vectors
+        Grid of wind vectors
     """
 
+    # check inputs
+    assert 0 < max_wind_speed_ms < 130
+    assert 0 < radius_to_max_winds_m < 1500000
+    assert 75000 < min_pressure_pa < 102000
+    assert 0 <= eye_speed_ms < 40
+
+    # clip eye speed to a maximum of 30ms-1
+    # greater than this is non-physical, and probably the result of a data error
+    # we do not want to propagate such an error to our advective wind field
+    adv_vector: np.complex128 = advective_vector(
+        advection_azimuth_deg,
+        eye_speed_ms,
+        np.sign(eye_lat),
+    )
+
+    # maximum wind speed, less advective component
+    # this is the maximum tangential wind speed in the eye's non-rotating reference frame
+    max_wind_speed_relative_to_eye_ms: float = max_wind_speed_ms - np.abs(adv_vector)
+
     X, Y = np.meshgrid(longitude, latitude)
     grid_shape = X.shape  # or Y.shape
 
@@ -225,21 +266,29 @@ def rotational_field(
         np.full(len(Y.ravel()), eye_lat),
     )
 
+    distance_to_eye_grid_m = radius_m.reshape(grid_shape)
+
+    # decay effect of advective field from maximum at storm eye out to zero at 1,000km radius
+    # we shouldn't claim any authority on winds outside the vicinity of the storm
+    adv_field: np.ndarray = adv_vector * sigmoid_decay(distance_to_eye_grid_m / 1_000, 500, 0.004)
+
     # magnitude of rotational wind component
     mag_v_r: np.ndarray = holland_wind_model(
         radius_to_max_winds_m,
-        max_wind_speed_ms,
+        max_wind_speed_relative_to_eye_ms,
         min_pressure_pa,
         env_pressure_pa,
-        radius_m.reshape(grid_shape),
+        distance_to_eye_grid_m,
         eye_lat
     )
 
     # azimuth of rotational component is tangent to radius, with direction set by hemisphere
     phi_r: np.ndarray = np.radians(grid_to_eye_azimuth_deg.reshape(grid_shape) + np.sign(eye_lat) * 90)
 
     # find components of vector at each pixel
-    return mag_v_r * np.sin(phi_r) + mag_v_r * np.cos(phi_r) * 1j
+    rot_field = mag_v_r * np.sin(phi_r) + mag_v_r * np.cos(phi_r) * 1j
+
+    return adv_field + rot_field
 
 
 def interpolate_track(track: gpd.GeoDataFrame, frequency: str = "1H") -> gpd.GeoDataFrame:
@@ -257,7 +306,6 @@ def interpolate_track(track: gpd.GeoDataFrame, frequency: str = "1H") -> gpd.Geo
         gpd.GeoDataFrame: Track with min_pressure_hpa, max_wind_speed_ms,
             radius_to_max_winds_km and geometry columns interpolated.
     """
-
     if len(track) == 0:
         raise ValueError("No track data")
     elif len(track) == 1:
@@ -294,8 +342,19 @@ def interpolate_track(track: gpd.GeoDataFrame, frequency: str = "1H") -> gpd.Geo
     # merge dataframes, on index collision, keep the non-NaN values from track
     interp_track = track.combine_first(interp_domain).sort_index()
 
+    # don't interpolate over some sensible duration
+    # note that we limit by an integer number of timesteps, not a time
+    # so our implicit assumption is that the input index is equally spaced
+    max_steps_to_fill: int = np.round(pd.Timedelta("6H") / pd.Timedelta(frequency)).astype(int)
+
     # interpolate over numeric value of index
-    interp_track.loc[:, interp_cols] = interp_track.loc[:, interp_cols].interpolate(method=interp_method)
+    interp_track.loc[:, interp_cols] = interp_track.loc[:, interp_cols].interpolate(
+        method=interp_method,
+        limit=max_steps_to_fill
+    )
+
+    # fail if we still have NaN values (probably the time gaps exceed `max_steps_to_fill`
+    assert not interp_track.loc[:, interp_cols].isnull().values.any()
 
     interp_track["geometry"] = gpd.points_from_xy(
         interp_track.x, interp_track.y, crs="EPSG:4326"