From b9d53c92bd1921dec7b93a6724f779a67b51f8e8 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 2 Mar 2026 20:24:44 +0000 Subject: [PATCH 1/8] fix ship trajectory by using geod.fwd_intermediate() to get all points across the shortest path on globe --- .../expedition/simulate_schedule.py | 128 ++++++++++-------- 1 file changed, 70 insertions(+), 58 deletions(-) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 5b76fb52..80d3f23e 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/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,6 +140,7 @@ 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: @@ -150,90 +151,101 @@ def _progress_time_traveling_towards(self, location: Location) -> None: 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 - - self._measurements_to_simulate.ship_underwater_sts.append( + 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.adcps.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, + distance_to_move, + time_to_reach: timedelta, + ): + """Get the times and locations of measurements between a waypoint (or the start) 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: 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): + 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: From 9818ae713b69b24667980dbe8227095e808f0672 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 2 Mar 2026 20:25:15 +0000 Subject: [PATCH 2/8] add missing time costs associated with argo float and drifter deployments --- src/virtualship/expedition/simulate_schedule.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 80d3f23e..42c17bed 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/src/virtualship/expedition/simulate_schedule.py @@ -280,6 +280,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( @@ -314,6 +320,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( @@ -327,5 +337,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) From c2eccf215c86e35fa9b72695f6a10002c1d35d27 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 07:02:11 +0000 Subject: [PATCH 3/8] fix, underwater st instead of adcp --- src/virtualship/expedition/simulate_schedule.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 42c17bed..7d043bd3 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/src/virtualship/expedition/simulate_schedule.py @@ -179,7 +179,7 @@ def _progress_time_traveling_towards(self, location: Location) -> None: for time, lon, lat in zip(st_times, st_lons, st_lats, strict=False): location = Location(latitude=lat, longitude=lon) - self._measurements_to_simulate.adcps.append( + self._measurements_to_simulate.ship_underwater_sts.append( Spacetime(location=location, time=time) ) From 17d770d20fa432088e63c2fdaeaa51f7a36a3adc Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 07:03:53 +0000 Subject: [PATCH 4/8] typechecking --- src/virtualship/expedition/simulate_schedule.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 7d043bd3..a04f8bae 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/src/virtualship/expedition/simulate_schedule.py @@ -189,8 +189,8 @@ def _progress_time_traveling_towards(self, location: Location) -> None: def _get_underway_measurements( self, underway_instrument_config, - azimuth, - distance_to_move, + azimuth: float, + distance_to_move: float, time_to_reach: timedelta, ): """Get the times and locations of measurements between a waypoint (or the start) and the next waypoint, for underway instruments.""" @@ -236,7 +236,9 @@ def _progress_time_stationary(self, time_passed: timedelta) -> None: self._time = end_time - def _get_underway_stationary_times(self, underway_instrument_config, time_passed): + def _get_underway_stationary_times( + self, underway_instrument_config, time_passed: timedelta + ): npts = ( time_passed.total_seconds() / underway_instrument_config.period.total_seconds() From 18598528ce15e4c3dee9a4c1db15121481a6be99 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 08:36:37 +0100 Subject: [PATCH 5/8] update/add docstrings --- src/virtualship/expedition/simulate_schedule.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index a04f8bae..d96efc91 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/src/virtualship/expedition/simulate_schedule.py @@ -144,6 +144,7 @@ def simulate(self) -> ScheduleOk | ScheduleProblem: 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, @@ -193,7 +194,7 @@ def _get_underway_measurements( distance_to_move: float, time_to_reach: timedelta, ): - """Get the times and locations of measurements between a waypoint (or the start) and the next waypoint, for underway instruments.""" + """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)] @@ -210,6 +211,7 @@ def _get_underway_measurements( 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 (stationary at wp) From 61bff39a4138c04f7291238b2a8c1c3d075702a5 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 10:52:22 +0100 Subject: [PATCH 6/8] add new test for ship track in waypoint/fieldset domain --- tests/expedition/test_simulate_schedule.py | 53 ++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py index ff5a3597..aa71e28b 100644 --- a/tests/expedition/test_simulate_schedule.py +++ b/tests/expedition/test_simulate_schedule.py @@ -65,3 +65,56 @@ 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 (and therefore where underway measurements will be taken) 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 + + # waypoints with enough distance where curvature is clear + expedition.schedule = Schedule( + waypoints=[ + Waypoint(location=Location(-63.0, -57.7), time=base_time), + Waypoint( + location=Location(-55.4, -66.2), time=base_time + timedelta(days=5) + ), + Waypoint( + location=Location(-61.8, -73.2), time=base_time + timedelta(days=10) + ), + Waypoint( + location=Location(-57.3, -51.8), 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), + ) + + 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 From c522eb168e921010f6e162d0eb84a42160b8c200 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 10:58:44 +0100 Subject: [PATCH 7/8] add more checks to the test --- tests/expedition/test_simulate_schedule.py | 26 +++++++++++++--------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py index aa71e28b..b71bd9f6 100644 --- a/tests/expedition/test_simulate_schedule.py +++ b/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 ( @@ -75,19 +76,18 @@ def test_ship_path_inside_domain() -> None: 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=Location(-63.0, -57.7), time=base_time), - Waypoint( - location=Location(-55.4, -66.2), time=base_time + timedelta(days=5) - ), - Waypoint( - location=Location(-61.8, -73.2), time=base_time + timedelta(days=10) - ), - Waypoint( - location=Location(-57.3, -51.8), time=base_time + timedelta(days=15) - ), + 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)), ] ) @@ -118,3 +118,9 @@ def test_ship_path_inside_domain() -> None: 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) From a86cfa2e1415d6738d723cd50651837eca190708 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 3 Mar 2026 11:00:35 +0100 Subject: [PATCH 8/8] update docstring --- tests/expedition/test_simulate_schedule.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py index b71bd9f6..f90eecc2 100644 --- a/tests/expedition/test_simulate_schedule.py +++ b/tests/expedition/test_simulate_schedule.py @@ -69,7 +69,7 @@ def test_time_in_minutes_in_ship_schedule() -> None: def test_ship_path_inside_domain() -> None: - """Test that the ship path (and therefore where underway measurements will be taken) is inside the domain defined by the waypoints (which determines the fieldset bounds).""" + """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") @@ -114,6 +114,7 @@ def test_ship_path_inside_domain() -> None: 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