diff --git a/tests/test_sailship.py b/tests/test_sailship.py index f7e42027..f3691d49 100644 --- a/tests/test_sailship.py +++ b/tests/test_sailship.py @@ -4,10 +4,12 @@ from datetime import timedelta import numpy as np +import pyproj +import pytest from parcels import Field, FieldSet from virtual_ship import InstrumentType, Location, Waypoint -from virtual_ship.sailship import sailship +from virtual_ship.sailship import PlanningError, _verify_waypoints, sailship from virtual_ship.virtual_ship_config import ( ADCPConfig, ArgoFloatConfig, @@ -19,10 +21,10 @@ def _make_ctd_fieldset(base_time: datetime) -> FieldSet: - u = np.zeros((2, 2, 2, 2)) - v = np.zeros((2, 2, 2, 2)) - t = np.zeros((2, 2, 2, 2)) - s = np.zeros((2, 2, 2, 2)) + u = np.full((2, 2, 2, 2), 1.0) + v = np.full((2, 2, 2, 2), 1.0) + t = np.full((2, 2, 2, 2), 1.0) + s = np.full((2, 2, 2, 2), 1.0) fieldset = FieldSet.from_data( {"V": v, "U": u, "T": t, "S": s}, @@ -64,12 +66,12 @@ def test_sailship() -> None: base_time = datetime.datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S") adcp_fieldset = FieldSet.from_data( - {"U": 0, "V": 0}, + {"U": 1, "V": 1}, {"lon": 0, "lat": 0}, ) ship_underwater_st_fieldset = FieldSet.from_data( - {"U": 0, "V": 0, "salinity": 0, "temperature": 0}, + {"U": 1, "V": 1, "salinity": 0, "temperature": 0}, {"lon": 0, "lat": 0}, ) @@ -78,7 +80,7 @@ def test_sailship() -> None: drifter_fieldset = _make_drifter_fieldset(base_time) argo_float_fieldset = FieldSet.from_data( - {"U": 0, "V": 0, "T": 0, "S": 0}, + {"U": 1, "V": 1, "T": 0, "S": 0}, { "lon": 0, "lat": 0, @@ -155,3 +157,78 @@ def test_sailship() -> None: ) sailship(config) + + +def test_verify_waypoints() -> None: + # arbitrary cruise start time + BASE_TIME = datetime.datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S") + PROJECTION = pyproj.Geod(ellps="WGS84") + + # the sets of waypoints to test + WAYPOINTS = [ + [], # require at least one waypoint + [Waypoint(Location(0.0, 0.0))], # first waypoint must have time + [ + Waypoint(Location(0.0, 0.0), BASE_TIME + datetime.timedelta(days=1)), + Waypoint(Location(0.0, 0.0), BASE_TIME), + ], # waypoint times must be in ascending order + [ + Waypoint(Location(0.0, 0.0), BASE_TIME), + ], # 0 uv points are on land + [ + Waypoint(Location(0.1, 0.1), BASE_TIME), + Waypoint(Location(1.0, 1.0), BASE_TIME + datetime.timedelta(seconds=1)), + ], # waypoints must be reachable in time + [ + Waypoint(Location(0.1, 0.1), BASE_TIME), + Waypoint(Location(1.0, 1.0), BASE_TIME + datetime.timedelta(days=1)), + ], # a valid schedule + ] + + # the expected errors for the schedules, or None if expected to be valid + EXPECT_MATCH = [ + "^At least one waypoint must be provided.$", + "^First waypoint must have a specified time.$", + "^Each waypoint should be timed after all previous waypoints$", + "^The following waypoints are on land: .*$", + "^Waypoint planning is not valid: would arrive too late at a waypoint number .*$", + None, + ] + + # create a fieldset matching the test waypoints + u = np.full((1, 1, 2, 2), 1.0) + v = np.full((1, 1, 2, 2), 1.0) + u[0, 0, 0, 0] = 0.0 + v[0, 0, 0, 0] = 0.0 + + fieldset = FieldSet.from_data( + {"V": v, "U": u}, + { + "time": [np.datetime64(BASE_TIME)], + "depth": [0], + "lat": [0, 1], + "lon": [0, 1], + }, + ) + + # dummy configs + ctd_config = CTDConfig(None, fieldset, None, None) + drifter_config = DrifterConfig(None, None, None) + argo_float_config = ArgoFloatConfig(None, None, None, None, None, None, None) + + # test each set of waypoints and verify the raised errors (or none if valid) + for waypoints, expect_match in zip(WAYPOINTS, EXPECT_MATCH, strict=True): + config = VirtualShipConfig( + ship_speed=5.14, + waypoints=waypoints, + argo_float_config=argo_float_config, + adcp_config=None, + ship_underwater_st_config=None, + ctd_config=ctd_config, + drifter_config=drifter_config, + ) + if expect_match is not None: + with pytest.raises(PlanningError, match=expect_match): + _verify_waypoints(PROJECTION, config) + else: + _verify_waypoints(PROJECTION, config) diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py index ea92fa24..4569135b 100644 --- a/virtual_ship/sailship.py +++ b/virtual_ship/sailship.py @@ -9,6 +9,7 @@ from typing import Generator import pyproj +from parcels import FieldSet from sortedcontainers import SortedList from .costs import costs @@ -36,7 +37,7 @@ def sailship(config: VirtualShipConfig): # projection used to sail between waypoints projection = pyproj.Geod(ellps="WGS84") - _verify_waypoints(config.waypoints, config.ship_speed, projection=projection) + _verify_waypoints(projection=projection, config=config) # simulate the sailing and aggregate what measurements should be simulated schedule_results = _simulate_schedule( @@ -386,22 +387,26 @@ def _argo_float_task( def _verify_waypoints( - waypoints: list[Waypoint], ship_speed: float, projection: pyproj.Geod + projection: pyproj.Geod, + config: VirtualShipConfig, ) -> None: """ Verify waypoints are ordered by time, first waypoint has a start time, and that schedule is feasible in terms of time if no unexpected events happen. - :param waypoints: The waypoints to check. - :param ship_speed: Speed of the ship. :param projection: projection used to sail between waypoints. + :param config: The cruise configuration. :raises PlanningError: If waypoints are not feasible or incorrect. + :raises ValueError: If there are no fieldsets in the config, which are needed to verify all waypoints are on water. """ + if len(config.waypoints) == 0: + raise PlanningError("At least one waypoint must be provided.") + # check first waypoint has a time - if waypoints[0].time is None: + if config.waypoints[0].time is None: raise PlanningError("First waypoint must have a specified time.") # check waypoint times are in ascending order - timed_waypoints = [wp for wp in waypoints if wp.time is not None] + timed_waypoints = [wp for wp in config.waypoints if wp.time is not None] if not all( [ next.time >= cur.time @@ -412,9 +417,47 @@ def _verify_waypoints( "Each waypoint should be timed after all previous waypoints" ) + # check if all waypoints are in water + # this is done by picking an arbitrary provided fieldset and checking if UV is not zero + + # get all available fieldsets + available_fieldsets = [ + fs + for fs in [ + config.adcp_config.fieldset if config.adcp_config is not None else None, + config.argo_float_config.fieldset, + config.ctd_config.fieldset, + config.drifter_config.fieldset, + ( + config.ship_underwater_st_config.fieldset + if config.ship_underwater_st_config is not None + else None + ), + ] + if fs is not None + ] + # check if there are any fieldsets, else its an error + if len(available_fieldsets) == 0: + raise ValueError( + "No fieldsets provided to check if waypoints are on land. Assuming no provided fieldsets is an error." + ) + # pick any + fieldset = available_fieldsets[0] + # get waypoints with 0 UV + land_waypoints = [ + (wp_i, wp) + for wp_i, wp in enumerate(config.waypoints) + if _is_on_land_zero_uv(fieldset, wp) + ] + # raise an error if there are any + if len(land_waypoints) > 0: + raise PlanningError( + f"The following waypoints are on land: {['#' + str(wp_i) + ' ' + str(wp) for (wp_i, wp) in land_waypoints]}" + ) + # check that ship will arrive on time at each waypoint (in case no unexpected event happen) - time = waypoints[0].time - for wp_i, (wp, wp_next) in enumerate(zip(waypoints, waypoints[1:])): + time = config.waypoints[0].time + for wp_i, (wp, wp_next) in enumerate(zip(config.waypoints, config.waypoints[1:])): if wp.instrument is InstrumentType.CTD: time += timedelta(minutes=20) @@ -423,7 +466,7 @@ def _verify_waypoints( ) distance = geodinv[2] - time_to_reach = timedelta(seconds=distance / ship_speed) + time_to_reach = timedelta(seconds=distance / config.ship_speed) arrival_time = time + time_to_reach if wp_next.time is None: @@ -434,3 +477,16 @@ def _verify_waypoints( ) else: time = wp_next.time + + +def _is_on_land_zero_uv(fieldset: FieldSet, waypoint: Waypoint) -> bool: + """ + Check if waypoint is on land by assuming zero velocity means land. + + :param fieldset: The fieldset to sample the velocity from. + :param waypoint: The waypoint to check. + :returns: If the waypoint is on land. + """ + return fieldset.UV.eval( + 0, 0, waypoint.location.lat, waypoint.location.lon, applyConversion=False + ) == (0.0, 0.0)