verify waypoints are on water

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)

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)