diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 5b76fb52..d96efc91 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,9 +140,11 @@ 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, @@ -150,90 +152,104 @@ 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 - + 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) diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py index ff5a3597..f90eecc2 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 ( @@ -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)