Fix ship/underway instrument sampling trajectories

src/virtualship/expedition/simulate_schedule.py

@@ -116,7 +116,7 @@ def __init__(self, projection: pyproj.Geod, expedition: Expedition) -> None:
         self._next_ship_underwater_st_time = self._time
 
     def simulate(self) -> ScheduleOk | ScheduleProblem:
-        # TODO: instrument config mapping (as introduced in #269) should be helpful for refactoring here...
+        # TODO: instrument config mapping (as introduced in #269) should be helpful for refactoring here (i.e. #236)...
 
         for wp_i, waypoint in enumerate(self._expedition.schedule.waypoints):
             # sail towards waypoint
@@ -140,100 +140,116 @@ def simulate(self) -> ScheduleOk | ScheduleProblem:
 
             # wait while measurements are being done
             self._progress_time_stationary(time_passed)
+
         return ScheduleOk(self._time, self._measurements_to_simulate)
 
     def _progress_time_traveling_towards(self, location: Location) -> None:
+        """Travel from current location/waypoint to next waypoint, also mark locations and times for underway instrument measurements."""
         time_to_reach, azimuth1, ship_speed_meter_per_second = _calc_sail_time(
             self._location,
             location,
             self._expedition.ship_config.ship_speed_knots,
             self._projection,
         )
         end_time = self._time + time_to_reach
+        distance_to_move = ship_speed_meter_per_second * time_to_reach.total_seconds()
+
         # note all ADCP measurements
         if self._expedition.instruments_config.adcp_config is not None:
-            location = self._location
-            time = self._time
-            while self._next_adcp_time <= end_time:
-                time_to_sail = self._next_adcp_time - time
-                distance_to_move = (
-                    ship_speed_meter_per_second * time_to_sail.total_seconds()
-                )
-                geodfwd: tuple[float, float, float] = self._projection.fwd(
-                    lons=location.lon,
-                    lats=location.lat,
-                    az=azimuth1,
-                    dist=distance_to_move,
-                )
-                location = Location(latitude=geodfwd[1], longitude=geodfwd[0])
-                time = time + time_to_sail
-
+            adcp_times, adcp_lons, adcp_lats = self._get_underway_measurements(
+                self._expedition.instruments_config.adcp_config,
+                azimuth1,
+                distance_to_move,
+                time_to_reach,
+            )
+
+            for time, lon, lat in zip(adcp_times, adcp_lons, adcp_lats, strict=False):
+                location = Location(latitude=lat, longitude=lon)
                 self._measurements_to_simulate.adcps.append(
                     Spacetime(location=location, time=time)
                 )
 
-                self._next_adcp_time = (
-                    self._next_adcp_time
-                    + self._expedition.instruments_config.adcp_config.period
-                )
-
         # note all ship underwater ST measurements
         if self._expedition.instruments_config.ship_underwater_st_config is not None:
-            location = self._location
-            time = self._time
-            while self._next_ship_underwater_st_time <= end_time:
-                time_to_sail = self._next_ship_underwater_st_time - time
-                distance_to_move = (
-                    ship_speed_meter_per_second * time_to_sail.total_seconds()
-                )
-                geodfwd: tuple[float, float, float] = self._projection.fwd(
-                    lons=location.lon,
-                    lats=location.lat,
-                    az=azimuth1,
-                    dist=distance_to_move,
-                )
-                location = Location(latitude=geodfwd[1], longitude=geodfwd[0])
-                time = time + time_to_sail
-
+            st_times, st_lons, st_lats = self._get_underway_measurements(
+                self._expedition.instruments_config.ship_underwater_st_config,
+                azimuth1,
+                distance_to_move,
+                time_to_reach,
+            )
+
+            for time, lon, lat in zip(st_times, st_lons, st_lats, strict=False):
+                location = Location(latitude=lat, longitude=lon)
                 self._measurements_to_simulate.ship_underwater_sts.append(
                     Spacetime(location=location, time=time)
                 )
 
