diff --git a/pyproject.toml b/pyproject.toml index 033d157a..c9867834 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,6 +66,8 @@ dev = [ "pytest == 8.2.0", "pytest-cov == 5.0.0", "codecov == 2.1.13", + "sortedcontainers == 2.4.0", + "sortedcontainers-stubs == 2.4.2", ] [tool.isort] diff --git a/tests/test_sailship.py b/tests/test_sailship.py index 2f3dd974..f7e42027 100644 --- a/tests/test_sailship.py +++ b/tests/test_sailship.py @@ -1,15 +1,19 @@ """Performs a complete cruise with virtual ship.""" import datetime +from datetime import timedelta import numpy as np from parcels import Field, FieldSet -from virtual_ship import Location +from virtual_ship import InstrumentType, Location, Waypoint from virtual_ship.sailship import sailship -from virtual_ship.virtual_ship_configuration import ( +from virtual_ship.virtual_ship_config import ( ADCPConfig, ArgoFloatConfig, + CTDConfig, + DrifterConfig, + ShipUnderwaterSTConfig, VirtualShipConfig, ) @@ -65,7 +69,7 @@ def test_sailship() -> None: ) ship_underwater_st_fieldset = FieldSet.from_data( - {"U": 0, "V": 0, "S": 0, "T": 0}, + {"U": 0, "V": 0, "salinity": 0, "temperature": 0}, {"lon": 0, "lat": 0}, ) @@ -84,6 +88,7 @@ def test_sailship() -> None: argo_float_config = ArgoFloatConfig( fieldset=argo_float_fieldset, + min_depth=-argo_float_fieldset.U.depth[0], max_depth=-2000, drift_depth=-1000, vertical_speed=-0.10, @@ -91,28 +96,62 @@ def test_sailship() -> None: drift_days=9, ) - adcp_config = ADCPConfig(max_depth=-1000, bin_size_m=24) + adcp_config = ADCPConfig( + max_depth=-1000, + bin_size_m=24, + period=timedelta(minutes=5), + fieldset=adcp_fieldset, + ) - config = VirtualShipConfig( - start_time=datetime.datetime.strptime( - "2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S" + ship_underwater_st_config = ShipUnderwaterSTConfig( + period=timedelta(minutes=5), fieldset=ship_underwater_st_fieldset + ) + + ctd_config = CTDConfig( + stationkeeping_time=timedelta(minutes=20), + fieldset=ctd_fieldset, + min_depth=ctd_fieldset.U.depth[0], + max_depth=ctd_fieldset.U.depth[-1], + ) + + drifter_config = DrifterConfig( + fieldset=drifter_fieldset, + depth=-drifter_fieldset.U.depth[0], + lifetime=timedelta(weeks=4), + ) + + waypoints = [ + Waypoint( + location=Location(latitude=-23.071289, longitude=63.743631), + time=base_time, + ), + Waypoint( + location=Location(latitude=-23.081289, longitude=63.743631), + instrument=InstrumentType.CTD, + ), + Waypoint( + location=Location(latitude=-23.181289, longitude=63.743631), + time=base_time + datetime.timedelta(hours=1), + instrument=InstrumentType.CTD, + ), + Waypoint( + location=Location(latitude=-23.281289, longitude=63.743631), + instrument=InstrumentType.DRIFTER, ), - route_coordinates=[ - Location(latitude=-23.071289, longitude=63.743631), - Location(latitude=-23.081289, longitude=63.743631), - Location(latitude=-23.191289, longitude=63.743631), - ], - adcp_fieldset=adcp_fieldset, - ship_underwater_st_fieldset=ship_underwater_st_fieldset, - ctd_fieldset=ctd_fieldset, - drifter_fieldset=drifter_fieldset, - argo_float_deploy_locations=[ - Location(latitude=-23.081289, longitude=63.743631) - ], - drifter_deploy_locations=[Location(latitude=-23.081289, longitude=63.743631)], - ctd_deploy_locations=[Location(latitude=-23.081289, longitude=63.743631)], + Waypoint( + location=Location(latitude=-23.381289, longitude=63.743631), + instrument=InstrumentType.ARGO_FLOAT, + ), + ] + + config = VirtualShipConfig( + ship_speed=5.14, + waypoints=waypoints, argo_float_config=argo_float_config, adcp_config=adcp_config, + ship_underwater_st_config=ship_underwater_st_config, + ctd_config=ctd_config, + drifter_config=drifter_config, ) sailship(config) diff --git a/virtual_ship/__init__.py b/virtual_ship/__init__.py index dbf175e3..31a64875 100644 --- a/virtual_ship/__init__.py +++ b/virtual_ship/__init__.py @@ -1,7 +1,18 @@ """Code for the Virtual Ship Classroom, where Marine Scientists can combine Copernicus Marine Data with an OceanParcels ship to go on a virtual expedition.""" from . import instruments, sailship +from .instrument_type import InstrumentType from .location import Location +from .planning_error import PlanningError from .spacetime import Spacetime +from .waypoint import Waypoint -__all__ = ["Location", "Spacetime", "instruments", "sailship"] +__all__ = [ + "InstrumentType", + "Location", + "PlanningError", + "Spacetime", + "Waypoint", + "instruments", + "sailship", +] diff --git a/virtual_ship/costs.py b/virtual_ship/costs.py index 459dba67..4339eca4 100644 --- a/virtual_ship/costs.py +++ b/virtual_ship/costs.py @@ -2,7 +2,8 @@ from datetime import timedelta -from .virtual_ship_configuration import VirtualShipConfig +from .instrument_type import InstrumentType +from .virtual_ship_config import VirtualShipConfig def costs(config: VirtualShipConfig, total_time: timedelta): @@ -18,8 +19,22 @@ def costs(config: VirtualShipConfig, total_time: timedelta): argo_deploy_cost = 15000 ship_cost = ship_cost_per_day / 24 * total_time.total_seconds() // 3600 - argo_cost = len(config.argo_float_deploy_locations) * argo_deploy_cost - drifter_cost = len(config.drifter_deploy_locations) * drifter_deploy_cost + num_argos = len( + [ + waypoint + for waypoint in config.waypoints + if waypoint.instrument is InstrumentType.ARGO_FLOAT + ] + ) + argo_cost = num_argos * argo_deploy_cost + num_drifters = len( + [ + waypoint + for waypoint in config.waypoints + if waypoint.instrument is InstrumentType.DRIFTER + ] + ) + drifter_cost = num_drifters * drifter_deploy_cost cost = ship_cost + argo_cost + drifter_cost return cost diff --git a/virtual_ship/instrument_type.py b/virtual_ship/instrument_type.py new file mode 100644 index 00000000..9a19b814 --- /dev/null +++ b/virtual_ship/instrument_type.py @@ -0,0 +1,11 @@ +"""InstrumentType Enum.""" + +from enum import Enum, auto + + +class InstrumentType(Enum): + """Types of instruments.""" + + CTD = auto() + DRIFTER = auto() + ARGO_FLOAT = auto() diff --git a/virtual_ship/instruments/adcp.py b/virtual_ship/instruments/adcp.py index 096facf1..bca4ab94 100644 --- a/virtual_ship/instruments/adcp.py +++ b/virtual_ship/instruments/adcp.py @@ -31,7 +31,7 @@ def simulate_adcp( sample_points: list[Spacetime], ) -> None: """ - Use parcels to simulate an ADCP in a fieldset. + Use Parcels to simulate an ADCP in a fieldset. :param fieldset: The fieldset to simulate the ADCP in. :param out_path: The path to write the results to. diff --git a/virtual_ship/instruments/argo_float.py b/virtual_ship/instruments/argo_float.py index 031c8309..ad4a4871 100644 --- a/virtual_ship/instruments/argo_float.py +++ b/virtual_ship/instruments/argo_float.py @@ -123,7 +123,7 @@ def simulate_argo_floats( endtime: datetime | None, ) -> None: """ - Use parcels to simulate a set of Argo floats in a fieldset. + Use Parcels to simulate a set of Argo floats in a fieldset. :param fieldset: The fieldset to simulate the Argo floats in. :param out_path: The path to write the results to. @@ -133,6 +133,13 @@ def simulate_argo_floats( """ DT = 10.0 # dt of Argo float simulation integrator + if len(argo_floats) == 0: + print( + "No Argo floats provided. Parcels currently crashes when providing an empty particle set, so no argo floats simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + # define parcel particles argo_float_particleset = ParticleSet( fieldset=fieldset, diff --git a/virtual_ship/instruments/ctd.py b/virtual_ship/instruments/ctd.py index cc8734c9..97057f46 100644 --- a/virtual_ship/instruments/ctd.py +++ b/virtual_ship/instruments/ctd.py @@ -60,7 +60,7 @@ def simulate_ctd( outputdt: timedelta, ) -> None: """ - Use parcels to simulate a set of CTDs in a fieldset. + Use Parcels to simulate a set of CTDs in a fieldset. :param fieldset: The fieldset to simulate the CTDs in. :param out_path: The path to write the results to. @@ -71,6 +71,13 @@ def simulate_ctd( WINCH_SPEED = 1.0 # sink and rise speed in m/s DT = 10.0 # dt of CTD simulation integrator + if len(ctds) == 0: + print( + "No CTDs provided. Parcels currently crashes when providing an empty particle set, so no CTD simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + fieldset_starttime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[0]) fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) diff --git a/virtual_ship/instruments/drifter.py b/virtual_ship/instruments/drifter.py index 3a5be7ab..7854d4cb 100644 --- a/virtual_ship/instruments/drifter.py +++ b/virtual_ship/instruments/drifter.py @@ -49,7 +49,7 @@ def simulate_drifters( endtime: datetime | None, ) -> None: """ - Use parcels to simulate a set of drifters in a fieldset. + Use Parcels to simulate a set of drifters in a fieldset. :param fieldset: The fieldset to simulate the Drifters in. :param out_path: The path to write the results to. @@ -58,6 +58,13 @@ def simulate_drifters( :param dt: Dt for integration. :param endtime: Stop at this time, or if None, continue until the end of the fieldset or until all drifters ended. If this is earlier than the last drifter ended or later than the end of the fieldset, a warning will be printed. """ + if len(drifters) == 0: + print( + "No drifters provided. Parcels currently crashes when providing an empty particle set, so no drifter simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + # define parcel particles drifter_particleset = ParticleSet( fieldset=fieldset, diff --git a/virtual_ship/instruments/ship_underwater_st.py b/virtual_ship/instruments/ship_underwater_st.py index b2226381..3b4ac59e 100644 --- a/virtual_ship/instruments/ship_underwater_st.py +++ b/virtual_ship/instruments/ship_underwater_st.py @@ -37,7 +37,7 @@ def simulate_ship_underwater_st( sample_points: list[Spacetime], ) -> None: """ - Use parcels to simulate underway data, measuring salinity and temperature at the given depth along the ship track in a fieldset. + Use Parcels to simulate underway data, measuring salinity and temperature at the given depth along the ship track in a fieldset. :param fieldset: The fieldset to simulate the sampling in. :param out_path: The path to write the results to. diff --git a/virtual_ship/planning_error.py b/virtual_ship/planning_error.py new file mode 100644 index 00000000..5d599dd2 --- /dev/null +++ b/virtual_ship/planning_error.py @@ -0,0 +1,7 @@ +"""PlanningError Exception.""" + + +class PlanningError(RuntimeError): + """An error when checking the schedule or during sailing.""" + + pass diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py index 9a02c828..ea92fa24 100644 --- a/virtual_ship/sailship.py +++ b/virtual_ship/sailship.py @@ -1,219 +1,86 @@ """sailship function.""" +from __future__ import annotations + import os -from datetime import timedelta +from contextlib import contextmanager +from dataclasses import dataclass, field +from datetime import datetime, timedelta +from typing import Generator -import numpy as np import pyproj +from sortedcontainers import SortedList from .costs import costs +from .instrument_type import InstrumentType from .instruments.adcp import simulate_adcp from .instruments.argo_float import ArgoFloat, simulate_argo_floats from .instruments.ctd import CTD, simulate_ctd from .instruments.drifter import Drifter, simulate_drifters from .instruments.ship_underwater_st import simulate_ship_underwater_st from .location import Location -from .postprocess import postprocess +from .planning_error import PlanningError from .spacetime import Spacetime -from .virtual_ship_configuration import VirtualShipConfig +from .virtual_ship_config import VirtualShipConfig +from .waypoint import Waypoint def sailship(config: VirtualShipConfig): """ - Use parcels to simulate the ship, take ctd_instruments and measure ADCP and underwaydata. + Use Parcels to simulate a virtual ship expedition. - :param config: The cruise configuration. + :param config: The expedition configuration. """ config.verify() - # combine identical instrument deploy location - argo_locations = set(config.argo_float_deploy_locations) - if len(argo_locations) != len(config.argo_float_deploy_locations): - print( - "WARN: Some argo float deployment locations are identical and have been combined." - ) + # projection used to sail between waypoints + projection = pyproj.Geod(ellps="WGS84") - drifter_locations = set(config.drifter_deploy_locations) - if len(drifter_locations) != len(config.drifter_deploy_locations): - print( - "WARN: Some drifter deployment locations are identical and have been combined." - ) + _verify_waypoints(config.waypoints, config.ship_speed, projection=projection) - ctd_locations = set(config.ctd_deploy_locations) - if len(drifter_locations) != len(config.ctd_deploy_locations): - print("WARN: Some CTD locations are identical and have been combined.") - - # get discrete points along the ships route were sampling and deployments will be performed - route_dt = timedelta(minutes=5) - route_points = shiproute(config=config, dt=route_dt) - - # adcp objects to be used in simulation - adcps: list[Spacetime] = [] - - # ship st objects to be used in simulation - ship_underwater_sts: list[Spacetime] = [] - - # argo float deployment locations that have been visited - argo_locations_visited: set[Location] = set() - # argo float objects to be used in simulation - argo_floats: list[ArgoFloat] = [] - - # drifter deployment locations that have been visited - drifter_locations_visited: set[Location] = set() - # drifter objects to be used in simulation - drifters: list[Drifter] = [] - - # ctd cast locations that have been visited - ctd_locations_visited: set[Location] = set() - # ctd cast objects to be used in simulation - ctds: list[CTD] = [] - - # iterate over each discrete route point, find deployment and measurement locations and times, and measure how much time this took - # TODO between drifters, argo floats, ctd there is quite some repetition - print("Traversing ship route.") - time_past = timedelta() - for i, route_point in enumerate(route_points): - if i % 96 == 0: - print(f"Gathered data {time_past} hours since start.") - - # find drifter deployments to be done at this location - drifters_here = set( - [ - drifter - for drifter in drifter_locations - drifter_locations_visited - if all( - np.isclose( - [drifter.lat, drifter.lon], [route_point.lat, route_point.lon] - ) - ) - ] - ) - if len(drifters_here) > 1: - print( - "WARN: Multiple drifter deployments match the current location. Only a single deployment will be performed." - ) - drifters.append( - Drifter( - spacetime=Spacetime( - location=route_point, time=time_past.total_seconds() - ), - depth=-config.drifter_fieldset.U.depth[0], - lifetime=timedelta(weeks=4), - ) - ) - drifter_locations_visited = drifter_locations_visited.union(drifters_here) - - # find argo float deployments to be done at this location - argos_here = set( - [ - argo - for argo in argo_locations - argo_locations_visited - if all( - np.isclose([argo.lat, argo.lon], [route_point.lat, route_point.lon]) - ) - ] - ) - if len(argos_here) > 1: - print( - "WARN: Multiple argo float deployments match the current location. Only a single deployment will be performed." - ) - argo_floats.append( - ArgoFloat( - spacetime=Spacetime( - location=route_point, time=time_past.total_seconds() - ), - min_depth=-config.argo_float_config.fieldset.U.depth[0], - max_depth=config.argo_float_config.max_depth, - drift_depth=config.argo_float_config.drift_depth, - vertical_speed=config.argo_float_config.vertical_speed, - cycle_days=config.argo_float_config.cycle_days, - drift_days=config.argo_float_config.drift_days, - ) - ) - argo_locations_visited = argo_locations_visited.union(argos_here) - - # find CTD casts to be done at this location - ctds_here = set( - [ - ctd - for ctd in ctd_locations - ctd_locations_visited - if all( - np.isclose([ctd.lat, ctd.lon], [route_point.lat, route_point.lon]) - ) - ] - ) - if len(ctds_here) > 1: - print( - "WARN: Multiple CTD casts match the current location. Only a single cast will be performed." - ) + # simulate the sailing and aggregate what measurements should be simulated + schedule_results = _simulate_schedule( + waypoints=config.waypoints, + projection=projection, + config=config, + ) - ctds.append( - CTD( - spacetime=Spacetime( - location=route_point, - time=config.start_time + time_past, - ), - min_depth=config.ctd_fieldset.U.depth[0], - max_depth=config.ctd_fieldset.U.depth[-1], - ) - ) - ctd_locations_visited = ctd_locations_visited.union(ctds_here) - # add 20 minutes to sailing time for deployment - if len(ctds_here) != 0: - time_past += timedelta(minutes=20) - - # add time it takes to move to the next route point - time_past += route_dt - # remove the last one, because no sailing to the next point was needed - time_past -= route_dt - - # check if all drifter, argo float, ctd locations were visited - if len(drifter_locations_visited) != len(drifter_locations): - print( - "WARN: some drifter deployments were not planned along the route and have not been performed." - ) + # simulate the measurements - if len(argo_locations_visited) != len(argo_locations): - print( - "WARN: some argo float deployments were not planned along the route and have not been performed." + if config.ship_underwater_st_config is not None: + print("Simulating onboard salinity and temperature measurements.") + simulate_ship_underwater_st( + fieldset=config.ship_underwater_st_config.fieldset, + out_path=os.path.join("results", "ship_underwater_st.zarr"), + depth=-2, + sample_points=schedule_results.measurements_to_simulate.ship_underwater_sts, ) - if len(ctd_locations_visited) != len(ctd_locations): - print( - "WARN: some CTD casts were not planned along the route and have not been performed." + if config.adcp_config is not None: + print("Simulating onboard ADCP.") + simulate_adcp( + fieldset=config.adcp_config.fieldset, + out_path=os.path.join("results", "adcp.zarr"), + max_depth=config.adcp_config.max_depth, + min_depth=-5, + num_bins=(-5 - config.adcp_config.max_depth) + // config.adcp_config.bin_size_m, + sample_points=schedule_results.measurements_to_simulate.adcps, ) - print("Simulating onboard salinity and temperature measurements.") - simulate_ship_underwater_st( - fieldset=config.ship_underwater_st_fieldset, - out_path=os.path.join("results", "ship_underwater_st.zarr"), - depth=-2, - sample_points=ship_underwater_sts, - ) - - print("Simulating onboard ADCP.") - simulate_adcp( - fieldset=config.adcp_fieldset, - out_path=os.path.join("results", "adcp.zarr"), - max_depth=config.adcp_config.max_depth, - min_depth=-5, - num_bins=(-5 - config.adcp_config.max_depth) // config.adcp_config.bin_size_m, - sample_points=adcps, - ) - print("Simulating CTD casts.") simulate_ctd( out_path=os.path.join("results", "ctd.zarr"), - fieldset=config.ctd_fieldset, - ctds=ctds, + fieldset=config.ctd_config.fieldset, + ctds=schedule_results.measurements_to_simulate.ctds, outputdt=timedelta(seconds=10), ) print("Simulating drifters") simulate_drifters( out_path=os.path.join("results", "drifters.zarr"), - fieldset=config.drifter_fieldset, - drifters=drifters, + fieldset=config.drifter_config.fieldset, + drifters=schedule_results.measurements_to_simulate.drifters, outputdt=timedelta(hours=5), dt=timedelta(minutes=5), endtime=None, @@ -222,68 +89,348 @@ def sailship(config: VirtualShipConfig): print("Simulating argo floats") simulate_argo_floats( out_path=os.path.join("results", "argo_floats.zarr"), - argo_floats=argo_floats, + argo_floats=schedule_results.measurements_to_simulate.argo_floats, fieldset=config.argo_float_config.fieldset, outputdt=timedelta(minutes=5), endtime=None, ) - # convert CTD data to CSV - print("Postprocessing..") - postprocess() - - print("All data has been gathered and postprocessed, returning home.") - + # calculate cruise cost + assert ( + config.waypoints[0].time is not None + ), "First waypoints cannot have None time as this has been verified before during config verification." + time_past = schedule_results.end_spacetime.time - config.waypoints[0].time cost = costs(config, time_past) print(f"This cruise took {time_past} and would have cost {cost:,.0f} euros.") -def shiproute(config: VirtualShipConfig, dt: timedelta) -> list[Location]: +def _simulate_schedule( + waypoints: list[Waypoint], + projection: pyproj.Geod, + config: VirtualShipConfig, +) -> _ScheduleResults: """ - Take in route coordinates and return lat and lon points within region of interest to sample. + Simulate the sailing and aggregate the virtual measurements that should be taken. + :param waypoints: The waypoints. + :param projection: Projection used to sail between waypoints. :param config: The cruise configuration. - :param dt: Sailing time between each discrete route point. - :returns: lat and lon points within region of interest to sample. + :returns: Results from the simulation. + :raises NotImplementedError: When unsupported instruments are encountered. + :raises RuntimeError: When schedule appears infeasible. This should not happen in this version of virtual ship as the schedule is verified beforehand. """ - CRUISE_SPEED = 5.14 - - # discrete points the ship will pass - sample_points: list[Location] = [] - - # projection used to get discrete locations - geod = pyproj.Geod(ellps="WGS84") - - # loop over station coordinates and calculate intermediate points along great circle path - for startloc, endloc in zip(config.route_coordinates, config.route_coordinates[1:]): - # iterate over each coordinate and the next coordinate - # last coordinate has no next coordinate and is skipped - - # get locations between start and end, seperate by 5 minutes of cruising - # excludes final point, but this is added explicitly after this loop - int_points = geod.inv_intermediate( - startloc.lon, - startloc.lat, - endloc.lon, - endloc.lat, - del_s=CRUISE_SPEED * dt.total_seconds(), - initial_idx=0, - return_back_azimuth=False, + cruise = _Cruise(Spacetime(waypoints[0].location, waypoints[0].time)) + measurements = _MeasurementsToSimulate() + + # add recurring tasks to task list + waiting_tasks = SortedList[_WaitingTask]() + if config.ship_underwater_st_config is not None: + waiting_tasks.add( + _WaitingTask( + task=_ship_underwater_st_loop( + config.ship_underwater_st_config.period, cruise, measurements + ), + wait_until=cruise.spacetime.time, + ) + ) + if config.adcp_config is not None: + waiting_tasks.add( + _WaitingTask( + task=_adcp_loop(config.adcp_config.period, cruise, measurements), + wait_until=cruise.spacetime.time, + ) ) - sample_points.extend( - [ - Location(latitude=lat, longitude=lon) - for lat, lon in zip(int_points.lats, int_points.lons, strict=True) - ] + # sail to each waypoint while executing tasks + for waypoint in waypoints: + if waypoint.time is not None and cruise.spacetime.time > waypoint.time: + raise RuntimeError( + "Could not reach waypoint in time. This should not happen in this version of virtual ship as the schedule is verified beforehand." + ) + + # add task to the task queue for the instrument at the current waypoint + if waypoint.instrument is InstrumentType.ARGO_FLOAT: + _argo_float_task(cruise, measurements, config=config) + elif waypoint.instrument is InstrumentType.DRIFTER: + _drifter_task(cruise, measurements, config=config) + elif waypoint.instrument is InstrumentType.CTD: + waiting_tasks.add( + _WaitingTask( + _ctd_task( + config.ctd_config.stationkeeping_time, + config.ctd_config.min_depth, + config.ctd_config.max_depth, + cruise, + measurements, + ), + cruise.spacetime.time, + ) + ) + elif waypoint.instrument is None: + pass + else: + raise NotImplementedError() + + # sail to the next waypoint + waypoint_reached = False + while not waypoint_reached: + # execute all tasks planned for current time + while ( + len(waiting_tasks) > 0 + and waiting_tasks[0].wait_until <= cruise.spacetime.time + ): + task = waiting_tasks.pop(0) + try: + wait_for = next(task.task) + waiting_tasks.add( + _WaitingTask(task.task, cruise.spacetime.time + wait_for.time) + ) + except StopIteration: + pass + + # if sailing is prevented by a current task, just let time pass until the next task + if cruise.sail_is_locked: + cruise.spacetime = Spacetime( + cruise.spacetime.location, waiting_tasks[0].wait_until + ) + # else, let time pass while sailing + else: + # calculate time at which waypoint would be reached if simply sailing + geodinv: tuple[float, float, float] = projection.inv( + lons1=cruise.spacetime.location.lon, + lats1=cruise.spacetime.location.lat, + lons2=waypoint.location.lon, + lats2=waypoint.location.lat, + ) + azimuth1 = geodinv[0] + distance_to_next_waypoint = geodinv[2] + time_to_reach = timedelta( + seconds=distance_to_next_waypoint / config.ship_speed + ) + arrival_time = cruise.spacetime.time + time_to_reach + + # if waypoint is reached before next task, sail to the waypoint + if ( + len(waiting_tasks) == 0 + or arrival_time <= waiting_tasks[0].wait_until + ): + cruise.spacetime = Spacetime(waypoint.location, arrival_time) + waypoint_reached = True + # else, sail until task starts + else: + time_to_sail = waiting_tasks[0].wait_until - cruise.spacetime.time + distance_to_move = config.ship_speed * time_to_sail.total_seconds() + geodfwd: tuple[float, float, float] = projection.fwd( + lons=cruise.spacetime.location.lon, + lats=cruise.spacetime.location.lat, + az=azimuth1, + dist=distance_to_move, + ) + lon = geodfwd[0] + lat = geodfwd[1] + cruise.spacetime = Spacetime( + Location(latitude=lat, longitude=lon), + cruise.spacetime.time + time_to_sail, + ) + + cruise.finish() + + # don't sail anymore, but let tasks finish + while len(waiting_tasks) > 0: + task = waiting_tasks.pop(0) + try: + wait_for = next(task.task) + waiting_tasks.add( + _WaitingTask(task.task, cruise.spacetime.time + wait_for.time) + ) + except StopIteration: + pass + + return _ScheduleResults( + measurements_to_simulate=measurements, end_spacetime=cruise.spacetime + ) + + +class _Cruise: + _finished: bool # if last waypoint has been reached + _sail_lock_count: int # if sailing should be paused because of tasks; number of tasks that requested a pause; 0 means good to go sail + spacetime: Spacetime # current location and time + + def __init__(self, spacetime: Spacetime) -> None: + self._finished = False + self._sail_lock_count = 0 + self.spacetime = spacetime + + @property + def finished(self) -> bool: + return self._finished + + @contextmanager + def do_not_sail(self) -> Generator[None, None, None]: + try: + self._sail_lock_count += 1 + yield + finally: + self._sail_lock_count -= 1 + + def finish(self) -> None: + self._finished = True + + @property + def sail_is_locked(self) -> bool: + return self._sail_lock_count > 0 + + +@dataclass +class _MeasurementsToSimulate: + adcps: list[Spacetime] = field(default_factory=list, init=False) + ship_underwater_sts: list[Spacetime] = field(default_factory=list, init=False) + argo_floats: list[ArgoFloat] = field(default_factory=list, init=False) + drifters: list[Drifter] = field(default_factory=list, init=False) + ctds: list[CTD] = field(default_factory=list, init=False) + + +@dataclass +class _ScheduleResults: + measurements_to_simulate: _MeasurementsToSimulate + end_spacetime: Spacetime + + +@dataclass +class _WaitFor: + time: timedelta + + +class _WaitingTask: + _task: Generator[_WaitFor, None, None] + _wait_until: datetime + + def __init__( + self, task: Generator[_WaitFor, None, None], wait_until: datetime + ) -> None: + self._task = task + self._wait_until = wait_until + + def __lt__(self, other: _WaitingTask): + return self._wait_until < other._wait_until + + @property + def task(self) -> Generator[_WaitFor, None, None]: + return self._task + + @property + def wait_until(self) -> datetime: + return self._wait_until + + +def _ship_underwater_st_loop( + sample_period: timedelta, cruise: _Cruise, measurements: _MeasurementsToSimulate +) -> Generator[_WaitFor, None, None]: + while not cruise.finished: + measurements.ship_underwater_sts.append(cruise.spacetime) + yield _WaitFor(sample_period) + + +def _adcp_loop( + sample_period: timedelta, cruise: _Cruise, measurements: _MeasurementsToSimulate +) -> Generator[_WaitFor, None, None]: + while not cruise.finished: + measurements.adcps.append(cruise.spacetime) + yield _WaitFor(sample_period) + + +def _ctd_task( + stationkeeping_time: timedelta, + min_depth: float, + max_depth: float, + cruise: _Cruise, + measurements: _MeasurementsToSimulate, +) -> Generator[_WaitFor, None, None]: + with cruise.do_not_sail(): + measurements.ctds.append( + CTD( + spacetime=cruise.spacetime, + min_depth=min_depth, + max_depth=max_depth, + ) + ) + yield _WaitFor(stationkeeping_time) + + +def _drifter_task( + cruise: _Cruise, measurements: _MeasurementsToSimulate, config: VirtualShipConfig +) -> None: + measurements.drifters.append( + Drifter( + cruise.spacetime, + depth=config.drifter_config.depth, + lifetime=config.drifter_config.lifetime, ) + ) + - # explitly include final point which is not added by the previous loop - sample_points.append( - Location( - latitude=config.route_coordinates[-1].lat, - longitude=config.route_coordinates[-1].lon, +def _argo_float_task( + cruise: _Cruise, measurements: _MeasurementsToSimulate, config: VirtualShipConfig +) -> None: + measurements.argo_floats.append( + ArgoFloat( + spacetime=cruise.spacetime, + min_depth=config.argo_float_config.min_depth, + max_depth=config.argo_float_config.max_depth, + drift_depth=config.argo_float_config.drift_depth, + vertical_speed=config.argo_float_config.vertical_speed, + cycle_days=config.argo_float_config.cycle_days, + drift_days=config.argo_float_config.drift_days, ) ) - return sample_points + +def _verify_waypoints( + waypoints: list[Waypoint], ship_speed: float, projection: pyproj.Geod +) -> 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. + :raises PlanningError: If waypoints are not feasible or incorrect. + """ + # check first waypoint has a time + if 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] + if not all( + [ + next.time >= cur.time + for cur, next in zip(timed_waypoints, timed_waypoints[1:]) + ] + ): + raise PlanningError( + "Each waypoint should be timed after all previous 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:])): + if wp.instrument is InstrumentType.CTD: + time += timedelta(minutes=20) + + geodinv: tuple[float, float, float] = projection.inv( + wp.location.lon, wp.location.lat, wp_next.location.lon, wp_next.location.lat + ) + distance = geodinv[2] + + time_to_reach = timedelta(seconds=distance / ship_speed) + arrival_time = time + time_to_reach + + if wp_next.time is None: + time = arrival_time + elif arrival_time > wp_next.time: + raise PlanningError( + f"Waypoint planning is not valid: would arrive too late at a waypoint number {wp_i}. location: {wp.location} time: {wp.time} instrument: {wp.instrument}" + ) + else: + time = wp_next.time diff --git a/virtual_ship/virtual_ship_config.py b/virtual_ship/virtual_ship_config.py new file mode 100644 index 00000000..7dba815e --- /dev/null +++ b/virtual_ship/virtual_ship_config.py @@ -0,0 +1,117 @@ +"""VirtualShipConfig class.""" + +from dataclasses import dataclass +from datetime import timedelta + +from parcels import FieldSet + +from .location import Location +from .waypoint import Waypoint + + +@dataclass +class ArgoFloatConfig: + """Configuration for argos floats.""" + + fieldset: FieldSet + min_depth: float + max_depth: float + drift_depth: float + vertical_speed: float + cycle_days: float + drift_days: float + + +@dataclass +class ADCPConfig: + """Configuration for ADCP instrument.""" + + max_depth: float + bin_size_m: int + period: timedelta + fieldset: FieldSet + + +@dataclass +class CTDConfig: + """Configuration for CTD instrument.""" + + stationkeeping_time: timedelta + fieldset: FieldSet + min_depth: float + max_depth: float + + +@dataclass +class ShipUnderwaterSTConfig: + """Configuration for underwater ST.""" + + period: timedelta + fieldset: FieldSet + + +@dataclass +class DrifterConfig: + """Configuration for drifters.""" + + fieldset: FieldSet + depth: float + lifetime: timedelta + + +@dataclass +class VirtualShipConfig: + """Configuration of the virtual ship.""" + + ship_speed: float # m/s + + waypoints: list[Waypoint] + + argo_float_config: ArgoFloatConfig + adcp_config: ADCPConfig | None # if None, ADCP is disabled + ctd_config: CTDConfig + ship_underwater_st_config: ( + ShipUnderwaterSTConfig | None + ) # if None, ship underwater st is disabled + drifter_config: DrifterConfig + + def verify(self) -> None: + """ + Verify this configuration is valid. + + :raises ValueError: If not valid. + """ + if len(self.waypoints) < 2: + raise ValueError("Waypoints require at least a start and an end.") + + if not all( + [self._is_valid_location(waypoint.location) for waypoint in self.waypoints] + ): + raise ValueError("Invalid location for waypoint.") + + if self.argo_float_config.max_depth > 0: + raise ValueError("Argo float max depth must be negative or zero.") + + if self.argo_float_config.drift_depth > 0: + raise ValueError("Argo float drift depth must be negative or zero.") + + if self.argo_float_config.vertical_speed >= 0: + raise ValueError("Argo float vertical speed must be negative.") + + if self.argo_float_config.cycle_days <= 0: + raise ValueError("Argo float cycle days must be larger than zero.") + + if self.argo_float_config.drift_days <= 0: + raise ValueError("Argo drift cycle days must be larger than zero.") + + if self.adcp_config is not None and self.adcp_config.max_depth > 0: + raise ValueError("ADCP max depth must be negative.") + + @staticmethod + def _is_valid_location(location: Location) -> bool: + return ( + location.lat >= -90 + and location.lat <= 90 + and location.lon >= -180 + and location.lon <= 360 + ) diff --git a/virtual_ship/virtual_ship_configuration.py b/virtual_ship/virtual_ship_configuration.py deleted file mode 100644 index 8d48d4c0..00000000 --- a/virtual_ship/virtual_ship_configuration.py +++ /dev/null @@ -1,156 +0,0 @@ -"""VirtualShipConfig class.""" - -from dataclasses import dataclass -from datetime import datetime - -import numpy as np -from parcels import FieldSet - -from .location import Location - - -@dataclass -class ArgoFloatConfig: - """Configuration for argos floats.""" - - fieldset: FieldSet - max_depth: float - drift_depth: float - vertical_speed: float - cycle_days: float - drift_days: float - - -@dataclass -class ADCPConfig: - """Configuration for ADCP instrument.""" - - max_depth: float - bin_size_m: int - - -@dataclass -class VirtualShipConfig: - """Configuration of the virtual ship.""" - - start_time: datetime - route_coordinates: list[Location] - - adcp_fieldset: FieldSet - ship_underwater_st_fieldset: FieldSet - ctd_fieldset: FieldSet - drifter_fieldset: FieldSet - - argo_float_deploy_locations: list[Location] - drifter_deploy_locations: list[Location] - ctd_deploy_locations: list[Location] - - argo_float_config: ArgoFloatConfig - adcp_config: ADCPConfig - - def verify(self) -> None: - """ - Verify this configuration is valid. - - :raises ValueError: If not valid. - """ - if len(self.route_coordinates) < 2: - raise ValueError("Route needs to consist of at least locations.") - - if not all( - [self._is_valid_location(coord) for coord in self.route_coordinates] - ): - raise ValueError("Invalid coordinates in route.") - - if not all( - [ - self._is_valid_location(coord) - for coord in self.argo_float_deploy_locations - ] - ): - raise ValueError("Argo float deploy locations are not valid coordinates.") - - if not all( - [ - any( - [ - np.isclose(deploy.lat, coord.lat) - and np.isclose(deploy.lon, coord.lon) - for coord in self.route_coordinates - ] - ) - for deploy in self.argo_float_deploy_locations - ] - ): - raise ValueError( - "Argo float deploy locations are not exactly on route coordinates." - ) - - if not all( - [self._is_valid_location(coord) for coord in self.drifter_deploy_locations] - ): - raise ValueError("Drifter deploy locations are not valid coordinates.") - - if not all( - [ - any( - [ - np.isclose(deploy.lat, coord.lat) - and np.isclose(deploy.lon, coord.lon) - for coord in self.route_coordinates - ] - ) - for deploy in self.drifter_deploy_locations - ] - ): - raise ValueError( - "Drifter deploy locations are not exactly on route coordinates." - ) - - if not all( - [self._is_valid_location(coord) for coord in self.ctd_deploy_locations] - ): - raise ValueError("CTD deploy locations are not valid coordinates.") - - if not all( - [ - any( - [ - np.isclose(deploy.lat, coord.lat) - and np.isclose(deploy.lon, coord.lon) - for coord in self.route_coordinates - ] - ) - for deploy in self.ctd_deploy_locations - ] - ): - raise ValueError( - "CTD deploy locations are not exactly on route coordinates." - ) - - if self.argo_float_config.max_depth > 0: - raise ValueError("Argo float max depth must be negative or zero.") - - if self.argo_float_config.drift_depth > 0: - raise ValueError("Argo float drift depth must be negative or zero.") - - if self.argo_float_config.vertical_speed >= 0: - raise ValueError("Argo float vertical speed must be negative.") - - if self.argo_float_config.cycle_days <= 0: - raise ValueError("Argo float cycle days must be larger than zero.") - - if self.argo_float_config.drift_days <= 0: - raise ValueError("Argo drift cycle days must be larger than zero.") - - if self.adcp_config.max_depth > 0: - raise ValueError("ADCP max depth must be negative.") - - @staticmethod - def _is_valid_location(location: Location) -> bool: - return ( - location.lat >= -90 - and location.lat <= 90 - and location.lon >= -180 - and location.lon <= 360 - ) diff --git a/virtual_ship/waypoint.py b/virtual_ship/waypoint.py new file mode 100644 index 00000000..bf3d33b4 --- /dev/null +++ b/virtual_ship/waypoint.py @@ -0,0 +1,16 @@ +"""Waypoint class.""" + +from dataclasses import dataclass +from datetime import datetime + +from .instrument_type import InstrumentType +from .location import Location + + +@dataclass +class Waypoint: + """A Waypoint to sail to with an optional time and an optional instrument.""" + + location: Location + time: datetime | None = None + instrument: InstrumentType | None = None