diff --git a/meta.yaml b/meta.yaml index f38a2c79..6073e7c7 100644 --- a/meta.yaml +++ b/meta.yaml @@ -5,6 +5,10 @@ package: source: path: virtual_ship +build: + entry_points: + - do_expedition = virtual_ship.cli.do_expedition:main + requirements: run: - python >=3.8 diff --git a/pyproject.toml b/pyproject.toml index d87416a4..f585c3f7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,7 @@ dependencies = [ "sortedcontainers == 2.4.0", "opensimplex == 0.4.5", "numpy >=1, < 2", + "pydantic >=2, <3", ] [project.urls] @@ -61,6 +62,9 @@ dev = [ "sortedcontainers-stubs == 2.4.2", ] +[project.scripts] +do_expedition = "virtual_ship.cli.do_expedition:main" + [tool.isort] profile = "black" skip_gitignore = true diff --git a/scripts/download_data.py b/scripts/download_data.py new file mode 100644 index 00000000..d89edba8 --- /dev/null +++ b/scripts/download_data.py @@ -0,0 +1,136 @@ +""" +Download data required to run expeditions. + +This is a very crude script, here just as long as we do not properly incorporate it into the library. +""" + +import copernicusmarine +import datetime + +if __name__ == "__main__": + datadir = "input_data" + username = input("username: ") + password = input("password: ") + + copernicusmarine.subset( + dataset_id="cmems_mod_glo_phy_my_0.083deg_static", + force_dataset_part="bathy", + variables=["deptho"], + minimum_longitude=-0.01, + maximum_longitude=0.01, + minimum_latitude=-0.01, + maximum_latitude=0.01, + minimum_depth=0.49402499198913574, + maximum_depth=5727.9169921875, + output_filename="bathymetry.nc", + output_directory=datadir, + username=username, + password=password, + force_download=True, + ) + + download_dict = { + "UVdata": { + "dataset_id": "cmems_mod_glo_phy-cur_anfc_0.083deg_PT6H-i", + "variables": ["uo", "vo"], + "output_filename": "default_uv.nc", + }, + "Sdata": { + "dataset_id": "cmems_mod_glo_phy-so_anfc_0.083deg_PT6H-i", + "variables": ["so"], + "output_filename": "default_s.nc", + }, + "Tdata": { + "dataset_id": "cmems_mod_glo_phy-thetao_anfc_0.083deg_PT6H-i", + "variables": ["thetao"], + "output_filename": "default_t.nc", + }, + } + + for dataset in download_dict: + copernicusmarine.subset( + dataset_id=download_dict[dataset]["dataset_id"], + variables=download_dict[dataset]["variables"], + minimum_longitude=-0.01, + maximum_longitude=0.01, + minimum_latitude=-0.01, + maximum_latitude=0.01, + start_datetime=datetime.datetime.strptime("2023-01-01", "%Y-%m-%d"), + end_datetime=datetime.datetime.strptime("2023-01-02", "%Y-%m-%d"), + minimum_depth=0.49402499198913574, + maximum_depth=5727.9169921875, + output_filename=download_dict[dataset]["output_filename"], + output_directory=datadir, + username=username, + password=password, + force_download=True, + ) + + download_dict = { + "UVdata": { + "dataset_id": "cmems_mod_glo_phy-cur_anfc_0.083deg_PT6H-i", + "variables": ["uo", "vo"], + "output_filename": "drifter_uv.nc", + }, + "Tdata": { + "dataset_id": "cmems_mod_glo_phy-thetao_anfc_0.083deg_PT6H-i", + "variables": ["thetao"], + "output_filename": "drifter_t.nc", + }, + } + + for dataset in download_dict: + copernicusmarine.subset( + dataset_id=download_dict[dataset]["dataset_id"], + variables=download_dict[dataset]["variables"], + minimum_longitude=-0.01, + maximum_longitude=0.01, + minimum_latitude=-0.01, + maximum_latitude=0.01, + start_datetime=datetime.datetime.strptime("2023-01-01", "%Y-%m-%d"), + end_datetime=datetime.datetime.strptime("2023-01-02", "%Y-%m-%d"), + minimum_depth=0.49402499198913574, + maximum_depth=0.49402499198913574, + output_filename=download_dict[dataset]["output_filename"], + output_directory=datadir, + username=username, + password=password, + force_download=True, + ) + + download_dict = { + "UVdata": { + "dataset_id": "cmems_mod_glo_phy-cur_anfc_0.083deg_PT6H-i", + "variables": ["uo", "vo"], + "output_filename": "argo_float_uv.nc", + }, + "Sdata": { + "dataset_id": "cmems_mod_glo_phy-so_anfc_0.083deg_PT6H-i", + "variables": ["so"], + "output_filename": "argo_float_s.nc", + }, + "Tdata": { + "dataset_id": "cmems_mod_glo_phy-thetao_anfc_0.083deg_PT6H-i", + "variables": ["thetao"], + "output_filename": "argo_float_t.nc", + }, + } + + for dataset in download_dict: + copernicusmarine.subset( + dataset_id=download_dict[dataset]["dataset_id"], + variables=download_dict[dataset]["variables"], + minimum_longitude=-0.01, + maximum_longitude=0.01, + minimum_latitude=-0.01, + maximum_latitude=0.01, + start_datetime=datetime.datetime.strptime("2023-01-01", "%Y-%m-%d"), + end_datetime=datetime.datetime.strptime("2023-01-02", "%Y-%m-%d"), + minimum_depth=0.49402499198913574, + maximum_depth=5727.9169921875, + output_filename=download_dict[dataset]["output_filename"], + output_directory=datadir, + username=username, + password=password, + force_download=True, + ) diff --git a/tests/expedition/expedition_dir/.gitignore b/tests/expedition/expedition_dir/.gitignore new file mode 100644 index 00000000..fbca2253 --- /dev/null +++ b/tests/expedition/expedition_dir/.gitignore @@ -0,0 +1 @@ +results/ diff --git a/tests/expedition/expedition_dir/input_data/.gitignore b/tests/expedition/expedition_dir/input_data/.gitignore new file mode 100644 index 00000000..0f4bf3bf --- /dev/null +++ b/tests/expedition/expedition_dir/input_data/.gitignore @@ -0,0 +1,9 @@ +!argo_float_s.nc +!argo_float_t.nc +!argo_float_uv.nc +!bathymetry.nc +!default_s.nc +!default_t.nc +!default_uv.nc +!drifter_t.nc +!drifter_uv.nc diff --git a/tests/expedition/expedition_dir/input_data/argo_float_s.nc b/tests/expedition/expedition_dir/input_data/argo_float_s.nc new file mode 100644 index 00000000..7743124b Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/argo_float_s.nc differ diff --git a/tests/expedition/expedition_dir/input_data/argo_float_t.nc b/tests/expedition/expedition_dir/input_data/argo_float_t.nc new file mode 100644 index 00000000..cee4478d Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/argo_float_t.nc differ diff --git a/tests/expedition/expedition_dir/input_data/argo_float_uv.nc b/tests/expedition/expedition_dir/input_data/argo_float_uv.nc new file mode 100644 index 00000000..6012603a Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/argo_float_uv.nc differ diff --git a/tests/expedition/expedition_dir/input_data/bathymetry.nc b/tests/expedition/expedition_dir/input_data/bathymetry.nc new file mode 100644 index 00000000..002eb339 Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/bathymetry.nc differ diff --git a/tests/expedition/expedition_dir/input_data/default_s.nc b/tests/expedition/expedition_dir/input_data/default_s.nc new file mode 100644 index 00000000..7743124b Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/default_s.nc differ diff --git a/tests/expedition/expedition_dir/input_data/default_t.nc b/tests/expedition/expedition_dir/input_data/default_t.nc new file mode 100644 index 00000000..cee4478d Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/default_t.nc differ diff --git a/tests/expedition/expedition_dir/input_data/default_uv.nc b/tests/expedition/expedition_dir/input_data/default_uv.nc new file mode 100644 index 00000000..908cebaf Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/default_uv.nc differ diff --git a/tests/expedition/expedition_dir/input_data/drifter_t.nc b/tests/expedition/expedition_dir/input_data/drifter_t.nc new file mode 100644 index 00000000..83cad2cf Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/drifter_t.nc differ diff --git a/tests/expedition/expedition_dir/input_data/drifter_uv.nc b/tests/expedition/expedition_dir/input_data/drifter_uv.nc new file mode 100644 index 00000000..5a14b9d4 Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/drifter_uv.nc differ diff --git a/tests/expedition/expedition_dir/schedule.yaml b/tests/expedition/expedition_dir/schedule.yaml new file mode 100644 index 00000000..6d55f6c5 --- /dev/null +++ b/tests/expedition/expedition_dir/schedule.yaml @@ -0,0 +1,16 @@ +waypoints: +- instrument: CTD + location: + latitude: 0 + longitude: 0 + time: 2023-01-01 00:00:00 +- instrument: DRIFTER + location: + latitude: 0.01 + longitude: 0.01 + time: 2023-01-01 01:00:00 +- instrument: ARGO_FLOAT + location: + latitude: 0.02 + longitude: 0.02 + time: 2023-01-01 02:00:00 diff --git a/tests/expedition/expedition_dir/ship_config.yaml b/tests/expedition/expedition_dir/ship_config.yaml new file mode 100644 index 00000000..c057d6b5 --- /dev/null +++ b/tests/expedition/expedition_dir/ship_config.yaml @@ -0,0 +1,21 @@ +ship_speed_meter_per_second: 5.14 +adcp_config: + num_bins: 40 + max_depth_meter: -1000.0 + period_minutes: 5.0 +argo_float_config: + cycle_days: 10.0 + drift_days: 9.0 + drift_depth_meter: -1000.0 + max_depth_meter: -2000.0 + min_depth_meter: 0.0 + vertical_speed_meter_per_second: -0.1 +ctd_config: + max_depth_meter: -2000.0 + min_depth_meter: -11.0 + stationkeeping_time_minutes: 20.0 +drifter_config: + depth_meter: 0.0 + lifetime_minutes: 40320.0 +ship_underwater_st_config: + period_minutes: 5.0 diff --git a/tests/expedition/test_do_expedition.py b/tests/expedition/test_do_expedition.py new file mode 100644 index 00000000..11e6b862 --- /dev/null +++ b/tests/expedition/test_do_expedition.py @@ -0,0 +1,9 @@ +from pytest import CaptureFixture + +from virtual_ship.expedition import do_expedition + + +def test_do_expedition(capfd: CaptureFixture) -> None: + do_expedition("expedition_dir") + out, _ = capfd.readouterr() + assert "This expedition took" in out, "Expedition did not complete successfully." diff --git a/tests/expedition/test_schedule.py b/tests/expedition/test_schedule.py new file mode 100644 index 00000000..37ee4a47 --- /dev/null +++ b/tests/expedition/test_schedule.py @@ -0,0 +1,26 @@ +from datetime import datetime, timedelta + +import py + +from virtual_ship import Location +from virtual_ship.expedition import Schedule, Waypoint + + +def test_schedule(tmpdir: py.path.LocalPath) -> None: + out_path = tmpdir.join("schedule.yaml") + + # arbitrary time for testing + base_time = datetime.strptime("1950-01-01", "%Y-%m-%d") + + schedule = Schedule( + waypoints=[ + Waypoint(Location(0, 0), time=base_time, instrument=None), + Waypoint( + Location(1, 1), time=base_time + timedelta(hours=1), instrument=None + ), + ] + ) + schedule.to_yaml(out_path) + + schedule2 = Schedule.from_yaml(out_path) + assert schedule == schedule2 diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py new file mode 100644 index 00000000..e27fed3f --- /dev/null +++ b/tests/expedition/test_simulate_schedule.py @@ -0,0 +1,48 @@ +from datetime import datetime, timedelta + +import pyproj + +from virtual_ship import Location +from virtual_ship.expedition import Schedule, ShipConfig, Waypoint +from virtual_ship.expedition.simulate_schedule import ( + ScheduleOk, + ScheduleProblem, + simulate_schedule, +) + + +def test_simulate_schedule_feasible() -> None: + """Test schedule with two waypoints that can be reached within time is OK.""" + base_time = datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S") + + projection = pyproj.Geod(ellps="WGS84") + ship_config = ShipConfig.from_yaml("expedition_dir/ship_config.yaml") + ship_config.ship_speed_meter_per_second = 5.14 + schedule = Schedule( + waypoints=[ + Waypoint(Location(0, 0), base_time), + Waypoint(Location(0.01, 0), base_time + timedelta(days=1)), + ] + ) + + result = simulate_schedule(projection, ship_config, schedule) + + assert isinstance(result, ScheduleOk) + + +def test_simulate_schedule_too_far() -> None: + """Test schedule with two waypoints that are very far away and cannot be reached in time is not OK.""" + base_time = datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S") + + projection = pyproj.Geod(ellps="WGS84") + ship_config = ShipConfig.from_yaml("expedition_dir/ship_config.yaml") + schedule = Schedule( + waypoints=[ + Waypoint(Location(0, 0), base_time), + Waypoint(Location(1.0, 0), base_time + timedelta(minutes=1)), + ] + ) + + result = simulate_schedule(projection, ship_config, schedule) + + assert isinstance(result, ScheduleProblem) diff --git a/tests/instruments/test_ship_underwater_st.py b/tests/instruments/test_ship_underwater_st.py index a0535c8d..54b8bdf0 100644 --- a/tests/instruments/test_ship_underwater_st.py +++ b/tests/instruments/test_ship_underwater_st.py @@ -27,15 +27,15 @@ def test_simulate_ship_underwater_st(tmpdir: py.path.LocalPath) -> None: # expected observations at sample points expected_obs = [ { - "salinity": 5, - "temperature": 6, + "S": 5, + "T": 6, "lat": sample_points[0].location.lat, "lon": sample_points[0].location.lon, "time": base_time + datetime.timedelta(seconds=0), }, { - "salinity": 7, - "temperature": 8, + "S": 7, + "T": 8, "lat": sample_points[1].location.lat, "lon": sample_points[1].location.lon, "time": base_time + datetime.timedelta(seconds=1), @@ -45,19 +45,19 @@ def test_simulate_ship_underwater_st(tmpdir: py.path.LocalPath) -> None: # create fieldset based on the expected observations # indices are time, latitude, longitude salinity = np.zeros((2, 2, 2)) - salinity[0, 0, 0] = expected_obs[0]["salinity"] - salinity[1, 1, 1] = expected_obs[1]["salinity"] + salinity[0, 0, 0] = expected_obs[0]["S"] + salinity[1, 1, 1] = expected_obs[1]["S"] temperature = np.zeros((2, 2, 2)) - temperature[0, 0, 0] = expected_obs[0]["temperature"] - temperature[1, 1, 1] = expected_obs[1]["temperature"] + temperature[0, 0, 0] = expected_obs[0]["T"] + temperature[1, 1, 1] = expected_obs[1]["T"] fieldset = FieldSet.from_data( { "V": np.zeros((2, 2, 2)), "U": np.zeros((2, 2, 2)), - "salinity": salinity, - "temperature": temperature, + "S": salinity, + "T": temperature, }, { "lat": np.array([expected_obs[0]["lat"], expected_obs[1]["lat"]]), @@ -95,7 +95,7 @@ def test_simulate_ship_underwater_st(tmpdir: py.path.LocalPath) -> None: zip(results.sel(trajectory=traj).obs, expected_obs, strict=True) ): obs = results.sel(trajectory=traj, obs=obs_i) - for var in ["salinity", "temperature", "lat", "lon"]: + for var in ["S", "T", "lat", "lon"]: obs_value = obs[var].values.item() exp_value = exp[var] assert np.isclose( diff --git a/tests/test_sailship.py b/tests/test_sailship.py deleted file mode 100644 index 0cdcd521..00000000 --- a/tests/test_sailship.py +++ /dev/null @@ -1,236 +0,0 @@ -"""Performs a complete cruise with virtual ship.""" - -import datetime -from datetime import timedelta - -import numpy as np -import pyproj -import pytest -from parcels import Field, FieldSet - -from virtual_ship import InstrumentType, Location, Schedule, Waypoint -from virtual_ship.sailship import PlanningError, _verify_waypoints, sailship -from virtual_ship.virtual_ship_config import ( - ADCPConfig, - ArgoFloatConfig, - CTDConfig, - DrifterConfig, - ShipUnderwaterSTConfig, - VirtualShipConfig, -) - - -def _make_ctd_fieldset(base_time: datetime) -> FieldSet: - 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}, - { - "time": [ - np.datetime64(base_time + datetime.timedelta(seconds=0)), - np.datetime64(base_time + datetime.timedelta(minutes=200)), - ], - "depth": [0, -1000], - "lat": [-40, 90], - "lon": [-90, 90], - }, - ) - fieldset.add_field(Field("bathymetry", [-1000], lon=0, lat=0)) - return fieldset - - -def _make_drifter_fieldset(base_time: datetime) -> FieldSet: - v = np.full((2, 2, 2), 1.0) - u = np.full((2, 2, 2), 1.0) - t = np.full((2, 2, 2), 1.0) - - fieldset = FieldSet.from_data( - {"V": v, "U": u, "T": t}, - { - "time": [ - np.datetime64(base_time + datetime.timedelta(seconds=0)), - np.datetime64(base_time + datetime.timedelta(weeks=10)), - ], - "lat": [-40, 90], - "lon": [-90, 90], - }, - ) - return fieldset - - -def test_sailship() -> None: - # arbitrary time offset for the dummy fieldsets - base_time = datetime.datetime.strptime("2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S") - - adcp_fieldset = FieldSet.from_data( - {"U": 1, "V": 1}, - {"lon": 0, "lat": 0}, - ) - - ship_underwater_st_fieldset = FieldSet.from_data( - {"U": 1, "V": 1, "salinity": 0, "temperature": 0}, - {"lon": 0, "lat": 0}, - ) - - ctd_fieldset = _make_ctd_fieldset(base_time) - - drifter_fieldset = _make_drifter_fieldset(base_time) - - argo_float_fieldset = FieldSet.from_data( - {"U": 1, "V": 1, "T": 0, "S": 0}, - { - "lon": 0, - "lat": 0, - "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], - }, - ) - - 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, - cycle_days=10, - drift_days=9, - ) - - adcp_config = ADCPConfig( - max_depth=-1000, - bin_size_m=24, - period=timedelta(minutes=5), - fieldset=adcp_fieldset, - ) - - 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, - ), - Waypoint( - location=Location(latitude=-23.381289, longitude=63.743631), - instrument=InstrumentType.ARGO_FLOAT, - ), - ] - - schedule = Schedule(waypoints=waypoints) - - config = VirtualShipConfig( - ship_speed=5.14, - schedule=schedule, - 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) - - -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, - schedule=Schedule(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/tests/test_schedule.py b/tests/test_schedule.py deleted file mode 100644 index 371c786f..00000000 --- a/tests/test_schedule.py +++ /dev/null @@ -1,18 +0,0 @@ -import py - -from virtual_ship import Location, Schedule, Waypoint - - -def test_schedule(tmpdir: py.path.LocalPath) -> None: - out_path = tmpdir.join("schedule.yaml") - - schedule = Schedule( - [ - Waypoint(Location(0, 0), time=0, instrument=None), - Waypoint(Location(1, 1), time=1, instrument=None), - ] - ) - schedule.to_yaml(out_path) - - schedule2 = Schedule.from_yaml(out_path) - assert schedule == schedule2 diff --git a/virtual_ship/__init__.py b/virtual_ship/__init__.py index 717f5c38..6e74f1c8 100644 --- a/virtual_ship/__init__.py +++ b/virtual_ship/__init__.py @@ -1,20 +1,9 @@ """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 .schedule import Schedule from .spacetime import Spacetime -from .waypoint import Waypoint __all__ = [ - "InstrumentType", "Location", - "PlanningError", - "Schedule", "Spacetime", - "Waypoint", - "instruments", - "sailship", ] diff --git a/virtual_ship/cli/__init__.py b/virtual_ship/cli/__init__.py new file mode 100644 index 00000000..d9f62d70 --- /dev/null +++ b/virtual_ship/cli/__init__.py @@ -0,0 +1 @@ +"""Command line interface tools.""" diff --git a/virtual_ship/cli/do_expedition.py b/virtual_ship/cli/do_expedition.py new file mode 100644 index 00000000..868a1aa3 --- /dev/null +++ b/virtual_ship/cli/do_expedition.py @@ -0,0 +1,30 @@ +""" +Command line interface tool for virtualship.expedition.do_expedition:do_expedition function. + +See --help for usage. +""" + +import argparse +from pathlib import Path + +from virtual_ship.expedition.do_expedition import do_expedition + + +def main() -> None: + """Entrypoint for the tool.""" + parser = argparse.ArgumentParser( + prog="do_expedition", + description="Perform an expedition based on a provided schedule.", + ) + parser.add_argument( + "dir", + type=str, + help="Directory for the expedition. This should contain all required configuration files, and the result will be saved here as well.", + ) + args = parser.parse_args() + + do_expedition(Path(args.dir)) + + +if __name__ == "__main__": + main() diff --git a/virtual_ship/costs.py b/virtual_ship/costs.py deleted file mode 100644 index 9e4c20e6..00000000 --- a/virtual_ship/costs.py +++ /dev/null @@ -1,40 +0,0 @@ -"""costs function.""" - -from datetime import timedelta - -from .instrument_type import InstrumentType -from .virtual_ship_config import VirtualShipConfig - - -def costs(config: VirtualShipConfig, total_time: timedelta): - """ - Calculate the cost of the virtual ship (in US$). - - :param config: The cruise configuration. - :param total_time: Time cruised. - :returns: The calculated cost of the cruise. - """ - ship_cost_per_day = 30000 - drifter_deploy_cost = 2500 - argo_deploy_cost = 15000 - - ship_cost = ship_cost_per_day / 24 * total_time.total_seconds() // 3600 - num_argos = len( - [ - waypoint - for waypoint in config.schedule.waypoints - if waypoint.instrument is InstrumentType.ARGO_FLOAT - ] - ) - argo_cost = num_argos * argo_deploy_cost - num_drifters = len( - [ - waypoint - for waypoint in config.schedule.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/expedition/__init__.py b/virtual_ship/expedition/__init__.py new file mode 100644 index 00000000..c755d33b --- /dev/null +++ b/virtual_ship/expedition/__init__.py @@ -0,0 +1,28 @@ +"""Everything for simulating an expedition.""" + +from .do_expedition import do_expedition +from .instrument_type import InstrumentType +from .schedule import Schedule +from .ship_config import ( + ADCPConfig, + ArgoFloatConfig, + CTDConfig, + DrifterConfig, + ShipConfig, + ShipUnderwaterSTConfig, +) +from .waypoint import Waypoint + +__all__ = [ + "ADCPConfig", + "ArgoFloatConfig", + "CTDConfig", + "DrifterConfig", + "InstrumentType", + "Schedule", + "ShipConfig", + "ShipUnderwaterSTConfig", + "Waypoint", + "do_expedition", + "instruments", +] diff --git a/virtual_ship/expedition/checkpoint.py b/virtual_ship/expedition/checkpoint.py new file mode 100644 index 00000000..d3bb08a4 --- /dev/null +++ b/virtual_ship/expedition/checkpoint.py @@ -0,0 +1,51 @@ +"""Checkpoint class.""" + +from __future__ import annotations + +from pathlib import Path + +import pydantic +import yaml + +from .instrument_type import InstrumentType +from .schedule import Schedule + + +class _YamlDumper(yaml.SafeDumper): + pass + + +_YamlDumper.add_representer( + InstrumentType, lambda dumper, data: dumper.represent_data(data.value) +) + + +class Checkpoint(pydantic.BaseModel): + """ + A checkpoint of schedule simulation. + + Copy of the schedule until where the simulation proceeded without troubles. + """ + + past_schedule: Schedule + + def to_yaml(self, file_path: str | Path) -> None: + """ + Write checkpoint to yaml file. + + :param file_path: Path to the file to write to. + """ + with open(file_path, "w") as file: + yaml.dump(self.model_dump(by_alias=True), file, Dumper=_YamlDumper) + + @classmethod + def from_yaml(cls, file_path: str | Path) -> Checkpoint: + """ + Load checkpoint from yaml file. + + :param file_path: Path to the file to load from. + :returns: The checkpoint. + """ + with open(file_path, "r") as file: + data = yaml.safe_load(file) + return Checkpoint(**data) diff --git a/virtual_ship/expedition/do_expedition.py b/virtual_ship/expedition/do_expedition.py new file mode 100644 index 00000000..0b43ee8e --- /dev/null +++ b/virtual_ship/expedition/do_expedition.py @@ -0,0 +1,150 @@ +"""do_expedition function.""" + +import os +import shutil +from pathlib import Path + +import pyproj + +from .checkpoint import Checkpoint +from .expedition_cost import expedition_cost +from .input_data import InputData +from .schedule import Schedule +from .ship_config import ShipConfig +from .simulate_measurements import simulate_measurements +from .simulate_schedule import ScheduleProblem, simulate_schedule +from .verify_schedule import verify_schedule + + +def do_expedition(expedition_dir: str | Path) -> None: + """ + Perform an expedition, providing terminal feedback and file output. + + :param expedition_dir: The base directory for the expedition. + """ + if isinstance(expedition_dir, str): + expedition_dir = Path(expedition_dir) + + # load ship configuration + ship_config = _get_ship_config(expedition_dir) + if ship_config is None: + return + + # load schedule + schedule = _get_schedule(expedition_dir) + if schedule is None: + return + + # load last checkpoint + checkpoint = _load_checkpoint(expedition_dir) + if checkpoint is None: + checkpoint = Checkpoint(past_schedule=Schedule(waypoints=[])) + + # verify that schedule and checkpoint match + if ( + not schedule.waypoints[: len(checkpoint.past_schedule.waypoints)] + == checkpoint.past_schedule.waypoints + ): + print( + "Past waypoints in schedule have been changed! Restore past schedule and only change future waypoints." + ) + return + + # projection used to sail between waypoints + projection = pyproj.Geod(ellps="WGS84") + + # load fieldsets + input_data = _load_input_data( + expedition_dir=expedition_dir, ship_config=ship_config + ) + + # verify schedule makes sense + verify_schedule(projection, ship_config, schedule, input_data) + + # simulate the schedule + schedule_results = simulate_schedule( + projection=projection, ship_config=ship_config, schedule=schedule + ) + if isinstance(schedule_results, ScheduleProblem): + print( + "Update your schedule and continue the expedition by running the tool again." + ) + _save_checkpoint( + Checkpoint( + past_schedule=Schedule( + waypoints=schedule.waypoints[: schedule_results.failed_waypoint_i] + ) + ), + expedition_dir, + ) + return + + # delete and create results directory + if os.path.exists(expedition_dir.joinpath("results")): + shutil.rmtree(expedition_dir.joinpath("results")) + os.makedirs(expedition_dir.joinpath("results")) + + # calculate expedition cost in US$ + assert ( + schedule.waypoints[0].time is not None + ), "First waypoint has no time. This should not be possible as it should have been verified before." + time_past = schedule_results.time - schedule.waypoints[0].time + cost = expedition_cost(schedule_results, time_past) + with open(expedition_dir.joinpath("results", "cost.txt"), "w") as file: + file.writelines(f"cost: {cost} US$") + print(f"This expedition took {time_past} and would have cost {cost:,.0f} US$.") + + # simulate measurements + print("Simulating measurements. This may take a while..") + simulate_measurements( + expedition_dir, + ship_config, + input_data, + schedule_results.measurements_to_simulate, + ) + print("Done simulating measurements.") + + print("Your expedition has concluded successfully!") + print("Your measurements can be found in the results directory.") + + +def _get_ship_config(expedition_dir: Path) -> ShipConfig | None: + file_path = expedition_dir.joinpath("ship_config.yaml") + try: + return ShipConfig.from_yaml(file_path) + except FileNotFoundError: + print(f'Schedule not found. Save it to "{file_path}".') + return None + + +def _load_input_data(expedition_dir: Path, ship_config: ShipConfig) -> InputData: + return InputData.load( + directory=expedition_dir.joinpath("input_data"), + load_adcp=ship_config.adcp_config is not None, + load_argo_float=ship_config.argo_float_config is not None, + load_ctd=ship_config.ctd_config is not None, + load_drifter=ship_config.drifter_config is not None, + load_ship_underwater_st=ship_config.ship_underwater_st_config is not None, + ) + + +def _get_schedule(expedition_dir: Path) -> Schedule | None: + file_path = expedition_dir.joinpath("schedule.yaml") + try: + return Schedule.from_yaml(file_path) + except FileNotFoundError: + print(f'Schedule not found. Save it to "{file_path}".') + return None + + +def _load_checkpoint(expedition_dir: Path) -> Checkpoint | None: + file_path = expedition_dir.joinpath("checkpoint.yaml") + try: + return Checkpoint.from_yaml(file_path) + except FileNotFoundError: + return None + + +def _save_checkpoint(checkpoint: Checkpoint, expedition_dir: Path) -> None: + file_path = expedition_dir.joinpath("checkpoint.yaml") + checkpoint.to_yaml(file_path) diff --git a/virtual_ship/expedition/expedition_cost.py b/virtual_ship/expedition/expedition_cost.py new file mode 100644 index 00000000..cab6ab7d --- /dev/null +++ b/virtual_ship/expedition/expedition_cost.py @@ -0,0 +1,27 @@ +"""expedition_cost function.""" + +from datetime import timedelta + +from .simulate_schedule import ScheduleOk + + +def expedition_cost(schedule_results: ScheduleOk, time_past: timedelta) -> float: + """ + Calculate the cost of the expedition in US$. + + :param schedule_results: Results from schedule simulation. + :param time_past: Time the expedition took. + :returns: The calculated cost of the expedition in US$. + """ + SHIP_COST_PER_DAY = 30000 + DRIFTER_DEPLOY_COST = 2500 + ARGO_DEPLOY_COST = 15000 + + ship_cost = SHIP_COST_PER_DAY / 24 * time_past.total_seconds() // 3600 + num_argos = len(schedule_results.measurements_to_simulate.argo_floats) + argo_cost = num_argos * ARGO_DEPLOY_COST + num_drifters = len(schedule_results.measurements_to_simulate.drifters) + drifter_cost = num_drifters * DRIFTER_DEPLOY_COST + + cost = ship_cost + argo_cost + drifter_cost + return cost diff --git a/virtual_ship/expedition/input_data.py b/virtual_ship/expedition/input_data.py new file mode 100644 index 00000000..7af9ef72 --- /dev/null +++ b/virtual_ship/expedition/input_data.py @@ -0,0 +1,177 @@ +"""InputData class.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +from parcels import Field, FieldSet + + +@dataclass +class InputData: + """A collection of fieldsets that function as input data for simulation.""" + + adcp_fieldset: FieldSet | None + argo_float_fieldset: FieldSet | None + ctd_fieldset: FieldSet | None + drifter_fieldset: FieldSet | None + ship_underwater_st_fieldset: FieldSet | None + + @classmethod + def load( + cls, + directory: str | Path, + load_adcp: bool, + load_argo_float: bool, + load_ctd: bool, + load_drifter: bool, + load_ship_underwater_st: bool, + ) -> InputData: + """ + Create an instance of this class from netCDF files. + + For now this function makes a lot of assumption about file location and contents. + + :param directory: Base directory of the expedition. + :param load_adcp: Whether to load the ADCP fieldset. + :param load_argo_float: Whether to load the argo float fieldset. + :param load_ctd: Whether to load the CTD fieldset. + :param load_drifter: Whether to load the drifter fieldset. + :param load_ship_underwater_st: Whether to load the ship underwater ST fieldset. + :returns: An instance of this class with loaded fieldsets. + """ + if load_drifter: + drifter_fieldset = cls._load_drifter_fieldset(directory) + else: + drifter_fieldset = None + if load_argo_float: + argo_float_fieldset = cls._load_argo_float_fieldset(directory) + else: + argo_float_fieldset = None + if load_adcp or load_ctd or load_ship_underwater_st: + default_fieldset = cls._load_default_fieldset(directory) + if load_adcp: + adcp_fieldset = default_fieldset + else: + adcp_fieldset = None + if load_ctd: + ctd_fieldset = default_fieldset + else: + ctd_fieldset = None + if load_ship_underwater_st: + ship_underwater_st_fieldset = default_fieldset + else: + ship_underwater_st_fieldset = None + + return InputData( + adcp_fieldset=adcp_fieldset, + argo_float_fieldset=argo_float_fieldset, + ctd_fieldset=ctd_fieldset, + drifter_fieldset=drifter_fieldset, + ship_underwater_st_fieldset=ship_underwater_st_fieldset, + ) + + @classmethod + def _load_default_fieldset(cls, directory: str | Path) -> FieldSet: + filenames = { + "U": directory.joinpath("default_uv.nc"), + "V": directory.joinpath("default_uv.nc"), + "S": directory.joinpath("default_s.nc"), + "T": directory.joinpath("default_t.nc"), + } + variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} + dimensions = { + "lon": "longitude", + "lat": "latitude", + "time": "time", + "depth": "depth", + } + + # create the fieldset and set interpolation methods + fieldset = FieldSet.from_netcdf( + filenames, variables, dimensions, allow_time_extrapolation=True + ) + fieldset.T.interp_method = "linear_invdist_land_tracer" + fieldset.S.interp_method = "linear_invdist_land_tracer" + + # make depth negative + for g in fieldset.gridset.grids: + g.depth = -g.depth + + # add bathymetry data + bathymetry_file = directory.joinpath("bathymetry.nc") + bathymetry_variables = ("bathymetry", "deptho") + bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"} + bathymetry_field = Field.from_netcdf( + bathymetry_file, bathymetry_variables, bathymetry_dimensions + ) + # make depth negative + bathymetry_field.data = -bathymetry_field.data + fieldset.add_field(bathymetry_field) + + # read in data already + fieldset.computeTimeChunk(0, 1) + + return fieldset + + @classmethod + def _load_drifter_fieldset(cls, directory: str | Path) -> FieldSet: + filenames = { + "U": directory.joinpath("drifter_uv.nc"), + "V": directory.joinpath("drifter_uv.nc"), + "T": directory.joinpath("drifter_t.nc"), + } + variables = {"U": "uo", "V": "vo", "T": "thetao"} + dimensions = { + "lon": "longitude", + "lat": "latitude", + "time": "time", + "depth": "depth", + } + + fieldset = FieldSet.from_netcdf( + filenames, variables, dimensions, allow_time_extrapolation=False + ) + fieldset.T.interp_method = "linear_invdist_land_tracer" + + # make depth negative + for g in fieldset.gridset.grids: + g.depth = -g.depth + + # read in data already + fieldset.computeTimeChunk(0, 1) + + return fieldset + + @classmethod + def _load_argo_float_fieldset(cls, directory: str | Path) -> FieldSet: + filenames = { + "U": directory.joinpath("argo_float_uv.nc"), + "V": directory.joinpath("argo_float_uv.nc"), + "S": directory.joinpath("argo_float_s.nc"), + "T": directory.joinpath("argo_float_t.nc"), + } + variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} + dimensions = { + "lon": "longitude", + "lat": "latitude", + "time": "time", + "depth": "depth", + } + + fieldset = FieldSet.from_netcdf( + filenames, variables, dimensions, allow_time_extrapolation=False + ) + fieldset.T.interp_method = "linear_invdist_land_tracer" + fieldset.S.interp_method = "linear_invdist_land_tracer" + + # make depth negative + for g in fieldset.gridset.grids: + if max(g.depth) > 0: + g.depth = -g.depth + + # read in data already + fieldset.computeTimeChunk(0, 1) + + return fieldset diff --git a/virtual_ship/expedition/instrument_type.py b/virtual_ship/expedition/instrument_type.py new file mode 100644 index 00000000..556d8464 --- /dev/null +++ b/virtual_ship/expedition/instrument_type.py @@ -0,0 +1,11 @@ +"""InstrumentType Enum.""" + +from enum import Enum + + +class InstrumentType(Enum): + """Types of instruments.""" + + CTD = "CTD" + DRIFTER = "DRIFTER" + ARGO_FLOAT = "ARGO_FLOAT" diff --git a/virtual_ship/expedition/schedule.py b/virtual_ship/expedition/schedule.py new file mode 100644 index 00000000..a92ba55b --- /dev/null +++ b/virtual_ship/expedition/schedule.py @@ -0,0 +1,44 @@ +"""Schedule class.""" + +from __future__ import annotations + +from pathlib import Path + +import pydantic +import yaml + +from .waypoint import Waypoint + + +class Schedule(pydantic.BaseModel): + """Schedule of the virtual ship.""" + + waypoints: list[Waypoint] + + model_config = pydantic.ConfigDict(extra="forbid") + + def to_yaml(self, file_path: str | Path) -> None: + """ + Write schedule to yaml file. + + :param file_path: Path to the file to write to. + """ + with open(file_path, "w") as file: + yaml.dump( + self.model_dump( + by_alias=True, + ), + file, + ) + + @classmethod + def from_yaml(cls, file_path: str | Path) -> Schedule: + """ + Load schedule from yaml file. + + :param file_path: Path to the file to load from. + :returns: The schedule. + """ + with open(file_path, "r") as file: + data = yaml.safe_load(file) + return Schedule(**data) diff --git a/virtual_ship/expedition/ship_config.py b/virtual_ship/expedition/ship_config.py new file mode 100644 index 00000000..8e4564fc --- /dev/null +++ b/virtual_ship/expedition/ship_config.py @@ -0,0 +1,156 @@ +"""ShipConfig and supporting classes.""" + +from __future__ import annotations + +from datetime import timedelta +from pathlib import Path + +import pydantic +import yaml + + +class ArgoFloatConfig(pydantic.BaseModel): + """Configuration for argos floats.""" + + min_depth_meter: float = pydantic.Field(le=0.0) + max_depth_meter: float = pydantic.Field(le=0.0) + drift_depth_meter: float = pydantic.Field(le=0.0) + vertical_speed_meter_per_second: float = pydantic.Field(lt=0.0) + cycle_days: float = pydantic.Field(gt=0.0) + drift_days: float = pydantic.Field(gt=0.0) + + +class ADCPConfig(pydantic.BaseModel): + """Configuration for ADCP instrument.""" + + max_depth_meter: float = pydantic.Field(le=0.0) + num_bins: int = pydantic.Field(gt=0.0) + period: timedelta = pydantic.Field( + serialization_alias="period_minutes", + validation_alias="period_minutes", + gt=timedelta(), + ) + + model_config = pydantic.ConfigDict(populate_by_name=True) + + @pydantic.field_serializer("period") + def _serialize_period(self, value: timedelta, _info): + return value.total_seconds() / 60.0 + + +class CTDConfig(pydantic.BaseModel): + """Configuration for CTD instrument.""" + + stationkeeping_time: timedelta = pydantic.Field( + serialization_alias="stationkeeping_time_minutes", + validation_alias="stationkeeping_time_minutes", + gt=timedelta(), + ) + min_depth_meter: float = pydantic.Field(le=0.0) + max_depth_meter: float = pydantic.Field(le=0.0) + + model_config = pydantic.ConfigDict(populate_by_name=True) + + @pydantic.field_serializer("stationkeeping_time") + def _serialize_stationkeeping_time(self, value: timedelta, _info): + return value.total_seconds() / 60.0 + + +class ShipUnderwaterSTConfig(pydantic.BaseModel): + """Configuration for underwater ST.""" + + period: timedelta = pydantic.Field( + serialization_alias="period_minutes", + validation_alias="period_minutes", + gt=timedelta(), + ) + + model_config = pydantic.ConfigDict(populate_by_name=True) + + @pydantic.field_serializer("period") + def _serialize_period(self, value: timedelta, _info): + return value.total_seconds() / 60.0 + + +class DrifterConfig(pydantic.BaseModel): + """Configuration for drifters.""" + + depth_meter: float = pydantic.Field(le=0.0) + lifetime: timedelta = pydantic.Field( + serialization_alias="lifetime_minutes", + validation_alias="lifetime_minutes", + gt=timedelta(), + ) + + model_config = pydantic.ConfigDict(populate_by_name=True) + + @pydantic.field_serializer("lifetime") + def _serialize_lifetime(self, value: timedelta, _info): + return value.total_seconds() / 60.0 + + +class ShipConfig(pydantic.BaseModel): + """Configuration of the virtual ship.""" + + ship_speed_meter_per_second: float = pydantic.Field(gt=0.0) + """ + Velocity of the ship in meters per second. + """ + + argo_float_config: ArgoFloatConfig | None = None + """ + Argo float configuration. + + If None, no argo floats can be deployed. + """ + + adcp_config: ADCPConfig | None = None + """ + ADCP configuration. + + If None, no ADCP measurements will be performed. + """ + + ctd_config: CTDConfig | None = None + """ + CTD configuration. + + If None, no CTDs can be cast. + """ + + ship_underwater_st_config: ShipUnderwaterSTConfig | None = None + """ + Ship underwater salinity temperature measurementconfiguration. + + If None, no ST measurements will be performed. + """ + + drifter_config: DrifterConfig | None = None + """ + Drifter configuration. + + If None, no drifters can be deployed. + """ + + model_config = pydantic.ConfigDict(extra="forbid") + + def to_yaml(self, file_path: str | Path) -> None: + """ + Write config to yaml file. + + :param file_path: Path to the file to write to. + """ + with open(file_path, "w") as file: + yaml.dump(self.model_dump(by_alias=True), file) + + @classmethod + def from_yaml(cls, file_path: str | Path) -> ShipConfig: + """ + Load config from yaml file. + + :param file_path: Path to the file to load from. + :returns: The config. + """ + with open(file_path, "r") as file: + data = yaml.safe_load(file) + return ShipConfig(**data) diff --git a/virtual_ship/expedition/simulate_measurements.py b/virtual_ship/expedition/simulate_measurements.py new file mode 100644 index 00000000..7ff7868f --- /dev/null +++ b/virtual_ship/expedition/simulate_measurements.py @@ -0,0 +1,104 @@ +"""simulate_measurements function.""" + +from datetime import timedelta +from pathlib import Path + +from ..instruments.adcp import simulate_adcp +from ..instruments.argo_float import simulate_argo_floats +from ..instruments.ctd import simulate_ctd +from ..instruments.drifter import simulate_drifters +from ..instruments.ship_underwater_st import simulate_ship_underwater_st +from .input_data import InputData +from .ship_config import ShipConfig +from .simulate_schedule import MeasurementsToSimulate + + +def simulate_measurements( + expedition_dir: str | Path, + ship_config: ShipConfig, + input_data: InputData, + measurements: MeasurementsToSimulate, +) -> None: + """ + Simulate measurements using Parcels. + + Saves everything in expedition_dir/results. + + :param expedition_dir: Base directory of the expedition. + :param ship_config: Ship configuration. + :param input_data: Input data for simulation. + :param measurements: The measurements to simulate. + :raises RuntimeError: In case fieldsets of configuration is not provided. Make sure to check this before calling this function. + """ + if isinstance(expedition_dir, str): + expedition_dir = Path(expedition_dir) + + if len(measurements.ship_underwater_sts) > 0: + print("Simulating onboard salinity and temperature measurements.") + if ship_config.ship_underwater_st_config is None: + raise RuntimeError("No configuration for ship underwater ST provided.") + if input_data.ship_underwater_st_fieldset is None: + raise RuntimeError("No fieldset for ship underwater ST provided.") + simulate_ship_underwater_st( + fieldset=input_data.ship_underwater_st_fieldset, + out_path=expedition_dir.joinpath("results", "ship_underwater_st.zarr"), + depth=-2, + sample_points=measurements.ship_underwater_sts, + ) + + if len(measurements.adcps) > 0: + print("Simulating onboard ADCP.") + if ship_config.adcp_config is None: + raise RuntimeError("No configuration for ADCP provided.") + if input_data.adcp_fieldset is None: + raise RuntimeError("No fieldset for ADCP provided.") + simulate_adcp( + fieldset=input_data.adcp_fieldset, + out_path=expedition_dir.joinpath("results", "adcp.zarr"), + max_depth=ship_config.adcp_config.max_depth_meter, + min_depth=-5, + num_bins=ship_config.adcp_config.num_bins, + sample_points=measurements.adcps, + ) + + if len(measurements.ctds) > 0: + print("Simulating CTD casts.") + if ship_config.ctd_config is None: + raise RuntimeError("No configuration for CTD provided.") + if input_data.ctd_fieldset is None: + raise RuntimeError("No fieldset for CTD provided.") + simulate_ctd( + out_path=expedition_dir.joinpath("results", "ctd.zarr"), + fieldset=input_data.ctd_fieldset, + ctds=measurements.ctds, + outputdt=timedelta(seconds=10), + ) + + if len(measurements.drifters) > 0: + print("Simulating drifters") + if ship_config.drifter_config is None: + raise RuntimeError("No configuration for drifters provided.") + if input_data.drifter_fieldset is None: + raise RuntimeError("No fieldset for drifters provided.") + simulate_drifters( + out_path=expedition_dir.joinpath("results", "drifters.zarr"), + fieldset=input_data.drifter_fieldset, + drifters=measurements.drifters, + outputdt=timedelta(hours=5), + dt=timedelta(minutes=5), + endtime=None, + ) + + if len(measurements.argo_floats) > 0: + print("Simulating argo floats") + if ship_config.argo_float_config is None: + raise RuntimeError("No configuration for argo floats provided.") + if input_data.argo_float_fieldset is None: + raise RuntimeError("No fieldset for argo floats provided.") + simulate_argo_floats( + out_path=expedition_dir.joinpath("results", "argo_floats.zarr"), + argo_floats=measurements.argo_floats, + fieldset=input_data.argo_float_fieldset, + outputdt=timedelta(minutes=5), + endtime=None, + ) diff --git a/virtual_ship/expedition/simulate_schedule.py b/virtual_ship/expedition/simulate_schedule.py new file mode 100644 index 00000000..85df8696 --- /dev/null +++ b/virtual_ship/expedition/simulate_schedule.py @@ -0,0 +1,264 @@ +"""simulate_schedule function and supporting classes.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from datetime import datetime, timedelta + +import pyproj + +from ..instruments.argo_float import ArgoFloat +from ..instruments.ctd import CTD +from ..instruments.drifter import Drifter +from ..location import Location +from ..spacetime import Spacetime +from .instrument_type import InstrumentType +from .schedule import Schedule +from .ship_config import ShipConfig +from .waypoint import Waypoint + + +@dataclass +class ScheduleOk: + """Result of schedule that could be completed.""" + + time: datetime + measurements_to_simulate: MeasurementsToSimulate + + +@dataclass +class ScheduleProblem: + """Result of schedule that could not be fully completed.""" + + time: datetime + failed_waypoint_i: int + + +@dataclass +class MeasurementsToSimulate: + """The measurements to simulate, as concluded from schedule simulation.""" + + 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) + + +def simulate_schedule( + projection: pyproj.Geod, ship_config: ShipConfig, schedule: Schedule +) -> ScheduleOk | ScheduleProblem: + """ + Simulate a schedule. + + :param projection: The projection to use for sailing. + :param ship_config: Ship configuration. + :param schedule: The schedule to simulate. + :returns: Either the results of a successfully simulated schedule, or information on where the schedule became infeasible. + """ + return _ScheduleSimulator(projection, ship_config, schedule).simulate() + + +class _ScheduleSimulator: + _projection: pyproj.Geod + _ship_config: ShipConfig + _schedule: Schedule + + _time: datetime + """Current time.""" + _location: Location + """Current ship location.""" + + _measurements_to_simulate: MeasurementsToSimulate + + _next_adcp_time: datetime + """Next moment ADCP measurement will be done.""" + _next_ship_underwater_st_time: datetime + """Next moment ship underwater ST measurement will be done.""" + + def __init__( + self, projection: pyproj.Geod, ship_config: ShipConfig, schedule: Schedule + ) -> None: + self._projection = projection + self._ship_config = ship_config + self._schedule = schedule + + assert ( + self._schedule.waypoints[0].time is not None + ), "First waypoint must have a time. This should have been verified before calling this function." + self._time = schedule.waypoints[0].time + self._location = schedule.waypoints[0].location + + self._measurements_to_simulate = MeasurementsToSimulate() + + self._next_adcp_time = self._time + self._next_ship_underwater_st_time = self._time + + def simulate(self) -> ScheduleOk | ScheduleProblem: + for wp_i, waypoint in enumerate(self._schedule.waypoints): + # sail towards waypoint + self._progress_time_traveling_towards(waypoint.location) + + # check if waypoint was reached in time + if waypoint.time is not None and self._time > waypoint.time: + print( + f"Waypoint {wp_i} could not be reached in time. Current time: {self._time}. Waypoint time: {waypoint.time}." + ) + return ScheduleProblem(self._time, wp_i) + + # note measurements made at waypoint + time_passed = self._make_measurements(waypoint) + + # 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: + geodinv: tuple[float, float, float] = self._projection.inv( + lons1=self._location.lon, + lats1=self._location.lat, + lons2=location.lon, + lats2=location.lat, + ) + azimuth1 = geodinv[0] + distance_to_next_waypoint = geodinv[2] + time_to_reach = timedelta( + seconds=distance_to_next_waypoint + / self._ship_config.ship_speed_meter_per_second + ) + end_time = self._time + time_to_reach + + # note all ADCP measurements + if self._ship_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 = ( + self._ship_config.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.adcps.append( + Spacetime(location=location, time=time) + ) + + self._next_adcp_time = ( + self._next_adcp_time + self._ship_config.adcp_config.period + ) + + # note all ship underwater ST measurements + if self._ship_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 = ( + self._ship_config.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( + Spacetime(location=location, time=time) + ) + + self._next_ship_underwater_st_time = ( + self._next_ship_underwater_st_time + + self._ship_config.ship_underwater_st_config.period + ) + + self._time = end_time + self._location = location + + def _progress_time_stationary(self, time_passed: timedelta) -> None: + end_time = self._time + time_passed + + # note all ADCP measurements + if self._ship_config.adcp_config is not None: + while self._next_adcp_time <= end_time: + self._measurements_to_simulate.adcps.append( + Spacetime(self._location, self._next_adcp_time) + ) + self._next_adcp_time = ( + self._next_adcp_time + self._ship_config.adcp_config.period + ) + + # note all ship underwater ST measurements + if self._ship_config.ship_underwater_st_config is not None: + while self._next_ship_underwater_st_time <= end_time: + 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._ship_config.ship_underwater_st_config.period + ) + + self._time = end_time + + def _make_measurements(self, waypoint: Waypoint) -> timedelta: + # if there are no instruments, there is no time cost + if waypoint.instrument is None: + return timedelta() + + # make instruments a list even if it's only a single one + instruments = ( + waypoint.instrument + if isinstance(waypoint.instrument, list) + else [waypoint.instrument] + ) + + # time costs of each measurement + time_costs = [timedelta()] + + for instrument in instruments: + if instrument is InstrumentType.ARGO_FLOAT: + self._measurements_to_simulate.argo_floats.append( + ArgoFloat( + spacetime=Spacetime(self._location, self._time), + min_depth=self._ship_config.argo_float_config.min_depth_meter, + max_depth=self._ship_config.argo_float_config.max_depth_meter, + drift_depth=self._ship_config.argo_float_config.drift_depth_meter, + vertical_speed=self._ship_config.argo_float_config.vertical_speed_meter_per_second, + cycle_days=self._ship_config.argo_float_config.cycle_days, + drift_days=self._ship_config.argo_float_config.drift_days, + ) + ) + elif instrument is InstrumentType.CTD: + self._measurements_to_simulate.ctds.append( + CTD( + spacetime=Spacetime(self._location, self._time), + min_depth=self._ship_config.ctd_config.min_depth_meter, + max_depth=self._ship_config.ctd_config.max_depth_meter, + ) + ) + time_costs.append(timedelta(minutes=20)) + elif instrument is InstrumentType.DRIFTER: + self._measurements_to_simulate.drifters.append( + Drifter( + spacetime=Spacetime(self._location, self._time), + depth=self._ship_config.drifter_config.depth_meter, + lifetime=self._ship_config.drifter_config.lifetime, + ) + ) + else: + raise NotImplementedError("Instrument type not supported.") + + # measurements are done in parallel, so return time of longest one + return max(time_costs) diff --git a/virtual_ship/expedition/verify_schedule.py b/virtual_ship/expedition/verify_schedule.py new file mode 100644 index 00000000..55c34587 --- /dev/null +++ b/virtual_ship/expedition/verify_schedule.py @@ -0,0 +1,166 @@ +"""verify_schedule function and supporting classes.""" + +from datetime import timedelta + +import pyproj +from parcels import FieldSet + +from .input_data import InputData +from .instrument_type import InstrumentType +from .schedule import Schedule +from .ship_config import ShipConfig +from .waypoint import Waypoint + + +def verify_schedule( + projection: pyproj.Geod, + ship_config: ShipConfig, + schedule: Schedule, + input_data: InputData, +) -> 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 projection: projection used to sail between waypoints. + :param ship_config: The cruise ship_configuration. + :param schedule: The schedule to verify. + :param input_data: Fieldsets that can be used to check for zero UV condition (is waypoint on land). + :raises PlanningError: If waypoints are not feasible or incorrect. + :raises NotImplementedError: If instrument is in schedule that is not implememented. + """ + if len(schedule.waypoints) == 0: + raise PlanningError("At least one waypoint must be provided.") + + # check first waypoint has a time + if schedule.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 schedule.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 if all waypoints are in water + # this is done by picking an arbitrary provided fieldset and checking if UV is not zero + + print("Verifying all waypoints are on water..") + + # get all available fieldsets + available_fieldsets = [ + fs + for fs in [ + input_data.adcp_fieldset, + input_data.argo_float_fieldset, + input_data.ctd_fieldset, + input_data.drifter_fieldset, + input_data.ship_underwater_st_fieldset, + ] + if fs is not None + ] + + # check if there are any fieldsets, else its an error + if len(available_fieldsets) == 0: + print( + "Cannot verify because no fieldsets have been loaded. This is probably because you are not using any instruments in your schedule. This is not a problem, but carefully check your waypoint locations manually.." + ) + + else: + # pick any + fieldset = available_fieldsets[0] + # get waypoints with 0 UV + land_waypoints = [ + (wp_i, wp) + for wp_i, wp in enumerate(schedule.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]}" + ) + print("Good, all waypoints are on water.") + + # check that ship will arrive on time at each waypoint (in case no unexpected event happen) + time = schedule.waypoints[0].time + for wp_i, (wp, wp_next) in enumerate( + zip(schedule.waypoints, schedule.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_config.ship_speed_meter_per_second + ) + 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 + 1}. location: {wp_next.location} time: {wp_next.time} instrument: {wp_next.instrument}" + ) + else: + time = wp_next.time + + # verify instruments in schedule have configuration + for wp in schedule.waypoints: + if wp.instrument is not None: + for instrument in ( + wp.instrument if isinstance(wp.instrument, list) else [wp.instrument] + ): + if instrument not in InstrumentType: + raise NotImplementedError("Instrument not supported.") + + if ( + instrument == InstrumentType.ARGO_FLOAT + and ship_config.argo_float_config is None + ): + raise PlanningError( + "Planning has waypoint with Argo float instrument, but configuration does not configure Argo floats." + ) + if instrument == InstrumentType.CTD and ship_config.ctd_config is None: + raise PlanningError( + "Planning has waypoint with CTD instrument, but configuration does not configure CTDs." + ) + if ( + instrument == InstrumentType.DRIFTER + and ship_config.drifter_config is None + ): + raise PlanningError( + "Planning has waypoint with drifter instrument, but configuration does not configure drifters." + ) + + +class PlanningError(RuntimeError): + """An error in the schedule.""" + + pass + + +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, + fieldset.gridset.grids[0].depth[0], + waypoint.location.lat, + waypoint.location.lon, + applyConversion=False, + ) == (0.0, 0.0) diff --git a/virtual_ship/waypoint.py b/virtual_ship/expedition/waypoint.py similarity index 75% rename from virtual_ship/waypoint.py rename to virtual_ship/expedition/waypoint.py index bf3d33b4..85e99181 100644 --- a/virtual_ship/waypoint.py +++ b/virtual_ship/expedition/waypoint.py @@ -3,8 +3,8 @@ from dataclasses import dataclass from datetime import datetime +from ..location import Location from .instrument_type import InstrumentType -from .location import Location @dataclass @@ -13,4 +13,4 @@ class Waypoint: location: Location time: datetime | None = None - instrument: InstrumentType | None = None + instrument: InstrumentType | list[InstrumentType] | None = None diff --git a/virtual_ship/instrument_type.py b/virtual_ship/instrument_type.py deleted file mode 100644 index 9a19b814..00000000 --- a/virtual_ship/instrument_type.py +++ /dev/null @@ -1,11 +0,0 @@ -"""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/ctd.py b/virtual_ship/instruments/ctd.py index 50a397af..0d017fa7 100644 --- a/virtual_ship/instruments/ctd.py +++ b/virtual_ship/instruments/ctd.py @@ -82,7 +82,9 @@ def simulate_ctd( fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) # deploy time for all ctds should be later than fieldset start time - if not all([ctd.spacetime.time >= fieldset_starttime for ctd in ctds]): + if not all( + [np.datetime64(ctd.spacetime.time) >= fieldset_starttime for ctd in ctds] + ): raise ValueError("CTD deployed before fieldset starts.") # depth the ctd will go to. shallowest between ctd max depth and bathymetry. diff --git a/virtual_ship/instruments/ship_underwater_st.py b/virtual_ship/instruments/ship_underwater_st.py index 3b4ac59e..3d2050d8 100644 --- a/virtual_ship/instruments/ship_underwater_st.py +++ b/virtual_ship/instruments/ship_underwater_st.py @@ -10,24 +10,20 @@ # there is some overhead with JITParticle and this ends up being significantly faster _ShipSTParticle = ScipyParticle.add_variables( [ - Variable("salinity", dtype=np.float32, initial=np.nan), - Variable("temperature", dtype=np.float32, initial=np.nan), + Variable("S", dtype=np.float32, initial=np.nan), + Variable("T", dtype=np.float32, initial=np.nan), ] ) # define function sampling Salinity def _sample_salinity(particle, fieldset, time): - particle.salinity = fieldset.salinity[ - time, particle.depth, particle.lat, particle.lon - ] + particle.S = fieldset.S[time, particle.depth, particle.lat, particle.lon] # define function sampling Temperature def _sample_temperature(particle, fieldset, time): - particle.temperature = fieldset.temperature[ - time, particle.depth, particle.lat, particle.lon - ] + particle.T = fieldset.T[time, particle.depth, particle.lat, particle.lon] def simulate_ship_underwater_st( diff --git a/virtual_ship/planning_error.py b/virtual_ship/planning_error.py deleted file mode 100644 index 5d599dd2..00000000 --- a/virtual_ship/planning_error.py +++ /dev/null @@ -1,7 +0,0 @@ -"""PlanningError Exception.""" - - -class PlanningError(RuntimeError): - """An error when checking the schedule or during sailing.""" - - pass diff --git a/virtual_ship/postprocess.py b/virtual_ship/postprocess.py deleted file mode 100644 index a234dab9..00000000 --- a/virtual_ship/postprocess.py +++ /dev/null @@ -1,71 +0,0 @@ -"""postprocess function.""" - -import os -import shutil - -import numpy as np -import xarray as xr -from scipy.ndimage import uniform_filter1d - - -def postprocess(): - """Postprocesses CTD data and writes to csv files.""" - if os.path.isdir(os.path.join("results", "CTDs")): - i = 0 - filenames = os.listdir(os.path.join("results", "CTDs")) - for filename in sorted(filenames): - if filename.endswith(".zarr"): - try: # too many errors, just skip the faulty zarr files - i += 1 - # Open output and read to x, y, z - ds = xr.open_zarr(os.path.join("results", "CTDs", filename)) - x = ds["lon"][:].squeeze() - y = ds["lat"][:].squeeze() - z = ds["z"][:].squeeze() - time = ds["time"][:].squeeze() - T = ds["temperature"][:].squeeze() - S = ds["salinity"][:].squeeze() - ds.close() - - random_walk = np.random.random() / 10 - z_norm = (z - np.min(z)) / (np.max(z) - np.min(z)) - t_norm = np.linspace(0, 1, num=len(time)) - # add smoothed random noise scaled with depth - # and random (reversed) diversion from initial through time scaled with depth - S = S + uniform_filter1d( - np.random.random(S.shape) / 5 * (1 - z_norm) - + random_walk - * (np.max(S).values - np.min(S).values) - * (1 - z_norm) - * t_norm - / 10, - max(int(len(time) / 40), 1), - ) - T = T + uniform_filter1d( - np.random.random(T.shape) * 5 * (1 - z_norm) - - random_walk - / 2 - * (np.max(T).values - np.min(T).values) - * (1 - z_norm) - * t_norm - / 10, - max(int(len(time) / 20), 1), - ) - - # reshaping data to export to csv - header = "pressure [dbar],temperature [degC],salinity [g kg-1]" - data = np.column_stack([-z, T, S]) - new_line = "\n" - np.savetxt( - f"{os.path.join('results', 'CTDs', 'CTD_station_')}{i}.csv", - data, - fmt="%.4f", - header=header, - delimiter=",", - comments=f'longitude,{x[0].values},"{x.attrs}"{new_line}latitude,{y[0].values},"{y.attrs}"{new_line}start time,{time[0].values}{new_line}end time,{time[-1].values}{new_line}', - ) - shutil.rmtree(os.path.join("results", "CTDs", filename)) - except TypeError: - print(f"CTD file {filename} seems faulty, skipping.") - continue - print("CTD data postprocessed.") diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py deleted file mode 100644 index 5da43dae..00000000 --- a/virtual_ship/sailship.py +++ /dev/null @@ -1,495 +0,0 @@ -"""sailship function.""" - -from __future__ import annotations - -import os -from contextlib import contextmanager -from dataclasses import dataclass, field -from datetime import datetime, timedelta -from typing import Generator - -import pyproj -from parcels import FieldSet -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 .planning_error import PlanningError -from .spacetime import Spacetime -from .virtual_ship_config import VirtualShipConfig -from .waypoint import Waypoint - - -def sailship(config: VirtualShipConfig): - """ - Use Parcels to simulate a virtual ship expedition. - - :param config: The expedition configuration. - """ - config.verify() - - # projection used to sail between waypoints - projection = pyproj.Geod(ellps="WGS84") - - _verify_waypoints(projection=projection, config=config) - - # simulate the sailing and aggregate what measurements should be simulated - schedule_results = _simulate_schedule( - projection=projection, - config=config, - ) - - # simulate the measurements - - 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 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 CTD casts.") - simulate_ctd( - out_path=os.path.join("results", "ctd.zarr"), - 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_config.fieldset, - drifters=schedule_results.measurements_to_simulate.drifters, - outputdt=timedelta(hours=5), - dt=timedelta(minutes=5), - endtime=None, - ) - - print("Simulating argo floats") - simulate_argo_floats( - out_path=os.path.join("results", "argo_floats.zarr"), - argo_floats=schedule_results.measurements_to_simulate.argo_floats, - fieldset=config.argo_float_config.fieldset, - outputdt=timedelta(minutes=5), - endtime=None, - ) - - # calculate cruise cost - assert ( - config.schedule.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.schedule.waypoints[0].time - cost = costs(config, time_past) - print(f"This cruise took {time_past} and would have cost {cost:,.0f} euros.") - - -def _simulate_schedule( - projection: pyproj.Geod, - config: VirtualShipConfig, -) -> _ScheduleResults: - """ - Simulate the sailing and aggregate the virtual measurements that should be taken. - - :param projection: Projection used to sail between waypoints. - :param config: The cruise configuration. - :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 = _Cruise( - Spacetime( - config.schedule.waypoints[0].location, config.schedule.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, - ) - ) - - # sail to each waypoint while executing tasks - for waypoint in config.schedule.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, - ) - ) - - -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, - ) - ) - - -def _verify_waypoints( - 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 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.schedule.waypoints) == 0: - raise PlanningError("At least one waypoint must be provided.") - - # check first waypoint has a time - if config.schedule.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 config.schedule.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 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.schedule.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 = config.schedule.waypoints[0].time - for wp_i, (wp, wp_next) in enumerate( - zip(config.schedule.waypoints, config.schedule.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 / config.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 - - -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) diff --git a/virtual_ship/schedule.py b/virtual_ship/schedule.py deleted file mode 100644 index 5418d2e4..00000000 --- a/virtual_ship/schedule.py +++ /dev/null @@ -1,60 +0,0 @@ -"""Schedule class.""" - -from __future__ import annotations - -from dataclasses import dataclass - -import yaml - -from .location import Location -from .waypoint import Waypoint - - -@dataclass -class Schedule: - """Schedule of the virtual ship.""" - - waypoints: list[Waypoint] - - @classmethod - def from_yaml(cls, path: str) -> Schedule: - """ - Load schedule from YAML file. - - :param path: The file to read from. - :returns: Schedule of waypoints from the YAML file. - """ - with open(path, "r") as in_file: - data = yaml.safe_load(in_file) - waypoints = [ - Waypoint( - location=Location(waypoint["lat"], waypoint["lon"]), - time=waypoint["time"], - instrument=waypoint["instrument"], - ) - for waypoint in data - ] - return Schedule(waypoints) - - def to_yaml(self, path: str) -> None: - """ - Save schedule to YAML file. - - :param path: The file to write to. - """ - with open(path, "w") as out_file: - print( - yaml.dump( - [ - { - "lat": waypoint.location.lat, - "lon": waypoint.location.lon, - "time": waypoint.time, - "instrument": waypoint.instrument, - } - for waypoint in self.waypoints - ], - out_file, - ) - ) - pass diff --git a/virtual_ship/virtual_ship_config.py b/virtual_ship/virtual_ship_config.py deleted file mode 100644 index 9dd20a20..00000000 --- a/virtual_ship/virtual_ship_config.py +++ /dev/null @@ -1,102 +0,0 @@ -"""VirtualShipConfig class.""" - -from dataclasses import dataclass -from datetime import timedelta - -from parcels import FieldSet - -from .schedule import Schedule - - -@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 - - schedule: Schedule - - 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.schedule.waypoints) < 2: - raise ValueError("Waypoints require at least a start and an end.") - - 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.")