-                self._next_ship_underwater_st_time = (
-                    self._next_ship_underwater_st_time
-                    + self._expedition.instruments_config.ship_underwater_st_config.period
-                )
-
         self._time = end_time
         self._location = location
 
+    def _get_underway_measurements(
+        self,
+        underway_instrument_config,
+        azimuth: float,
+        distance_to_move: float,
+        time_to_reach: timedelta,
+    ):
+        """Get the times and locations of measurements between current location/waypoint and the next waypoint, for underway instruments."""
+        period = underway_instrument_config.period
+        npts = (time_to_reach.total_seconds() / period.total_seconds()) + 1
+        times = [self._time + i * period for i in range(1, int(npts) + 1)]
+
+        geodfwd = self._projection.fwd_intermediate(
+            lon1=self._location.lon,
+            lat1=self._location.lat,
+            azi1=azimuth,
+            npts=npts,
+            del_s=distance_to_move / npts,
+            return_back_azimuth=False,
+        )
+
+        return times, geodfwd.lons, geodfwd.lats
+
     def _progress_time_stationary(self, time_passed: timedelta) -> None:
+        """Make ship stay at waypoint whilst instruments are deployed, also set the underway instrument measurements that are taken during this time whilst stationary."""
         end_time = self._time + time_passed
 
-        # note all ADCP measurements
+        # note all ADCP measurements (stationary at wp)
         if self._expedition.instruments_config.adcp_config is not None:
-            while self._next_adcp_time <= end_time:
+            adcp_times = self._get_underway_stationary_times(
+                self._expedition.instruments_config.adcp_config, time_passed
+            )
+
+            for time in adcp_times:
                 self._measurements_to_simulate.adcps.append(
-                    Spacetime(self._location, self._next_adcp_time)
-                )
-                self._next_adcp_time = (
-                    self._next_adcp_time
-                    + self._expedition.instruments_config.adcp_config.period
+                    Spacetime(location=self._location, time=time)
                 )
 
-        # note all ship underwater ST measurements
+        # note all underwater ST measurements (stationary at wp)
         if self._expedition.instruments_config.ship_underwater_st_config is not None:
-            while self._next_ship_underwater_st_time <= end_time:
+            st_times = self._get_underway_stationary_times(
+                self._expedition.instruments_config.ship_underwater_st_config,
+                time_passed,
+            )
+            for time in st_times:
                 self._measurements_to_simulate.ship_underwater_sts.append(
-                    Spacetime(self._location, self._next_ship_underwater_st_time)
-                )
-                self._next_ship_underwater_st_time = (
-                    self._next_ship_underwater_st_time
-                    + self._expedition.instruments_config.ship_underwater_st_config.period
+                    Spacetime(location=self._location, time=time)
                 )
 
         self._time = end_time
 
+    def _get_underway_stationary_times(
+        self, underway_instrument_config, time_passed: timedelta
+    ):
+        npts = (
+            time_passed.total_seconds()
+            / underway_instrument_config.period.total_seconds()
+        ) + 1
+        return [
+            self._time + i * underway_instrument_config.period
+            for i in range(1, int(npts) + 1)
+        ]
+
     def _make_measurements(self, waypoint: Waypoint) -> timedelta:
         # if there are no instruments, there is no time cost
         if waypoint.instrument is None:
@@ -268,6 +284,12 @@ def _make_measurements(self, waypoint: Waypoint) -> timedelta:
                         drift_days=self._expedition.instruments_config.argo_float_config.drift_days,
                     )
                 )
+                # TODO: would be good to avoid having to twice make sure that stationkeeping time is factored in; i.e. in schedule validity checks and here (and for CTDs and Drifters)
+                # TODO: makes it easy to forget to update both...
+                # TODO: this is likely to fall under refactoring simulate_schedule.py (i.e. #236)
+                time_costs.append(
+                    self._expedition.instruments_config.argo_float_config.stationkeeping_time
+                )
 
             elif instrument is InstrumentType.CTD:
                 self._measurements_to_simulate.ctds.append(
@@ -302,6 +324,10 @@ def _make_measurements(self, waypoint: Waypoint) -> timedelta:
                         lifetime=self._expedition.instruments_config.drifter_config.lifetime,
                     )
                 )
