diff --git a/config/config.yaml b/config/config.yaml
index f11d1213..f6680929 100644
--- a/config/config.yaml
+++ b/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
diff --git a/src/open_gira/tests/test_wind.py b/src/open_gira/tests/test_wind.py
index 01d8d13e..17729625 100644
--- a/src/open_gira/tests/test_wind.py
+++ b/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]]
+            )
+        )
\ No newline at end of file
diff --git a/src/open_gira/wind.py b/src/open_gira/wind.py
index 6c1659e5..21a34ef6 100644
--- a/src/open_gira/wind.py
+++ b/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,13 +266,19 @@ 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
     )
 
@@ -239,7 +286,9 @@ def rotational_field(
     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"
diff --git a/src/open_gira/wind_plotting.py b/src/open_gira/wind_plotting.py
index 7d3f196f..423347dc 100644
--- a/src/open_gira/wind_plotting.py
+++ b/src/open_gira/wind_plotting.py
@@ -4,21 +4,28 @@
 
 
 import geopandas as gpd
+import matplotlib
 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
 
 from open_gira.curves import logistic_min
 
+matplotlib.use("Agg")
+plt.set_loglevel("warning")
 
-WIND_CMAP = "turbo"
-MAX_SPEED = 80  # clip wind speeds above this value when plotting
-# origin lower so latitude indicies increasing northwards
-WIND_PLOT_ORIGIN = "lower"
+WIND_CMAP = "magma_r"
+WIND_PLOT_ORIGIN = "lower"  # origin lower so latitude indicies increasing northwards
 QUIVER_SCALE = 2_500  # larger values give shorter quiver arrows
 
+# discretise the wind field chloropleth maps
+MIN = 3
+MAX = 75
+STEP = 3
+LEVELS = int((MAX - MIN) / STEP + 1)
+
 
 def size_plot(i: int, j: int) -> tuple[float, float]:
     """
@@ -42,7 +49,7 @@ def plot_quivers(field: np.ndarray, title: str, colorbar_label: str, file_path:
     ax.quiver(field.real, field.imag, angles='xy', scale=QUIVER_SCALE, color='white')
 
     mag = np.abs(field)
-    img = ax.imshow(mag, vmin=0, vmax=MAX_SPEED, origin=WIND_PLOT_ORIGIN, cmap=WIND_CMAP)
+    img = ax.imshow(mag, vmin=0, vmax=80, origin=WIND_PLOT_ORIGIN, cmap=WIND_CMAP)
     fig.colorbar(img, ax=ax, location="right", label=colorbar_label, shrink=0.81)
 
     stats_str = fr"min: {mag.min():.2f}, max: {mag.max():.2f}, $\sigma:$ {mag.std():.2f}"
@@ -67,7 +74,19 @@ 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=LEVELS,
+        x="i",
+        y="j",
+        ax=ax,
+        vmin=MIN,
+        vmax=MAX,
+        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 +128,31 @@ 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(MIN, MAX, LEVELS), cmap.N),
+        origin=WIND_PLOT_ORIGIN,
+        cmap=cmap,
     )
 
-    def animate(i, image, quiver):
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size="5%", pad=0.1)
+    fig.colorbar(img, ax=ax, cax=cax, label="Wind speed [m/s]", extend="min")
+
+    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,
     )
 
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc
index 5813a3df..1e4892fa 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png
index 8bf1aca7..c70aca47 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png
index 91df67e7..653a6454 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png
index 3ef338a1..0eb93695 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png
index 5340fa20..ef6b29ea 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png
index 0a223a94..66440caa 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png
index dfb0683f..9abe0736 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png
index 7e2e6b27..d798497d 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png
index 28b41f15..55015e60 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png
index 3b9710c4..d4a8517b 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png
index 5bc8f760..e5b5070b 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png
index 2b09739c..6cc79dc3 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png
index 0378dbf0..67e63482 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png
index 02d1436a..1dc1d5ac 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png
index 46243093..263b2e8e 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png
index 6fdd1023..aa5e2a11 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png
index a07c272f..408511e1 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png
index a7e9fb24..cd7226a6 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png
index 085c7ce8..fd27174f 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png
index 33532afa..68297077 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png
index fd5ba576..f36ec485 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png
index 1e24fd2b..fad3329b 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png
index 0cd5f1a6..765f3e6d 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png
index 8364e322..ac3bce5d 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png
index 36df4d14..e47f5e01 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png
index 4ef24c66..4032a117 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png
index 5c187829..26a0f872 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png differ
diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png
index 07744eb9..0039a601 100644
Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png and b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png differ
diff --git a/workflow/scripts/intersect/estimate_wind_fields.py b/workflow/scripts/intersect/estimate_wind_fields.py
index be760bb9..fabc84d1 100644
--- a/workflow/scripts/intersect/estimate_wind_fields.py
+++ b/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
 
@@ -80,11 +80,14 @@ def process_track(
     basin: str = track.iloc[0, track.columns.get_loc("basin_id")]
 
     # interpolate track (avoid 'doughnut effect' of wind field from infrequent eye observations)
-    track: gpd.GeoDataFrame = interpolate_track(track)
-
-    geod_wgs84: pyproj.Geod = pyproj.CRS("epsg:4326").get_geod()
+    try:
+        track: gpd.GeoDataFrame = interpolate_track(track)
+    except AssertionError:
+        logging.warning(f"Could not successfully interpolate {track_id}")
+        return track_id, np.zeros_like(downscaling_factors)
 
     # forward azimuth angle and distances from track eye to next track eye
+    geod_wgs84: pyproj.Geod = pyproj.CRS("epsg:4326").get_geod()
     advection_azimuth_deg, _, eye_step_distance_m = geod_wgs84.inv(
         track.geometry.x.iloc[:-1],
         track.geometry.y.iloc[:-1],
@@ -101,39 +104,26 @@ def process_track(
     # calculate eye speed
     track["eye_speed_ms"] = eye_step_distance_m / period.seconds.values
 
-    # 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(
-            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.min_pressure_hpa * 100,  # convert to Pascals
-            ENV_PRESSURE[basin] * 100,  # convert to Pascals
-        )
-
-    # sum components of wind field, (timesteps, y, x)
-    wind_field: np.ndarray[complex] = adv_field + rot_field
+        try:
+            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
+                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,
+            )
+        except AssertionError:
+            logging.warning(f"{track_id} failed wind field estimation for {track_i + 1} of {len(track)}, writing zeros")
 
     # take factors calculated from surface roughness of region and use to downscale speeds
     downscaled_wind_field = downscaling_factors * wind_field
diff --git a/workflow/scripts/preprocess/parse_IBTrACS.py b/workflow/scripts/preprocess/parse_IBTrACS.py
index f6ec53ec..4bfa4749 100644
--- a/workflow/scripts/preprocess/parse_IBTrACS.py
+++ b/workflow/scripts/preprocess/parse_IBTrACS.py
@@ -172,12 +172,6 @@ def saffir_simpson_classifier(wind_speed_ms: float) -> Union[int, float]:
     # change geometry from 0-360 to -180-180
     df.lon = np.where(df.lon > 180, df.lon - 360, df.lon)
 
-    df["point_id"] = df.apply(lambda row: f"{row.track_id}_{row.lat}_{row.lon}", axis=1)
-    n_rows_raw = len(df)
-    logging.info(f"Collated {n_rows_raw} track points")
-    df = df.drop_duplicates(subset="point_id").drop(columns=["point_id"])
-    logging.info(f"Dropped {n_rows_raw - len(df)} track points as duplicates")
-
     # construct geometry from lat and long
     df = gpd.GeoDataFrame(
         data=df,
diff --git a/workflow/scripts/preprocess/parse_IRIS.py b/workflow/scripts/preprocess/parse_IRIS.py
index 2b1f7401..2f4225eb 100644
--- a/workflow/scripts/preprocess/parse_IRIS.py
+++ b/workflow/scripts/preprocess/parse_IRIS.py
@@ -46,7 +46,9 @@
     sample = snakemake.wildcards.SAMPLE
 
     data = []
-    for path in tqdm(natural_sort(glob(f"{csv_dir}/*_1000Y_n{sample}.txt"))):
+    for path in natural_sort(glob(f"{csv_dir}/*_1000Y_n{sample}.txt")):
+
+        logging.info(path)
 
         df = pd.read_csv(path, header=1, sep=r"\s+", names=IRIS_CSV_SCHEMA.keys(), dtype=IRIS_CSV_SCHEMA)
         df = df.rename(
@@ -82,6 +84,12 @@
             + df["tc_number"].astype(int).astype(str)
         )
 
+        # drop any duplicates
+        df["point_id"] = df.apply(lambda row: f"{row.track_id}_{row.timestep}", axis=1)
+        n_rows_raw = len(df)
+        df = df.drop_duplicates(subset="point_id").drop(columns=["point_id"])
+        logging.info(f"Collated {n_rows_raw} track points, dropped {n_rows_raw - len(df)} as duplicates")
+
         # we'll want to interpolate and then measure the speed of tracks later,
         # this is easiest when we have some temporal index (as in IBTrACS)
         # so make up an artificial one here based on the IRIS reporting frequency
@@ -97,12 +105,6 @@
 
     df = pd.concat(data)
 
-    df["point_id"] = df.apply(lambda row: f"{row.track_id}_{row.lat}_{row.lon}", axis=1)
-    n_rows_raw = len(df)
-    logging.info(f"Collated {n_rows_raw} track points")
-    df = df.drop_duplicates(subset="point_id").drop(columns=["point_id"])
-    logging.info(f"Dropped {n_rows_raw - len(df)} track points as duplicates")
-
     # construct geometry from lat and long
     df = gpd.GeoDataFrame(
         data=df,
diff --git a/workflow/scripts/preprocess/parse_STORM.py b/workflow/scripts/preprocess/parse_STORM.py
index 65d0960e..f8ea2b3d 100644
--- a/workflow/scripts/preprocess/parse_STORM.py
+++ b/workflow/scripts/preprocess/parse_STORM.py
@@ -55,7 +55,9 @@
 
     data = []
     # loop over basins for given sample, processing tracks into a common format
-    for path in tqdm(natural_sort(glob(f"{csv_dir}/*_1000_YEARS_{sample}*.csv"))):
+    for path in natural_sort(glob(f"{csv_dir}/*_1000_YEARS_{sample}*.csv")):
+
+        logging.info(path)
 
         df = pd.read_csv(path, names=STORM_CSV_SCHEMA.keys(), dtype=STORM_CSV_SCHEMA)
 
@@ -79,10 +81,15 @@
             + df["tc_number"].astype(int).astype(str)
         )
 
+        # STORM contains duplicate track points (duplicate timesteps for a given year-tc_number combination)
+        df["point_id"] = df.apply(lambda row: f"{row.track_id}_{row.timestep}", axis=1)
+        n_rows_raw = len(df)
+        df = df.drop_duplicates(subset="point_id").drop(columns=["point_id"])
+        logging.info(f"Collated {n_rows_raw} track points, dropped {n_rows_raw - len(df)} as duplicates")
+
         # we'll want to interpolate and then measure the speed of tracks later,
         # this is easiest when we have some temporal index (as in IBTrACS)
         # so make up an artificial one here based on the STORM reporting frequency
-
         track_datetimes: List[np.ndarray] = []
         track_lengths: np.ndarray = df.track_id.apply(hash).value_counts(sort=False).values
         for length in track_lengths:
@@ -100,12 +107,6 @@
     # rescale winds to 1-minutely
     df.max_wind_speed_ms /= STORM_1MIN_WIND_FACTOR
 
-    df["point_id"] = df.apply(lambda row: f"{row.track_id}_{row.lat}_{row.lon}", axis=1)
-    n_rows_raw = len(df)
-    logging.info(f"Collated {n_rows_raw} track points")
-    df = df.drop_duplicates(subset="point_id").drop(columns=["point_id"])
-    logging.info(f"Dropped {n_rows_raw - len(df)} track points as duplicates")
-
     # construct geometry from lat and long
     df = gpd.GeoDataFrame(
         data=df,
diff --git a/workflow/scripts/preprocess/slice_storm_tracks.py b/workflow/scripts/preprocess/slice_storm_tracks.py
index 0b3514d5..c8ea06cf 100644
--- a/workflow/scripts/preprocess/slice_storm_tracks.py
+++ b/workflow/scripts/preprocess/slice_storm_tracks.py
@@ -1,5 +1,6 @@
 """
-Subset IBTrACS (historic) or STORM (synthetic) storm tracks by a buffer polygon
+Subset storm tracks by an AOI buffer polygon, taking all points in the track from
+first arrival to last departure (including voyages outside the AOI).
 """
 
 import os
@@ -8,7 +9,10 @@
 import sys
 
 import geopandas as gpd
+import numpy as np
+import pandas as pd
 import shapely
+from tqdm import tqdm
 
 
 logging.basicConfig(format="%(asctime)s %(process)d %(filename)s %(message)s", level=logging.INFO)
@@ -41,8 +45,29 @@
     logging.info("Subsetting tracks to those which intersect with grid (+ buffer)")
     # adding the buffer reduces the likelihood we construct wind fields where a
     # storm suddenly appears well inside the domain
-    tracks = tracks[tracks.intersects(hull.buffer(track_slicing_buffer_deg))]
+    AOI_points = tracks[tracks.intersects(hull.buffer(track_slicing_buffer_deg))]
+    logging.info(f"Found {len(AOI_points.track_id.unique())} tracks passing within buffer")
+
+    # some storms impact a country, wander off outside our AOI, then return some days later
+    # want time continuity in our tracks (so we don't interpolate wildly, so our eye speed estimates are reasonable)
+    logging.info("Indexing tracks' first arrival to last departure (for time continuity)")
+    sliced_tracks_by_track: list[pd.DataFrame] = []
+    for track_id in tqdm(AOI_points.track_id.unique()):
+        df = AOI_points[AOI_points.track_id == track_id]
+
+        # we need at least two track points to calculate the eye velocity and advective winds
+        if len(df) > 1:
+            arrival, *_, departure = df.index
+            sliced = tracks[tracks.track_id == track_id].loc[arrival: departure]
+            assert sliced.timestep.is_monotonic_increasing
+            assert np.all(np.diff(sliced.timestep) == 1)
+            sliced_tracks_by_track.append(sliced)
+
+    sliced_tracks = pd.concat(sliced_tracks_by_track)
+    logging.info(
+        f"Filtered to {len(sliced_tracks.track_id.unique())} valid (n > 1) tracks, with {len(sliced_tracks)} points"
+    )
 
     logging.info("Writing tracks subset to disk")
     os.makedirs(os.path.dirname(sliced_tracks_path), exist_ok=True)
-    tracks.to_parquet(sliced_tracks_path)
+    sliced_tracks.to_parquet(sliced_tracks_path)