diff --git a/.github/workflows/codetools.yml b/.github/workflows/codetools.yml index b53a9d5b..7027fec2 100644 --- a/.github/workflows/codetools.yml +++ b/.github/workflows/codetools.yml @@ -4,6 +4,7 @@ on: [push, pull_request] env: PACKAGE: virtual_ship + TESTS: tests jobs: codetools: @@ -19,7 +20,7 @@ jobs: - name: install run: pip install ".[dev]" - name: flake8 - run: flake8 ./$PACKAGE + run: flake8 ./$PACKAGE ./$TESTS - name: pydocstyle run: pydocstyle ./$PACKAGE - name: sort-all @@ -28,6 +29,18 @@ jobs: [[ -z $(git status -s) ]] git checkout -- . - name: black - run: black --diff --check ./$PACKAGE + run: black --diff --check ./$PACKAGE ./$TESTS - name: isort - run: isort --check-only --diff ./$PACKAGE + run: isort --check-only --diff ./$PACKAGE ./$TESTS + + tests: + runs-on: ubuntu-20.04 + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5.1.0 + with: + python-version: 3.12 + - name: install + run: pip install ".[dev]" + - name: run_tests + run: pytest --cov=virtual_ship tests diff --git a/.gitignore b/.gitignore index 160e42a9..924a28fe 100644 --- a/.gitignore +++ b/.gitignore @@ -172,3 +172,5 @@ cython_debug/ # Auto generated by setuptools scm virtual_ship/_version_setup.py + +.vscode/ diff --git a/codetools.sh b/codetools.sh old mode 100644 new mode 100755 index 960a0eff..a359c5d6 --- a/codetools.sh +++ b/codetools.sh @@ -9,11 +9,12 @@ set -e cd "$(dirname "$0")" PACKAGE=virtual_ship +TESTS=tests echo "--------------" echo "flake8" echo "--------------" -flake8 ./$PACKAGE +flake8 ./$PACKAGE ./$TESTS # darglint is ran as a plugin for flake8. echo "--------------" @@ -29,10 +30,10 @@ find ./$PACKAGE -type f -name '__init__.py' -print0 | xargs -0 sort-all echo "--------------" echo "black" echo "--------------" -black ./$PACKAGE +black ./$PACKAGE ./$TESTS echo "--------------" echo "isort" echo "--------------" -isort ./$PACKAGE +isort ./$PACKAGE ./$TESTS diff --git a/pyproject.toml b/pyproject.toml index 69ce0bf2..033d157a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,16 @@ dependencies = [ "parcels >= 3, < 4", "scipy >= 1, < 2", "xarray >= 2023, < 2024", + # The packages below are parcels packages not included in their distribution requirements. + "numpy >= 1, < 2", + "cgen >= 2020, < 2021", + "dask >= 2023, < 2025", + "cftime >= 1, < 2", + "psutil >= 1, < 2", + "netCDF4 >= 1, < 2", + "zarr >= 2, < 3", + "tqdm >= 4, < 5", + "pymbolic >= 2022, < 2023", ] [project.urls] @@ -53,6 +63,9 @@ dev = [ "isort == 5.13.2", "pydocstyle == 6.3.0", "sort-all == 1.2.0", + "pytest == 8.2.0", + "pytest-cov == 5.0.0", + "codecov == 2.1.13", ] [tool.isort] diff --git a/tests.sh b/tests.sh new file mode 100755 index 00000000..a6203336 --- /dev/null +++ b/tests.sh @@ -0,0 +1,5 @@ +#!/bin/sh + +# Runs the tests and creates a code coverage report. + +pytest --cov=virtual_ship --cov-report=html tests diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 00000000..01e91653 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,7 @@ +import pytest + + +# Set the working directory for each test to the directory of that test. +@pytest.fixture(autouse=True) +def change_test_dir(request, monkeypatch): + monkeypatch.chdir(request.fspath.dirname) diff --git a/tests/instruments/test_argo.py b/tests/instruments/test_argo.py new file mode 100644 index 00000000..349ed051 --- /dev/null +++ b/tests/instruments/test_argo.py @@ -0,0 +1,48 @@ +"""Test the simulation of Argo floats.""" + +from datetime import timedelta + +import numpy as np +from parcels import FieldSet + +from virtual_ship.instruments import Location +from virtual_ship.instruments.argo_float import ArgoFloat, simulate_argo_floats + + +def test_simulate_argo_floats() -> None: + DRIFT_DEPTH = -1000 + MAX_DEPTH = -2000 + VERTICLE_SPEED = -0.10 + CYCLE_DAYS = 10 + DRIFT_DAYS = 9 + + fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "T": 0, "S": 0}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + min_depth = -fieldset.U.depth[0] + + argo_floats = [ + ArgoFloat( + location=Location(latitude=0, longitude=0), + deployment_time=0, + min_depth=min_depth, + max_depth=MAX_DEPTH, + drift_depth=DRIFT_DEPTH, + vertical_speed=VERTICLE_SPEED, + cycle_days=CYCLE_DAYS, + drift_days=DRIFT_DAYS, + ) + ] + + simulate_argo_floats( + argo_floats=argo_floats, + fieldset=fieldset, + out_file_name="test", + outputdt=timedelta(minutes=5), + ) diff --git a/tests/sailship_config.json b/tests/sailship_config.json new file mode 100644 index 00000000..23cf4c1a --- /dev/null +++ b/tests/sailship_config.json @@ -0,0 +1,40 @@ +{ + "region_of_interest": { + "North": 64, + "East": -23, + "South": 59, + "West": -43 + }, + "requested_ship_time": { + "start": "2022-01-01T00:00:00", + "end": "2022-01-2T00:00:00" + }, + "route_coordinates": [ + [-23.071289, 63.743631], [-23.081289, 63.743631], [-23.091289, 63.743631] + ], + "underway_data": true, + "ADCP_data": false, + "ADCP_settings": { + "max_depth": -1000, + "bin_size_m": 24 + }, + "CTD_locations": [ + [-23.071289, 63.743631] + ], + "CTD_settings": { + "max_depth": "max" + }, + "drifter_deploylocations": [ + + ], + "argo_deploylocations": [ + [-23.081289, 63.743631] + ], + "argo_characteristics": { + "driftdepth": -1000, + "maxdepth": -2000, + "vertical_speed": -0.10, + "cycle_days" : 10, + "drift_days": 9 + } + } diff --git a/tests/test_sailship.py b/tests/test_sailship.py new file mode 100644 index 00000000..afd8a95f --- /dev/null +++ b/tests/test_sailship.py @@ -0,0 +1,43 @@ +"""Performs a complete cruise with virtual ship.""" + +import numpy as np +from parcels import FieldSet + +from virtual_ship.sailship import sailship +from virtual_ship.virtual_ship_configuration import VirtualShipConfiguration + + +def test_sailship() -> None: + ctd_fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "S": 0, "T": 0, "bathymetry": 0}, + {"lon": 0, "lat": 0}, + ) + ctd_fieldset.add_constant("maxtime", ctd_fieldset.U.grid.time_full[-1]) + ctd_fieldset.add_constant("mindepth", -ctd_fieldset.U.depth[0]) + ctd_fieldset.add_constant("max_depth", -ctd_fieldset.U.depth[-1]) + + drifter_fieldset = FieldSet.from_data( + { + "U": 0, + "V": 0, + }, + {"lon": 0, "lat": 0}, + ) + + argo_float_fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "T": 0, "S": 0}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + config = VirtualShipConfiguration( + "sailship_config.json", + ctd_fieldset=ctd_fieldset, + drifter_fieldset=drifter_fieldset, + argo_float_fieldset=argo_float_fieldset, + ) + + sailship(config) diff --git a/virtual_ship/__init__.py b/virtual_ship/__init__.py index d885347c..bf152389 100644 --- a/virtual_ship/__init__.py +++ b/virtual_ship/__init__.py @@ -1 +1,5 @@ """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 + +__all__ = ["instruments", "sailship"] diff --git a/virtual_ship/argo_deployments.py b/virtual_ship/argo_deployments.py deleted file mode 100644 index 5a31870f..00000000 --- a/virtual_ship/argo_deployments.py +++ /dev/null @@ -1,189 +0,0 @@ -"""argo_deployments function.""" - -import datetime -import math -import os -from datetime import timedelta - -import numpy as np -from parcels import ( - AdvectionRK4, - FieldSet, - JITParticle, - ParticleSet, - StatusCode, - Variable, -) - - -def argo_deployments(config, argo_time): - """ - Deploy argo floats. - - :param config: The cruise configuration. - :param argo_time: TODO - """ - # particle_* such as particle_ddepth are local variables defined by parcels. - # See https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html#Background - - if len(config.argo_deploylocations) > 0: - - fieldset = create_argo_fieldset(config) - - # Define the new Kernel that mimics Argo vertical movement - def ArgoVerticalMovement(particle, fieldset, time): - - if particle.cycle_phase == 0: - # Phase 0: Sinking with vertical_speed until depth is driftdepth - particle_ddepth += ( # noqa See comment above about particle_* variables. - fieldset.vertical_speed * particle.dt - ) - if particle.depth + particle_ddepth <= fieldset.driftdepth: - particle_ddepth = fieldset.driftdepth - particle.depth - particle.cycle_phase = 1 - - elif particle.cycle_phase == 1: - # Phase 1: Drifting at depth for drifttime seconds - particle.drift_age += particle.dt - if particle.drift_age >= fieldset.drift_days * 86400: - particle.drift_age = 0 # reset drift_age for next cycle - particle.cycle_phase = 2 - - elif particle.cycle_phase == 2: - # Phase 2: Sinking further to maxdepth - particle_ddepth += fieldset.vertical_speed * particle.dt - if particle.depth + particle_ddepth <= fieldset.maxdepth: - particle_ddepth = fieldset.maxdepth - particle.depth - particle.cycle_phase = 3 - - elif particle.cycle_phase == 3: - # Phase 3: Rising with vertical_speed until at surface - particle_ddepth -= fieldset.vertical_speed * particle.dt - particle.cycle_age += ( - particle.dt - ) # solve issue of not updating cycle_age during ascent - if particle.depth + particle_ddepth >= fieldset.mindepth: - particle_ddepth = fieldset.mindepth - particle.depth - particle.temperature = ( - math.nan - ) # reset temperature to NaN at end of sampling cycle - particle.salinity = math.nan # idem - particle.cycle_phase = 4 - else: - particle.temperature = fieldset.T[ - time, particle.depth, particle.lat, particle.lon - ] - particle.salinity = fieldset.S[ - time, particle.depth, particle.lat, particle.lon - ] - - elif particle.cycle_phase == 4: - # Phase 4: Transmitting at surface until cycletime is reached - if particle.cycle_age > fieldset.cycle_days * 86400: - particle.cycle_phase = 0 - particle.cycle_age = 0 - - if particle.state == StatusCode.Evaluate: - particle.cycle_age += particle.dt # update cycle_age - - def KeepAtSurface(particle, fieldset, time): - # Prevent error when float reaches surface - if particle.state == StatusCode.ErrorThroughSurface: - particle.depth = fieldset.mindepth - particle.state = StatusCode.Success - - def CheckError(particle, fieldset, time): - if particle.state >= 50: # This captures all Errors - particle.delete() - - class ArgoParticle(JITParticle): - cycle_phase = Variable("cycle_phase", dtype=np.int32, initial=0.0) - cycle_age = Variable("cycle_age", dtype=np.float32, initial=0.0) - drift_age = Variable("drift_age", dtype=np.float32, initial=0.0) - salinity = Variable("salinity", initial=np.nan) - temperature = Variable("temperature", initial=np.nan) - - # initialize argo floats - lon = [] - lat = [] - for i in range(len(config.argo_deploylocations)): - lon.append(config.argo_deploylocations[i][0]) - lat.append(config.argo_deploylocations[i][1]) - time = argo_time - - # Create and execute argo particles - argoset = ParticleSet( - fieldset=fieldset, - pclass=ArgoParticle, - lon=lon, - lat=lat, - depth=np.repeat(fieldset.mindepth, len(time)), - time=time, - ) - argo_output_file = argoset.ParticleFile( - name=os.path.join("results", "Argos.zarr"), - outputdt=timedelta(minutes=5), - chunks=(1, 500), - ) - fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) - argo_endtime = np.array( - ( - datetime.datetime.strptime( - config.requested_ship_time["start"], "%Y-%m-%dT%H:%M:%S" - ) - + timedelta(weeks=6) - ) - ).astype("datetime64[ms]") - - argoset.execute( - [ - ArgoVerticalMovement, - AdvectionRK4, - KeepAtSurface, - CheckError, - ], # list of kernels to be executed - endtime=min(fieldset_endtime, argo_endtime), - dt=timedelta(minutes=5), - output_file=argo_output_file, - ) - - -def create_argo_fieldset(config): - """ - Create a fieldset from netcdf files for argo floats, return fieldset with negative depth values. - - :param config: The cruise configuration. - :returns: The fieldset. - """ - datadirname = os.path.dirname(__file__) - filenames = { - "U": os.path.join(datadirname, "argodata_UV.nc"), - "V": os.path.join(datadirname, "argodata_UV.nc"), - "S": os.path.join(datadirname, "argodata_S.nc"), - "T": os.path.join(datadirname, "argodata_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=False - ) - fieldset.T.interp_method = "linear_invdist_land_tracer" - for g in fieldset.gridset.grids: - if max(g.depth) > 0: - g.depth = -g.depth # make depth negative - fieldset.mindepth = -fieldset.U.depth[0] # uppermost layer in the hydrodynamic data - fieldset.add_constant("driftdepth", config.argo_characteristics["driftdepth"]) - fieldset.add_constant("maxdepth", config.argo_characteristics["maxdepth"]) - fieldset.add_constant( - "vertical_speed", config.argo_characteristics["vertical_speed"] - ) - fieldset.add_constant("cycle_days", config.argo_characteristics["cycle_days"]) - fieldset.add_constant("drift_days", config.argo_characteristics["drift_days"]) - return fieldset diff --git a/virtual_ship/drifter_deployments.py b/virtual_ship/drifter_deployments.py index 89e19039..69836bc3 100644 --- a/virtual_ship/drifter_deployments.py +++ b/virtual_ship/drifter_deployments.py @@ -7,8 +7,10 @@ import numpy as np from parcels import AdvectionRK4, FieldSet, JITParticle, ParticleSet, Variable +from .virtual_ship_configuration import VirtualShipConfiguration -def drifter_deployments(config, drifter_time): + +def drifter_deployments(config: VirtualShipConfiguration, drifter_time): """ Deploy drifters. @@ -17,7 +19,7 @@ def drifter_deployments(config, drifter_time): """ if len(config.drifter_deploylocations) > 0: - fieldset = create_drifter_fieldset(config) + fieldset = config.drifter_fieldset # Create particle to sample water underway class DrifterParticle(JITParticle): @@ -74,18 +76,18 @@ def CheckError(particle, fieldset, time): ) -def create_drifter_fieldset(config): +def create_drifter_fieldset(config, data_dir: str): """ Create a fieldset from netcdf files for drifters, returns fieldset with negative depth values. :param config: The cruise configuration. + :param data_dir: TODO :returns: The fieldset. """ - datadirname = os.path.dirname(__file__) filenames = { - "U": os.path.join(datadirname, "drifterdata_UV.nc"), - "V": os.path.join(datadirname, "drifterdata_UV.nc"), - "T": os.path.join(datadirname, "drifterdata_T.nc"), + "U": os.path.join(data_dir, "drifterdata_UV.nc"), + "V": os.path.join(data_dir, "drifterdata_UV.nc"), + "T": os.path.join(data_dir, "drifterdata_T.nc"), } variables = {"U": "uo", "V": "vo", "T": "thetao"} dimensions = { diff --git a/virtual_ship/instruments/__init__.py b/virtual_ship/instruments/__init__.py new file mode 100644 index 00000000..283e7cba --- /dev/null +++ b/virtual_ship/instruments/__init__.py @@ -0,0 +1,6 @@ +"""Measurement instrument that can be used with Parcels.""" + +from . import argo_float +from .location import Location + +__all__ = ["Location", "argo_float"] diff --git a/virtual_ship/instruments/argo_float.py b/virtual_ship/instruments/argo_float.py new file mode 100644 index 00000000..d9e7d50e --- /dev/null +++ b/virtual_ship/instruments/argo_float.py @@ -0,0 +1,173 @@ +"""Argo float instrument.""" + +import math +from dataclasses import dataclass +from datetime import timedelta + +import numpy as np +from parcels import ( + AdvectionRK4, + FieldSet, + JITParticle, + ParticleSet, + StatusCode, + Variable, +) + +from .location import Location + + +@dataclass +class ArgoFloat: + """Configuration for a single Argo float.""" + + location: Location + deployment_time: float + min_depth: float + max_depth: float + drift_depth: float + vertical_speed: float + cycle_days: float + drift_days: float + + +_ArgoParticle = JITParticle.add_variables( + [ + Variable("cycle_phase", dtype=np.int32, initial=0.0), + Variable("cycle_age", dtype=np.float32, initial=0.0), + Variable("drift_age", dtype=np.float32, initial=0.0), + Variable("salinity", initial=np.nan), + Variable("temperature", initial=np.nan), + Variable("min_depth", dtype=np.float32), + Variable("max_depth", dtype=np.float32), + Variable("drift_depth", dtype=np.float32), + Variable("vertical_speed", dtype=np.float32), + Variable("cycle_days", dtype=np.int32), + Variable("drift_days", dtype=np.int32), + ] +) + + +def _argo_float_vertical_movement(particle, fieldset, time): + if particle.cycle_phase == 0: + # Phase 0: Sinking with vertical_speed until depth is drift_depth + particle_ddepth += ( # noqa See comment above about particle_* variables. + particle.vertical_speed * particle.dt + ) + if particle.depth + particle_ddepth <= particle.drift_depth: + particle_ddepth = particle.drift_depth - particle.depth + particle.cycle_phase = 1 + + elif particle.cycle_phase == 1: + # Phase 1: Drifting at depth for drifttime seconds + particle.drift_age += particle.dt + if particle.drift_age >= particle.drift_days * 86400: + particle.drift_age = 0 # reset drift_age for next cycle + particle.cycle_phase = 2 + + elif particle.cycle_phase == 2: + # Phase 2: Sinking further to max_depth + particle_ddepth += particle.vertical_speed * particle.dt + if particle.depth + particle_ddepth <= particle.max_depth: + particle_ddepth = particle.max_depth - particle.depth + particle.cycle_phase = 3 + + elif particle.cycle_phase == 3: + # Phase 3: Rising with vertical_speed until at surface + particle_ddepth -= particle.vertical_speed * particle.dt + particle.cycle_age += ( + particle.dt + ) # solve issue of not updating cycle_age during ascent + if particle.depth + particle_ddepth >= particle.min_depth: + particle_ddepth = particle.min_depth - particle.depth + particle.temperature = ( + math.nan + ) # reset temperature to NaN at end of sampling cycle + particle.salinity = math.nan # idem + particle.cycle_phase = 4 + else: + particle.temperature = fieldset.T[ + time, particle.depth, particle.lat, particle.lon + ] + particle.salinity = fieldset.S[ + time, particle.depth, particle.lat, particle.lon + ] + + elif particle.cycle_phase == 4: + # Phase 4: Transmitting at surface until cycletime is reached + if particle.cycle_age > particle.cycle_days * 86400: + particle.cycle_phase = 0 + particle.cycle_age = 0 + + if particle.state == StatusCode.Evaluate: + particle.cycle_age += particle.dt # update cycle_age + + +def _keep_at_surface(particle, fieldset, time): + # Prevent error when float reaches surface + if particle.state == StatusCode.ErrorThroughSurface: + particle.depth = particle.min_depth + particle.state = StatusCode.Success + + +def _check_error(particle, fieldset, time): + if particle.state >= 50: # This captures all Errors + particle.delete() + + +def simulate_argo_floats( + argo_floats: list[ArgoFloat], + fieldset: FieldSet, + out_file_name: str, + outputdt: timedelta, +) -> None: + """ + Use parcels to simulate a set of Argo floats in a fieldset. + + :param argo_floats: A list of Argo floats to simulate. + :param fieldset: The fieldset to simulate the Argo floats in. + :param out_file_name: The file to write the results to. + :param outputdt: Interval which dictates the update frequency of file output during simulation + """ + lon = [argo.location.lon for argo in argo_floats] + lat = [argo.location.lat for argo in argo_floats] + time = [argo.deployment_time for argo in argo_floats] + + # define the parcels simulation + argo_float_particleset = ParticleSet( + fieldset=fieldset, + pclass=_ArgoParticle, + lon=lon, + lat=lat, + depth=[argo.min_depth for argo in argo_floats], + time=time, + min_depth=[argo.min_depth for argo in argo_floats], + max_depth=[argo.max_depth for argo in argo_floats], + drift_depth=[argo.drift_depth for argo in argo_floats], + vertical_speed=[argo.vertical_speed for argo in argo_floats], + cycle_days=[argo.cycle_days for argo in argo_floats], + drift_days=[argo.drift_days for argo in argo_floats], + ) + + # define output file for the simulation + out_file = argo_float_particleset.ParticleFile( + name=out_file_name, + outputdt=timedelta(minutes=5), + chunks=(1, 500), + ) + + # get time when the fieldset ends + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + + # execute simulation + argo_float_particleset.execute( + [ + _argo_float_vertical_movement, + AdvectionRK4, + _keep_at_surface, + _check_error, + ], + endtime=fieldset_endtime, + dt=outputdt, + output_file=out_file, + ) diff --git a/virtual_ship/instruments/location.py b/virtual_ship/instruments/location.py new file mode 100644 index 00000000..edeb4efc --- /dev/null +++ b/virtual_ship/instruments/location.py @@ -0,0 +1,29 @@ +"""Location class. See class description.""" + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class Location: + """A location on a sphere.""" + + latitude: float + longitude: float + + @property + def lat(self) -> float: + """ + Shorthand for latitude variable. + + :returns: Latitude variable. + """ + return self.latitude + + @property + def lon(self) -> float: + """ + Shorthand for longitude variable. + + :returns: Longitude variable. + """ + return self.longitude diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py index f753677d..51238069 100644 --- a/virtual_ship/sailship.py +++ b/virtual_ship/sailship.py @@ -8,16 +8,23 @@ from parcels import Field, FieldSet, JITParticle, ParticleSet, Variable from shapely.geometry import Point, Polygon +from .costs import costs +from .drifter_deployments import drifter_deployments +from .instruments.argo_float import ArgoFloat, simulate_argo_floats +from .instruments.location import Location +from .postprocess import postprocess +from .virtual_ship_configuration import VirtualShipConfiguration -def sailship(config): + +def sailship(config: VirtualShipConfiguration): """ Use parcels to simulate the ship, take CTDs and measure ADCP and underwaydata. :param config: The cruise configuration. - :returns: drifter_time, argo_time, total_time + :raises NotImplementedError: TODO """ # Create fieldset and retreive final schip route as sample_lons and sample_lats - fieldset = create_fieldset(config) + fieldset = config.ctd_fieldset # create_fieldset(config, data_dir) sample_lons, sample_lats = shiproute(config) print("Arrived in region of interest, starting to gather data.") @@ -132,7 +139,9 @@ def SampleT(particle, fieldset, time): drifter = 0 drifter_time = [] argo = 0 - argo_time = [] + argo_floats: list[ArgoFloat] = [] + + ARGO_MIN_DEPTH = -config.argo_float_fieldset.U.depth[0] # run the model for the length of the sample_lons list for i in range(len(sample_lons) - 1): @@ -155,7 +164,8 @@ def SampleT(particle, fieldset, time): print( "Ship time is over, waiting for drifters and/or Argo floats to finish." ) - return drifter_time, argo_time + raise NotImplementedError() + # return drifter_time, argo_time # check if virtual ship is at a CTD station if ctd < len(config.CTD_locations): @@ -214,7 +224,21 @@ def SampleT(particle, fieldset, time): abs(sample_lons[i] - config.argo_deploylocations[argo][0]) < 0.01 and abs(sample_lats[i] - config.argo_deploylocations[argo][1]) < 0.01 ): - argo_time.append(total_time) + argo_floats.append( + ArgoFloat( + location=Location( + latitude=config.argo_deploylocations[argo][0], + longitude=config.argo_deploylocations[argo][1], + ), + deployment_time=total_time, + min_depth=ARGO_MIN_DEPTH, + max_depth=config.argo_characteristics["maxdepth"], + drift_depth=config.argo_characteristics["driftdepth"], + vertical_speed=config.argo_characteristics["vertical_speed"], + cycle_days=config.argo_characteristics["cycle_days"], + drift_days=config.argo_characteristics["drift_days"], + ) + ) argo += 1 print(f"Argo {argo} deployed at {sample_lons[i]}, {sample_lats[i]}") if argo == len(config.argo_deploylocations): @@ -243,23 +267,42 @@ def SampleT(particle, fieldset, time): ) print("Cruise has ended. Please wait for drifters and/or Argo floats to finish.") - return drifter_time, argo_time, total_time + # simulate drifter deployments + drifter_deployments(config, drifter_time) + + # simulate argo deployments + simulate_argo_floats( + argo_floats=argo_floats, + fieldset=config.argo_float_fieldset, + out_file_name=os.path.join("results", "argo_floats.zarr"), + outputdt=timedelta(minutes=5), + ) + + # convert CTD data to CSV + postprocess() + + print("All data has been gathered and postprocessed, returning home.") + + cost = costs(config, total_time) + print( + f"This cruise took {timedelta(seconds=total_time)} and would have cost {cost:,.0f} euros." + ) -def create_fieldset(config): +def create_fieldset(config, data_dir: str): """ Create fieldset from netcdf files and adds bathymetry data for CTD cast, returns fieldset with negative depth values. :param config: The cruise configuration. + :param data_dir: TODO :returns: The fieldset. :raises ValueError: If downloaded data is not as expected. """ - datadirname = os.path.dirname(__file__) filenames = { - "U": os.path.join(datadirname, "studentdata_UV.nc"), - "V": os.path.join(datadirname, "studentdata_UV.nc"), - "S": os.path.join(datadirname, "studentdata_S.nc"), - "T": os.path.join(datadirname, "studentdata_T.nc"), + "U": os.path.join(data_dir, "studentdata_UV.nc"), + "V": os.path.join(data_dir, "studentdata_UV.nc"), + "S": os.path.join(data_dir, "studentdata_S.nc"), + "T": os.path.join(data_dir, "studentdata_T.nc"), } variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} dimensions = { @@ -286,7 +329,7 @@ def create_fieldset(config): fieldset.add_constant("maxtime", fieldset.U.grid.time_full[-1]) # add bathymetry data to the fieldset for CTD cast - bathymetry_file = os.path.join(datadirname, "GLO-MFC_001_024_mask_bathy.nc") + bathymetry_file = os.path.join(data_dir, "GLO-MFC_001_024_mask_bathy.nc") bathymetry_variables = ("bathymetry", "deptho") bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"} bathymetry_field = Field.from_netcdf( diff --git a/virtual_ship/virtual_ship_configuration.py b/virtual_ship/virtual_ship_configuration.py index fb796dc8..f74251b7 100644 --- a/virtual_ship/virtual_ship_configuration.py +++ b/virtual_ship/virtual_ship_configuration.py @@ -2,23 +2,40 @@ import datetime import json -import os from datetime import timedelta +from parcels import FieldSet from shapely.geometry import Point, Polygon class VirtualShipConfiguration: """Configuration of the virtual ship, initialized from a json file.""" - def __init__(self, json_file): + ctd_fieldset: FieldSet + drifter_fieldset: FieldSet + argo_float_fieldset: FieldSet + + def __init__( + self, + json_file, + ctd_fieldset: FieldSet, + drifter_fieldset: FieldSet, + argo_float_fieldset: FieldSet, + ): """ Initialize this object. :param json_file: Path to the JSON file to init from. + :param ctd_fieldset: Fieldset for CTD measurements. + :param drifter_fieldset: Fieldset for CTD measurements. + :param argo_float_fieldset: FIeldset for argo float measurements. :raises ValueError: If JSON file not valid. """ - with open(os.path.join(os.path.dirname(__file__), json_file), "r") as file: + self.ctd_fieldset = ctd_fieldset + self.drifter_fieldset = drifter_fieldset + self.argo_float_fieldset = argo_float_fieldset + + with open(json_file, "r") as file: json_input = json.loads(file.read()) for key in json_input: setattr(self, key, json_input[key]) diff --git a/virtualship.py b/virtualship.py index 08774699..6ac71b70 100644 --- a/virtualship.py +++ b/virtualship.py @@ -1,17 +1,6 @@ -from datetime import timedelta from virtual_ship.virtual_ship_configuration import VirtualShipConfiguration -from virtual_ship.costs import costs from virtual_ship.sailship import sailship -from virtual_ship.drifter_deployments import drifter_deployments -from virtual_ship.argo_deployments import argo_deployments -from virtual_ship.postprocess import postprocess -if __name__ == '__main__': - config = VirtualShipConfiguration('student_input.json') - drifter_time, argo_time, total_time = sailship(config) - drifter_deployments(config, drifter_time) - argo_deployments(config, argo_time) - postprocess() - print("All data has been gathered and postprocessed, returning home.") - cost = costs(config, total_time) - print(f"This cruise took {timedelta(seconds=total_time)} and would have cost {cost:,.0f} euros.") +if __name__ == "__main__": + config = VirtualShipConfiguration("student_input.json") + sailship(config)