+                time_costs.append(
+                    self._expedition.instruments_config.drifter_config.stationkeeping_time
+                )
+
             elif instrument is InstrumentType.XBT:
                 self._measurements_to_simulate.xbts.append(
                     XBT(
@@ -315,5 +341,7 @@ def _make_measurements(self, waypoint: Waypoint) -> timedelta:
             else:
                 raise NotImplementedError("Instrument type not supported.")
 
-        # measurements are done in parallel, so return time of longest one
+        # measurements are done simultaneously onboard, so return time of longest one
+        # TODO: docs suggest that add individual instrument stationkeeping times are cumulative, which is at odds with measurements being done simultaneously onboard here
+        # TODO: update one or the other?
         return max(time_costs)

tests/expedition/test_simulate_schedule.py

@@ -1,5 +1,6 @@
 from datetime import datetime, timedelta
 
+import numpy as np
 import pyproj
 
 from virtualship.expedition.simulate_schedule import (
@@ -65,3 +66,62 @@ def test_time_in_minutes_in_ship_schedule() -> None:
         minutes=20
     )
     assert instruments_config.ship_underwater_st_config.period == timedelta(minutes=5)
+
+
+def test_ship_path_inside_domain() -> None:
+    """Test that the ship path (here represented by underway ADCP measurement sites) is inside the domain defined by the waypoints (which determines the fieldset bounds)."""
+    base_time = datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S")
+
+    projection = pyproj.Geod(ellps="WGS84")
+    expedition = Expedition.from_yaml("expedition_dir/expedition.yaml")
+    expedition.ship_config.ship_speed_knots = 10.0
+
+    wp1 = Location(-63.0, -57.7)  # most southern
+    wp2 = Location(-55.4, -66.2)  # most northern
+    wp3 = Location(-61.8, -73.2)  # most western
+    wp4 = Location(-57.3, -51.8)  # most eastern
+
+    # waypoints with enough distance where curvature is clear
+    expedition.schedule = Schedule(
+        waypoints=[
+            Waypoint(location=wp1, time=base_time),
+            Waypoint(location=wp2, time=base_time + timedelta(days=5)),
+            Waypoint(location=wp3, time=base_time + timedelta(days=10)),
+            Waypoint(location=wp4, time=base_time + timedelta(days=15)),
+        ]
+    )
+
+    # get waypoint domain bounds
+    wp_max_lat, wp_min_lat, wp_max_lon, wp_min_lon = (
+        max(wp.location.lat for wp in expedition.schedule.waypoints),
+        min(wp.location.lat for wp in expedition.schedule.waypoints),
+        max(wp.location.lon for wp in expedition.schedule.waypoints),
+        min(wp.location.lon for wp in expedition.schedule.waypoints),
+    )
+
+    result = simulate_schedule(projection, expedition)
+    assert isinstance(result, ScheduleOk)
+
+    # adcp measurements path
+    adcp_measurements = result.measurements_to_simulate.adcps
+    adcp_lats = [m.location.lat for m in adcp_measurements]
+    adcp_lons = [m.location.lon for m in adcp_measurements]
+
+    adcp_max_lat, adcp_min_lat, adcp_max_lon, adcp_min_lon = (
+        max(adcp_lats),
+        min(adcp_lats),
+        max(adcp_lons),
+        min(adcp_lons),
+    )
+
+    # check adcp route is within wp bounds
+    assert adcp_max_lat <= wp_max_lat
+    assert adcp_min_lat >= wp_min_lat
+    assert adcp_max_lon <= wp_max_lon
+    assert adcp_min_lon >= wp_min_lon
+
+    # the adcp route extremes should also approximately match waypoints defined in this test
+    assert np.isclose(adcp_max_lat, wp2.lat, atol=0.1)
+    assert np.isclose(adcp_min_lat, wp1.lat, atol=0.1)
+    assert np.isclose(adcp_max_lon, wp4.lon, atol=0.1)
+    assert np.isclose(adcp_min_lon, wp3.lon, atol=0.1)