Taper advective wind field as function of radius to eye

The estimated wind field for each time step is a vector sum of two
components, the advective field due to the translational motion of the
storm and the rotational field. The rotational field intensity is a
function of radius from storm eye. Previously, the advective field
had no spatial dependence. This meant the entire spatial domain had some
sort of background wind field (typically(!) a few m/s). In the event of
a fast moving storm (say a decayed cyclone transitting the NA, this
speed may be quite significant. For large domains like the USA, this
could mean not insignificant background wind speeds many thousands of
miles from the storm. I have decided that preventing this is worth the
necessary assumptions and extra complexity.

We now take the advective wind vector and multiply it by some reduction
factor. This factor is the distance to storm eye, normalised by the
radius to maximum winds, then tapered to 0 using a hyperbolic tangent
function.

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, cutoff: float, steepness: float) -> np.ndarray:
+    """
+    Transform input by a sigmoid shape decay to return values in range [0, 1].
+
+    Args:
+        x: Input array to transform
+        cutoff: Approximate location in x where decay starts
+        steepness: How quickly the function decays
+    """
+    return 0.5 * (1 + np.tanh((cutoff / 2) - x / steepness))
+
+
+def estimate_wind_field(
     longitude: np.ndarray,
     latitude: np.ndarray,
     eye_long: float,
@@ -196,23 +209,43 @@ def rotational_field(
     max_wind_speed_ms: float,
     min_pressure_pa: float,
     env_pressure_pa: float,
+    advection_azimuth_deg: float,
+    eye_speed_ms: float,
+    hemisphere: int,
 ) -> 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
+        hemisphere: +1 for northern, -1 for southern
 
     Returns:
-        np.ndarray[complex]: Grid of wind vectors
+        Grid of wind vectors
     """
+    adv_vector: np.complex128 = advective_vector(
+        advection_azimuth_deg,
+        eye_speed_ms,
+        hemisphere,
+    )
+
+    # 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 +258,35 @@ def rotational_field(
         np.full(len(Y.ravel()), eye_lat),
     )
 
+    distance_to_eye_grid_m = radius_m.reshape(grid_shape)
+
+    # find the radius, r to eye for each pixel, normalised by radius to maximum winds, RMW
+    # multiply by a decay function, so that
+    # for r / RMW < 5, output ~ 1
+    # for 10 < r / RMW < 30, there is decay from 1 to 0
+    # for r / RMW > 30, output ~ 0
+
+    # including this decay hopefully prevents us from making nonsense estimates
+    # of wind speeds thousands of km from the storm eye
+    adv_field: np.ndarray = adv_vector * sigmoid_decay(distance_to_eye_grid_m / radius_to_max_winds_m, 6, 6)
+
     # 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:

src/open_gira/wind_plotting.py

@@ -5,7 +5,7 @@
 
 import geopandas as gpd
 import matplotlib.pyplot as plt
-from matplotlib import animation
+from matplotlib import animation, colors, colormaps
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 import numpy as np
 import xarray as xr
@@ -67,7 +67,9 @@ def plot_contours(field: np.ndarray, title: str, colorbar_label: str, file_path:
 
     da = xr.DataArray(field.T, coords={"i": range(x), "j": range(y)})
 
-    xr.plot.pcolormesh(da, levels=17, x="i", y="j", ax=ax, vmin=0, vmax=MAX_SPEED, cmap=WIND_CMAP, cbar_ax=cax)
+    cmap = colormaps[WIND_CMAP]
+    cmap.set_under("w")
+    xr.plot.pcolormesh(da, levels=16, x="i", y="j", ax=ax, vmin=5, vmax=MAX_SPEED, cmap=cmap, cbar_ax=cax)
 
     stats_str = fr"min: {field.min():.2f}, mean: {field.mean():.2f}, max: {field.max():.2f}"
     ax.set_title(title + "\n" + stats_str)
@@ -109,29 +111,28 @@ def animate_track(wind_field: np.ndarray, track: gpd.GeoDataFrame, file_path: st
     track_length, *grid_shape = wind_field.shape
 
     fig, ax = plt.subplots(figsize=size_plot(*grid_shape[::-1]))
+    cmap = colormaps[WIND_CMAP]
+    cmap.set_under("w")
 
     # origin lower so latitude indicies increasing northwards
-    img = ax.imshow(np.zeros(grid_shape), vmin=0, vmax=MAX_SPEED, origin=WIND_PLOT_ORIGIN, cmap=WIND_CMAP)
-    fig.colorbar(img, ax=ax, location="right", label="Wind speed [m/s]", shrink=0.81)
-    quiv = ax.quiver(
-        np.zeros(grid_shape),
+    img = ax.imshow(
         np.zeros(grid_shape),
-        angles='xy',
-        scale=QUIVER_SCALE,
-        color='white'
+        norm=colors.BoundaryNorm(np.linspace(5, 80, 16), cmap.N, extend="both"),
+        origin=WIND_PLOT_ORIGIN,
+        cmap=cmap,
     )
+    fig.colorbar(img, ax=ax, location="right", label="Wind speed [m/s]", shrink=0.81)
 
-    def animate(i, image, quiver):
+    def animate(i, image):
         image.set_array(np.abs(wind_field[i]))
-        quiver.set_UVC(wind_field[i].real, wind_field[i].imag)
         ax.set_title(f"{track_name} wind vector\n{i + 1} of {track_length}")
-        return img, quiv
+        return [img]
 
     anim = animation.FuncAnimation(
         fig,
         animate,
         frames=track_length,
-        fargs=(img, quiv),
+        fargs=(img,),
         blit=True,
     )
 

workflow/scripts/intersect/estimate_wind_fields.py

@@ -18,8 +18,8 @@
 
 from open_gira.io import bit_pack_dataarray_encoding
 from open_gira.wind import (
-    advective_vector, rotational_field, interpolate_track,
-    empty_wind_da, WIND_COORDS, ENV_PRESSURE
+    estimate_wind_field, interpolate_track, empty_wind_da, WIND_COORDS,
+    ENV_PRESSURE
 )
 from open_gira.wind_plotting import plot_contours, animate_track
 
@@ -104,37 +104,25 @@ def process_track(
     # hemisphere belongs to {-1, 1}
     track["hemisphere"] = np.sign(track.geometry.y)
 
-    adv_field: np.ndarray = np.zeros((len(track), *grid_shape), dtype=complex)
-    rot_field: np.ndarray = np.zeros((len(track), *grid_shape), dtype=complex)
+    # result array
+    wind_field: np.ndarray = np.zeros((len(track), *grid_shape), dtype=complex)
 
     for track_i, track_point in enumerate(track.itertuples()):
 
-        adv_vector: np.complex128 = advective_vector(
-            track_point.advection_azimuth_deg,
-            track_point.eye_speed_ms,
-            track_point.hemisphere,
-        )
-
-        adv_field[track_i, :] = np.full(grid_shape, adv_vector)
-
-        # 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 = track_point.max_wind_speed_ms - np.abs(adv_vector)
-
-        rot_field[track_i, :] = rotational_field(
+        wind_field[track_i, :] = estimate_wind_field(
             longitude,  # degrees
             latitude,  # degrees
             track_point.geometry.x,  # degrees
             track_point.geometry.y,  # degrees
             track_point.radius_to_max_winds_km * 1_000,  # convert to meters
-            max_wind_speed_relative_to_eye_ms,
+            track_point.max_wind_speed_ms,
             track_point.min_pressure_hpa * 100,  # convert to Pascals
             ENV_PRESSURE[basin] * 100,  # convert to Pascals
+            track_point.advection_azimuth_deg,
+            track_point.eye_speed_ms,
+            track_point.hemisphere,  # {+1|-1}
         )
 
-    # sum components of wind field, (timesteps, y, x)
-    wind_field: np.ndarray[complex] = adv_field + rot_field
-
     # take factors calculated from surface roughness of region and use to downscale speeds
     downscaled_wind_field = downscaling_factors * wind_field