diff --git a/codetools.sh b/codetools.sh index a359c5d6a..8e0a4accf 100755 --- a/codetools.sh +++ b/codetools.sh @@ -20,7 +20,7 @@ flake8 ./$PACKAGE ./$TESTS echo "--------------" echo "pydocstyle" echo "--------------" -pydocstyle ./$PACKAGE +pydocstyle ./$PACKAGE ./$TESTS echo "--------------" echo "sort-all" diff --git a/tests/conftest.py b/tests/conftest.py index 01e916532..f656c38f9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,7 +1,14 @@ +"""Test configuration that is ran for every test.""" + 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): + """ + Set the working directory for each test to the directory of that test. + + :param request: - + :param monkeypatch: - + """ monkeypatch.chdir(request.fspath.dirname) diff --git a/tests/instruments/test_adcp.py b/tests/instruments/test_adcp.py new file mode 100644 index 000000000..fdbd50e3a --- /dev/null +++ b/tests/instruments/test_adcp.py @@ -0,0 +1,33 @@ +"""Test the simulation of ADCP instruments.""" + +import numpy as np +from parcels import FieldSet + +from virtual_ship import Location, Spacetime +from virtual_ship.instruments.adcp import simulate_adcp + + +def test_simulate_adcp() -> None: + MAX_DEPTH = -1000 + MIN_DEPTH = -5 + BIN_SIZE = 24 + + fieldset = FieldSet.from_data( + {"U": 0, "V": 0}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + sample_points = [Spacetime(Location(0, 0), 0)] + + simulate_adcp( + fieldset=fieldset, + out_file_name="test", + max_depth=MAX_DEPTH, + min_depth=MIN_DEPTH, + bin_size=BIN_SIZE, + sample_points=sample_points, + ) diff --git a/tests/instruments/test_argo.py b/tests/instruments/test_argo.py index 349ed051a..1b6493ea8 100644 --- a/tests/instruments/test_argo.py +++ b/tests/instruments/test_argo.py @@ -5,14 +5,14 @@ import numpy as np from parcels import FieldSet -from virtual_ship.instruments import Location +from virtual_ship import Location, Spacetime 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 + VERTICAL_SPEED = -0.10 CYCLE_DAYS = 10 DRIFT_DAYS = 9 @@ -29,12 +29,11 @@ def test_simulate_argo_floats() -> None: argo_floats = [ ArgoFloat( - location=Location(latitude=0, longitude=0), - deployment_time=0, + spacetime=Spacetime(location=Location(latitude=0, longitude=0), time=0), min_depth=min_depth, max_depth=MAX_DEPTH, drift_depth=DRIFT_DEPTH, - vertical_speed=VERTICLE_SPEED, + vertical_speed=VERTICAL_SPEED, cycle_days=CYCLE_DAYS, drift_days=DRIFT_DAYS, ) diff --git a/tests/instruments/test_ctd.py b/tests/instruments/test_ctd.py new file mode 100644 index 000000000..92c456cd3 --- /dev/null +++ b/tests/instruments/test_ctd.py @@ -0,0 +1,38 @@ +"""Test the simulation of CTD instruments.""" + +from datetime import timedelta + +import numpy as np +from parcels import FieldSet + +from virtual_ship import Location, Spacetime +from virtual_ship.instruments.ctd import CTD, simulate_ctd + + +def test_simulate_ctds() -> None: + fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "T": 0, "S": 0, "bathymetry": 100}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + min_depth = -fieldset.U.depth[0] + max_depth = -fieldset.U.depth[-1] + + ctds = [ + CTD( + spacetime=Spacetime(location=Location(latitude=0, longitude=0), time=0), + min_depth=min_depth, + max_depth=max_depth, + ) + ] + + simulate_ctd( + ctds=ctds, + fieldset=fieldset, + out_file_name="test", + outputdt=timedelta(seconds=10), + ) diff --git a/tests/instruments/test_drifters.py b/tests/instruments/test_drifters.py new file mode 100644 index 000000000..421499d60 --- /dev/null +++ b/tests/instruments/test_drifters.py @@ -0,0 +1,36 @@ +"""Test the simulation of drifters.""" + +from datetime import timedelta + +import numpy as np +from parcels import FieldSet + +from virtual_ship import Location, Spacetime +from virtual_ship.instruments.drifter import Drifter, simulate_drifters + + +def test_simulate_drifters() -> None: + fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "T": 0}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + min_depth = -fieldset.U.depth[0] + + drifters = [ + Drifter( + spacetime=Spacetime(location=Location(latitude=0, longitude=0), time=0), + min_depth=min_depth, + ) + ] + + simulate_drifters( + drifters=drifters, + fieldset=fieldset, + out_file_name="test", + outputdt=timedelta(minutes=5), + ) diff --git a/tests/instruments/test_ship_underwater_st.py b/tests/instruments/test_ship_underwater_st.py new file mode 100644 index 000000000..05de1d496 --- /dev/null +++ b/tests/instruments/test_ship_underwater_st.py @@ -0,0 +1,29 @@ +"""Test the simulation of ship salinity temperature measurements.""" + +import numpy as np +from parcels import FieldSet + +from virtual_ship import Location, Spacetime +from virtual_ship.instruments.ship_underwater_st import simulate_ship_underwater_st + + +def test_simulate_ship_underwater_st() -> None: + DEPTH = -2 + + fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "S": 0, "T": 0}, + { + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], + }, + ) + + sample_points = [Spacetime(Location(0, 0), 0)] + + simulate_ship_underwater_st( + fieldset=fieldset, + out_file_name="test", + depth=DEPTH, + sample_points=sample_points, + ) diff --git a/tests/sailship_config.json b/tests/sailship_config.json index 23cf4c1a1..1f3805dcb 100644 --- a/tests/sailship_config.json +++ b/tests/sailship_config.json @@ -19,13 +19,13 @@ "bin_size_m": 24 }, "CTD_locations": [ - [-23.071289, 63.743631] + [-23.081289, 63.743631] ], "CTD_settings": { "max_depth": "max" }, "drifter_deploylocations": [ - + [-23.081289, 63.743631] ], "argo_deploylocations": [ [-23.081289, 63.743631] diff --git a/tests/test_sailship.py b/tests/test_sailship.py index afd8a95fb..6307af92a 100644 --- a/tests/test_sailship.py +++ b/tests/test_sailship.py @@ -8,20 +8,28 @@ def test_sailship() -> None: + adcp_fieldset = FieldSet.from_data( + {"U": 0, "V": 0}, + {"lon": 0, "lat": 0}, + ) + + ship_underwater_st_fieldset = FieldSet.from_data( + {"U": 0, "V": 0, "S": 0, "T": 0}, + {"lon": 0, "lat": 0}, + ) + 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, "T": 0}, { - "U": 0, - "V": 0, + "lon": 0, + "lat": 0, + "time": [np.datetime64("1950-01-01") + np.timedelta64(632160, "h")], }, - {"lon": 0, "lat": 0}, ) argo_float_fieldset = FieldSet.from_data( @@ -35,6 +43,8 @@ def test_sailship() -> None: config = VirtualShipConfiguration( "sailship_config.json", + adcp_fieldset=adcp_fieldset, + ship_underwater_st_fieldset=ship_underwater_st_fieldset, ctd_fieldset=ctd_fieldset, drifter_fieldset=drifter_fieldset, argo_float_fieldset=argo_float_fieldset, diff --git a/virtual_ship/__init__.py b/virtual_ship/__init__.py index bf152389c..dbf175e35 100644 --- a/virtual_ship/__init__.py +++ b/virtual_ship/__init__.py @@ -1,5 +1,7 @@ """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 .location import Location +from .spacetime import Spacetime -__all__ = ["instruments", "sailship"] +__all__ = ["Location", "Spacetime", "instruments", "sailship"] diff --git a/virtual_ship/costs.py b/virtual_ship/costs.py index 6c27f81f3..9ba1d9708 100644 --- a/virtual_ship/costs.py +++ b/virtual_ship/costs.py @@ -1,19 +1,23 @@ """costs function.""" +from datetime import timedelta -def costs(config, total_time): +from .virtual_ship_configuration import VirtualShipConfiguration + + +def costs(config: VirtualShipConfiguration, total_time: timedelta): """ Calculate the cost of the virtual ship (in US$). :param config: The cruise configuration. - :param total_time: Time cruised in seconds. + :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 // 3600 + ship_cost = ship_cost_per_day / 24 * total_time.total_seconds() // 3600 argo_cost = len(config.argo_deploylocations) * argo_deploy_cost drifter_cost = len(config.drifter_deploylocations) * drifter_deploy_cost diff --git a/virtual_ship/drifter_deployments.py b/virtual_ship/drifter_deployments.py deleted file mode 100644 index 69836bc33..000000000 --- a/virtual_ship/drifter_deployments.py +++ /dev/null @@ -1,109 +0,0 @@ -"""drifter_deployments function.""" - -import datetime -import os -from datetime import timedelta - -import numpy as np -from parcels import AdvectionRK4, FieldSet, JITParticle, ParticleSet, Variable - -from .virtual_ship_configuration import VirtualShipConfiguration - - -def drifter_deployments(config: VirtualShipConfiguration, drifter_time): - """ - Deploy drifters. - - :param config: The cruise configuration. - :param drifter_time: TODO - """ - if len(config.drifter_deploylocations) > 0: - - fieldset = config.drifter_fieldset - - # Create particle to sample water underway - class DrifterParticle(JITParticle): - """Define a new particle class that samples Temperature as a surface drifter.""" - - temperature = Variable("temperature", initial=np.nan) - - # define function sampling Temperature - def SampleT(particle, fieldset, time): - particle.temperature = fieldset.T[ - time, particle.depth, particle.lat, particle.lon - ] - - def CheckError(particle, fieldset, time): - if particle.state >= 50: # This captures all Errors - particle.delete() - - # initialize drifters - lon = [] - lat = [] - for i in range(len(config.drifter_deploylocations)): - lon.append(config.drifter_deploylocations[i][0]) - lat.append(config.drifter_deploylocations[i][1]) - time = drifter_time - - # Create and execute drifter particles - pset = ParticleSet( - fieldset=fieldset, - pclass=DrifterParticle, - lon=lon, - lat=lat, - depth=np.repeat(fieldset.mindepth, len(time)), - time=time, - ) - output_file = pset.ParticleFile( - name=os.path.join("results", "Drifters.zarr"), outputdt=timedelta(hours=1) - ) - - fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) - drifter_endtime = np.array( - ( - datetime.datetime.strptime( - config.requested_ship_time["start"], "%Y-%m-%dT%H:%M:%S" - ) - + timedelta(weeks=6) - ) - ).astype("datetime64[ms]") - - pset.execute( - [AdvectionRK4, SampleT, CheckError], - endtime=min(fieldset_endtime, drifter_endtime), - dt=timedelta(minutes=5), - output_file=output_file, - ) - - -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. - """ - filenames = { - "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 = { - "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 - return fieldset diff --git a/virtual_ship/instruments/__init__.py b/virtual_ship/instruments/__init__.py index 283e7cbac..edcaa8c17 100644 --- a/virtual_ship/instruments/__init__.py +++ b/virtual_ship/instruments/__init__.py @@ -1,6 +1,5 @@ """Measurement instrument that can be used with Parcels.""" -from . import argo_float -from .location import Location +from . import adcp, argo_float, ctd, drifter, ship_underwater_st -__all__ = ["Location", "argo_float"] +__all__ = ["adcp", "argo_float", "ctd", "drifter", "ship_underwater_st"] diff --git a/virtual_ship/instruments/adcp.py b/virtual_ship/instruments/adcp.py new file mode 100644 index 000000000..2d7d032c4 --- /dev/null +++ b/virtual_ship/instruments/adcp.py @@ -0,0 +1,66 @@ +"""ADCP instrument.""" + +import numpy as np +from parcels import FieldSet, JITParticle, ParticleSet, Variable + +from ..spacetime import Spacetime + +_ADCPParticle = JITParticle.add_variables( + [ + Variable("U", dtype=np.float32, initial=np.nan), + Variable("V", dtype=np.float32, initial=np.nan), + ] +) + + +def _sample_velocity(particle, fieldset, time): + particle.U, particle.V = fieldset.UV.eval( + time, particle.depth, particle.lat, particle.lon, applyConversion=False + ) + + +def simulate_adcp( + fieldset: FieldSet, + out_file_name: str, + max_depth: float, + min_depth: float, + bin_size: float, + sample_points: list[Spacetime], +) -> None: + """ + Use parcels to simulate an ADCP in a fieldset. + + :param fieldset: The fieldset to simulate the ADCP in. + :param out_file_name: The file to write the results to. + :param max_depth: Maximum depth the ADCP can measure. + :param min_depth: Minimum depth the ADCP can measure. + :param bin_size: How many samples to take in the complete range between max_depth and min_depth. + :param sample_points: The places and times to sample at. + """ + sample_points.sort(key=lambda p: p.time) + + bins = np.arange(max_depth, min_depth, bin_size) + num_particles = len(bins) + particleset = ParticleSet.from_list( + fieldset=fieldset, + pclass=_ADCPParticle, + lon=np.full( + num_particles, 0.0 + ), # initial lat/lon are irrelevant and will be overruled later. + lat=np.full(num_particles, 0.0), + depth=bins, + time=0, # same for time + ) + + # define output file for the simulation + out_file = particleset.ParticleFile( + name=out_file_name, + ) + + for point in sample_points: + particleset.lon_nextloop[:] = point.location.lon + particleset.lat_nextloop[:] = point.location.lat + particleset.time_nextloop[:] = point.time + + particleset.execute([_sample_velocity], dt=1, runtime=1, verbose_progress=False) + out_file.write(particleset, time=particleset[0].time) diff --git a/virtual_ship/instruments/argo_float.py b/virtual_ship/instruments/argo_float.py index d9e7d50e3..cd014671c 100644 --- a/virtual_ship/instruments/argo_float.py +++ b/virtual_ship/instruments/argo_float.py @@ -14,15 +14,14 @@ Variable, ) -from .location import Location +from ..spacetime import Spacetime @dataclass class ArgoFloat: """Configuration for a single Argo float.""" - location: Location - deployment_time: float + spacetime: Spacetime min_depth: float max_depth: float drift_depth: float @@ -36,8 +35,8 @@ class ArgoFloat: 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("salinity", dtype=np.float32, initial=np.nan), + Variable("temperature", dtype=np.float32, initial=np.nan), Variable("min_depth", dtype=np.float32), Variable("max_depth", dtype=np.float32), Variable("drift_depth", dtype=np.float32), @@ -129,11 +128,11 @@ def simulate_argo_floats( :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] + lon = [argo.spacetime.location.lon for argo in argo_floats] + lat = [argo.spacetime.location.lat for argo in argo_floats] + time = [argo.spacetime.time for argo in argo_floats] - # define the parcels simulation + # define parcel particles argo_float_particleset = ParticleSet( fieldset=fieldset, pclass=_ArgoParticle, @@ -152,7 +151,7 @@ def simulate_argo_floats( # define output file for the simulation out_file = argo_float_particleset.ParticleFile( name=out_file_name, - outputdt=timedelta(minutes=5), + outputdt=outputdt, chunks=(1, 500), ) diff --git a/virtual_ship/instruments/ctd.py b/virtual_ship/instruments/ctd.py new file mode 100644 index 000000000..899909702 --- /dev/null +++ b/virtual_ship/instruments/ctd.py @@ -0,0 +1,113 @@ +"""CTD instrument.""" + +from dataclasses import dataclass +from datetime import timedelta + +import numpy as np +from parcels import FieldSet, JITParticle, ParticleSet, Variable + +from ..spacetime import Spacetime + + +@dataclass +class CTD: + """Configuration for a single CTD.""" + + spacetime: Spacetime + min_depth: float + max_depth: float + + +_CTDParticle = JITParticle.add_variables( + [ + Variable("salinity", dtype=np.float32, initial=np.nan), + Variable("temperature", dtype=np.float32, initial=np.nan), + Variable("raising", dtype=np.int32, initial=0), + Variable("max_depth", dtype=np.float32), + ] +) + + +def _sample_temperature(particle, fieldset, time): + particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] + + +def _sample_salinity(particle, fieldset, time): + particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon] + + +def _ctd_cast(particle, fieldset, time): + # Lowering and raising CTD + if ( + -fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon] + > particle.max_depth + ): + maxdepth = ( + -fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon] + 20 + ) + else: + maxdepth = particle.max_depth + winch_speed = -1.0 # sink and rise speed in m/s + + if particle.raising == 0: + # Sinking with winch_speed until near seafloor + particle_ddepth = winch_speed * particle.dt + if particle.depth <= maxdepth: + particle.raising = 1 + + if particle.raising == 1: + # Rising with winch_speed until depth is -2 m + if particle.depth < -2: + particle_ddepth = -winch_speed * particle.dt + if particle.depth + particle_ddepth >= -2: + # to break the loop ... + particle.state = 41 + print("CTD cast finished.") + + +def simulate_ctd( + ctds: list[CTD], + fieldset: FieldSet, + out_file_name: str, + outputdt: timedelta, +) -> None: + """ + Use parcels to simulate a set of CTDs in a fieldset. + + :param ctds: A list of CTDs to simulate. + :param fieldset: The fieldset to simulate the CTDs 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 = [ctd.spacetime.location.lon for ctd in ctds] + lat = [ctd.spacetime.location.lat for ctd in ctds] + time = [ctd.spacetime.time for ctd in ctds] + + # define parcel particles + ctd_particleset = ParticleSet( + fieldset=fieldset, + pclass=_CTDParticle, + lon=lon, + lat=lat, + depth=[ctd.min_depth for ctd in ctds], + time=time, + max_depth=[ctd.max_depth for ctd in ctds], + ) + + # define output file for the simulation + out_file = ctd_particleset.ParticleFile( + name=out_file_name, + outputdt=outputdt, + chunks=(1, 500), + ) + + # get time when the fieldset ends + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + + # execute simulation + ctd_particleset.execute( + [_sample_salinity, _sample_temperature, _ctd_cast], + endtime=fieldset_endtime, + dt=outputdt, + output_file=out_file, + ) diff --git a/virtual_ship/instruments/drifter.py b/virtual_ship/instruments/drifter.py new file mode 100644 index 000000000..5bed859f3 --- /dev/null +++ b/virtual_ship/instruments/drifter.py @@ -0,0 +1,80 @@ +"""Drifter instrument.""" + +from dataclasses import dataclass +from datetime import timedelta + +import numpy as np +from parcels import AdvectionRK4, FieldSet, JITParticle, ParticleSet, Variable + +from ..spacetime import Spacetime + + +@dataclass +class Drifter: + """Configuration for a single Drifter.""" + + spacetime: Spacetime + min_depth: float + + +_DrifterParticle = JITParticle.add_variables( + [ + Variable("temperature", dtype=np.float32, initial=np.nan), + ] +) + + +def _sample_temperature(particle, fieldset, time): + particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] + + +def _check_error(particle, fieldset, time): + if particle.state >= 50: # This captures all Errors + particle.delete() + + +def simulate_drifters( + drifters: list[Drifter], + fieldset: FieldSet, + out_file_name: str, + outputdt: timedelta, +) -> None: + """ + Use parcels to simulate a set of drifters in a fieldset. + + :param drifters: A list of drifters to simulate. + :param fieldset: The fieldset to simulate the drifters 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 = [drifter.spacetime.location.lon for drifter in drifters] + lat = [drifter.spacetime.location.lat for drifter in drifters] + time = [drifter.spacetime.time for drifter in drifters] + + # define parcel particles + drifter_particleset = ParticleSet( + fieldset=fieldset, + pclass=_DrifterParticle, + lon=lon, + lat=lat, + depth=[drifter.min_depth for drifter in drifters], + time=time, + ) + + # define output file for the simulation + out_file = drifter_particleset.ParticleFile( + name=out_file_name, + outputdt=outputdt, + chunks=(1, 500), + ) + + # get time when the fieldset ends + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + + # execute simulation + drifter_particleset.execute( + [AdvectionRK4, _sample_temperature, _check_error], + endtime=fieldset_endtime, + dt=outputdt, + output_file=out_file, + ) diff --git a/virtual_ship/instruments/ship_underwater_st.py b/virtual_ship/instruments/ship_underwater_st.py new file mode 100644 index 000000000..389a1e8e4 --- /dev/null +++ b/virtual_ship/instruments/ship_underwater_st.py @@ -0,0 +1,67 @@ +"""Ship salinity and temperature.""" + +import numpy as np +from parcels import FieldSet, JITParticle, ParticleSet, Variable + +from ..spacetime import Spacetime + +_ShipSTParticle = JITParticle.add_variables( + [ + Variable("salinity", dtype=np.float32, initial=np.nan), + Variable("temperature", dtype=np.float32, initial=np.nan), + ] +) + + +# define function sampling Salinity +def _sample_salinity(particle, fieldset, time): + particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon] + + +# define function sampling Temperature +def _sample_temperature(particle, fieldset, time): + particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] + + +def simulate_ship_underwater_st( + fieldset: FieldSet, + out_file_name: str, + depth: float, + sample_points: list[Spacetime], +) -> None: + """ + Use parcels to simulate underway data, measuring salinity and temperature at the given depth along the ship track in a fieldset. + + :param fieldset: The fieldset to simulate the sampling in. + :param out_file_name: The file to write the results to. + :param depth: The depth at which to measure. 0 is water surface, negative is into the water. + :param sample_points: The places and times to sample at. + """ + sample_points.sort(key=lambda p: p.time) + + particleset = ParticleSet.from_list( + fieldset=fieldset, + pclass=_ShipSTParticle, + lon=0.0, # initial lat/lon are irrelevant and will be overruled later. + lat=0.0, + depth=depth, + time=0, # same for time + ) + + # define output file for the simulation + out_file = particleset.ParticleFile( + name=out_file_name, + ) + + for point in sample_points: + particleset.lon_nextloop[:] = point.location.lon + particleset.lat_nextloop[:] = point.location.lat + particleset.time_nextloop[:] = point.time + + particleset.execute( + [_sample_salinity, _sample_temperature], + dt=1, + runtime=1, + verbose_progress=False, + ) + out_file.write(particleset, time=particleset[0].time) diff --git a/virtual_ship/instruments/location.py b/virtual_ship/location.py similarity index 100% rename from virtual_ship/instruments/location.py rename to virtual_ship/location.py diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py index 51238069b..62bf57671 100644 --- a/virtual_ship/sailship.py +++ b/virtual_ship/sailship.py @@ -5,272 +5,228 @@ import numpy as np import pyproj -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.adcp import simulate_adcp from .instruments.argo_float import ArgoFloat, simulate_argo_floats -from .instruments.location import Location +from .instruments.ctd import CTD, simulate_ctd +from .instruments.drifter import Drifter, simulate_drifters +from .instruments.ship_underwater_st import simulate_ship_underwater_st +from .location import Location from .postprocess import postprocess +from .spacetime import Spacetime from .virtual_ship_configuration import VirtualShipConfiguration def sailship(config: VirtualShipConfiguration): """ - Use parcels to simulate the ship, take CTDs and measure ADCP and underwaydata. + Use parcels to simulate the ship, take ctd_instruments and measure ADCP and underwaydata. :param config: The cruise configuration. - :raises NotImplementedError: TODO """ - # Create fieldset and retreive final schip route as sample_lons and sample_lats - 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.") - - # Create Vessel Mounted ADCP like particles to sample the ocean - class VM_ADCPParticle(JITParticle): - """Define a new particle class that does Vessel Mounted ADCP like measurements.""" - - U = Variable("U", dtype=np.float32, initial=0.0) - V = Variable("V", dtype=np.float32, initial=0.0) - - # Create particle to sample water underway - class UnderwayDataParticle(JITParticle): - """Define a new particle class that samples water directly under the hull.""" + # TODO this will be in the config later, but for now we don't change the config structure + # from here ----- + argo_locations_list = [ + Location(latitude=argo[1], longitude=argo[0]) + for argo in config.argo_deploylocations + ] + argo_locations = set(argo_locations_list) + if len(argo_locations) != len(argo_locations_list): + print( + "WARN: Some argo float deployment locations are identical and have been combined." + ) - salinity = Variable("salinity", initial=np.nan) - temperature = Variable("temperature", initial=np.nan) + drifter_locations_list = [ + Location(latitude=drifter[1], longitude=drifter[0]) + for drifter in config.drifter_deploylocations + ] + drifter_locations = set(drifter_locations_list) + if len(drifter_locations) != len(drifter_locations_list): + print( + "WARN: Some drifter deployment locations are identical and have been combined." + ) - # Create CTD like particles to sample the ocean - class CTDParticle(JITParticle): - """Define a new particle class that does CTD like measurements.""" + ctd_locations_list = [ + Location(latitude=ctd[1], longitude=ctd[0]) for ctd in config.CTD_locations + ] + ctd_locations = set(ctd_locations_list) + if len(ctd_locations) != len(ctd_locations_list): + print("WARN: Some CTD locations are identical and have been combined.") + # until here ---- - salinity = Variable("salinity", initial=np.nan) - temperature = Variable("temperature", initial=np.nan) - raising = Variable("raising", dtype=np.int32, initial=0.0) + # get discrete points along the ships route were sampling and deployments will be performed + route_dt = timedelta(minutes=5) + route_points = shiproute(config=config, dt=route_dt) - # define ADCP sampling function without conversion (because of A grid) - def SampleVel(particle, fieldset, time): - particle.U, particle.V = fieldset.UV.eval( - time, particle.depth, particle.lat, particle.lon, applyConversion=False - ) + # adcp objects to be used in simulation + adcps: list[Spacetime] = [] - # define function lowering and raising CTD - def CTDcast(particle, fieldset, time): - # TODO question: if is executed every time... move outside function? Not if "drifting" now possible - if ( - -fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon] - > fieldset.max_depth - ): - maxdepth = ( - -fieldset.bathymetry[time, particle.depth, particle.lat, particle.lon] - + 20 - ) - else: - maxdepth = fieldset.max_depth - winch_speed = -1.0 # sink and rise speed in m/s - - if particle.raising == 0: - # Sinking with winch_speed until near seafloor - particle_ddepth = winch_speed * particle.dt - if particle.depth <= maxdepth: - particle.raising = 1 - - if particle.raising == 1: - # Rising with winch_speed until depth is -2 m - if particle.depth < -2: - particle_ddepth = -winch_speed * particle.dt - if particle.depth + particle_ddepth >= -2: - # to break the loop ... - particle.state = 41 - print("CTD cast finished.") - - # define function sampling Salinity - def SampleS(particle, fieldset, time): - particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon] - - # define function sampling Temperature - def SampleT(particle, fieldset, time): - particle.temperature = fieldset.T[ - time, particle.depth, particle.lat, particle.lon - ] - - # Create ADCP like particleset and output file - ADCP_bins = np.arange( - config.ADCP_settings["max_depth"], -5, config.ADCP_settings["bin_size_m"] - ) - vert_particles = len(ADCP_bins) - pset_ADCP = ParticleSet.from_list( - fieldset=fieldset, - pclass=VM_ADCPParticle, - lon=np.full(vert_particles, sample_lons[0]), - lat=np.full(vert_particles, sample_lats[0]), - depth=ADCP_bins, - time=0, - ) - adcp_output_file = pset_ADCP.ParticleFile(name=os.path.join("results", "ADCP.zarr")) - adcp_dt = timedelta( - minutes=5 - ).total_seconds() # timestep of ADCP output, every 5 min - - # Create underway particle - pset_UnderwayData = ParticleSet.from_list( - fieldset=fieldset, - pclass=UnderwayDataParticle, - lon=sample_lons[0], - lat=sample_lats[0], - depth=-2, - time=0, - ) - UnderwayData_output_file = pset_UnderwayData.ParticleFile( - name=os.path.join("results", "UnderwayData.zarr") - ) + # ship st objects to be used in simulation + ship_underwater_sts: list[Spacetime] = [] - # initialize CTD station number and time - total_time = timedelta(hours=0).total_seconds() - ctd = 0 - ctd_dt = timedelta( - seconds=10 - ) # timestep of CTD output reflecting post-process binning into 10m bins - - # initialize drifters and argo floats - drifter = 0 - drifter_time = [] - argo = 0 + # argo float deployment locations that have been visited + argo_locations_visited: set[Location] = set() + # argo float objects to be used in simulation argo_floats: list[ArgoFloat] = [] - 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): - + # drifter deployment locations that have been visited + drifter_locations_visited: set[Location] = set() + # drifter objects to be used in simulation + drifters: list[Drifter] = [] + + # ctd cast locations that have been visited + ctd_locations_visited: set[Location] = set() + # ctd cast objects to be used in simulation + ctds: list[CTD] = [] + + # iterate over each discrete route point, find deployment and measurement locations and times, and measure how much time this took + # TODO between drifters, argo floats, ctd there is quite some repetition + print("Traversing ship route.") + time_past = timedelta() + for i, route_point in enumerate(route_points): if i % 96 == 0: - print(f"Gathered data {timedelta(seconds=total_time)} hours since start.") - - # execute the ADCP kernels to sample U and V and underway T and S - if config.ADCP_data: - pset_ADCP.execute( - [SampleVel], dt=adcp_dt, runtime=1, verbose_progress=False + print(f"Gathered data {time_past} hours since start.") + + # find drifter deployments to be done at this location + drifters_here = set( + [ + drifter + for drifter in drifter_locations - drifter_locations_visited + if all( + np.isclose( + [drifter.lat, drifter.lon], [route_point.lat, route_point.lon] + ) + ) + ] + ) + if len(drifters_here) > 1: + print( + "WARN: Multiple drifter deployments match the current location. Only a single deployment will be performed." ) - adcp_output_file.write(pset_ADCP, time=pset_ADCP[0].time) - if config.underway_data: - pset_UnderwayData.execute( - [SampleS, SampleT], dt=adcp_dt, runtime=1, verbose_progress=False + drifters.append( + Drifter( + spacetime=Spacetime( + location=route_point, time=time_past.total_seconds() + ), + min_depth=-config.drifter_fieldset.U.depth[0], ) - UnderwayData_output_file.write(pset_UnderwayData, time=pset_ADCP[0].time) - if pset_ADCP[0].time > fieldset.maxtime: + ) + drifter_locations_visited = drifter_locations_visited.union(drifters_here) + + # find argo float deployments to be done at this location + argos_here = set( + [ + argo + for argo in argo_locations - argo_locations_visited + if all( + np.isclose([argo.lat, argo.lon], [route_point.lat, route_point.lon]) + ) + ] + ) + if len(argos_here) > 1: print( - "Ship time is over, waiting for drifters and/or Argo floats to finish." + "WARN: Multiple argo float deployments match the current location. Only a single deployment will be performed." ) - raise NotImplementedError() - # return drifter_time, argo_time - - # check if virtual ship is at a CTD station - if ctd < len(config.CTD_locations): - if ( - abs(sample_lons[i] - config.CTD_locations[ctd][0]) < 0.001 - and abs(sample_lats[i] - config.CTD_locations[ctd][1]) < 0.001 - ): - ctd += 1 - - # release CTD particle - pset_CTD = ParticleSet( - fieldset=fieldset, - pclass=CTDParticle, - lon=sample_lons[i], - lat=sample_lats[i], - depth=fieldset.mindepth, - time=total_time, - ) - - # create a ParticleFile to store the CTD output - ctd_output_file = pset_CTD.ParticleFile( - name=f"{os.path.join('results', 'CTDs', 'CTD_')}{ctd:03d}.zarr", - outputdt=ctd_dt, + argo_floats.append( + ArgoFloat( + spacetime=Spacetime( + location=route_point, time=time_past.total_seconds() + ), + min_depth=-config.argo_float_fieldset.U.depth[0], + 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_locations_visited = argo_locations_visited.union(argos_here) + + # find CTD casts to be done at this location + ctds_here = set( + [ + ctd + for ctd in ctd_locations - ctd_locations_visited + if all( + np.isclose([ctd.lat, ctd.lon], [route_point.lat, route_point.lon]) ) + ] + ) + if len(ctds_here) > 1: + print( + "WARN: Multiple CTD casts match the current location. Only a single cast will be performed." + ) + ctds.append( + CTD( + spacetime=Spacetime( + location=route_point, time=time_past.total_seconds() + ), + min_depth=-config.ctd_fieldset.U.depth[0], + max_depth=-config.ctd_fieldset.U.depth[-1], + ) + ) + ctd_locations_visited = ctd_locations_visited.union(ctds_here) + # add 20 minutes to sailing time for deployment + if len(ctds_here) != 0: + time_past += timedelta(minutes=20) + + # add time it takes to move to the next route point + time_past += route_dt + # remove the last one, because no sailing to the next point was needed + time_past -= route_dt + + # check if all drifter, argo float, ctd locations were visited + if len(drifter_locations_visited) != len(drifter_locations): + print( + "WARN: some drifter deployments were not planned along the route and have not been performed." + ) - # record the temperature and salinity of the particle - pset_CTD.execute( - [SampleS, SampleT, CTDcast], - runtime=timedelta(hours=8), - dt=ctd_dt, - output_file=ctd_output_file, - verbose_progress=False, - ) - total_time = ( - pset_CTD.time[0] + timedelta(minutes=20).total_seconds() - ) # add CTD time and 20 minutes for deployment - - # check if we are at a drifter deployment location - if drifter < len(config.drifter_deploylocations): - while ( - abs(sample_lons[i] - config.drifter_deploylocations[drifter][0]) < 0.01 - and abs(sample_lats[i] - config.drifter_deploylocations[drifter][1]) - < 0.01 - ): - drifter_time.append(total_time) - drifter += 1 - print( - f"Drifter {drifter} deployed at {sample_lons[i]}, {sample_lats[i]}" - ) - if drifter == len(config.drifter_deploylocations): - break - - # check if we are at a argo deployment location - if argo < len(config.argo_deploylocations): - while ( - abs(sample_lons[i] - config.argo_deploylocations[argo][0]) < 0.01 - and abs(sample_lats[i] - config.argo_deploylocations[argo][1]) < 0.01 - ): - 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): - break - - # update the particle time and location - pset_ADCP.lon_nextloop[:] = sample_lons[i + 1] - pset_ADCP.lat_nextloop[:] = sample_lats[i + 1] - pset_UnderwayData.lon_nextloop[:] = sample_lons[i + 1] - pset_UnderwayData.lat_nextloop[:] = sample_lats[i + 1] - - total_time += adcp_dt - pset_ADCP.time_nextloop[:] = total_time - pset_UnderwayData.time_nextloop[:] = total_time - - # write the final locations of the ADCP and Underway data particles - if config.ADCP_data: - pset_ADCP.execute(SampleVel, dt=adcp_dt, runtime=1, verbose_progress=False) - adcp_output_file.write_latest_locations(pset_ADCP, time=total_time) - if config.underway_data: - pset_UnderwayData.execute( - [SampleS, SampleT], dt=adcp_dt, runtime=1, verbose_progress=False + if len(argo_locations_visited) != len(argo_locations): + print( + "WARN: some argo float deployments were not planned along the route and have not been performed." ) - UnderwayData_output_file.write_latest_locations( - pset_UnderwayData, time=total_time + + if len(ctd_locations_visited) != len(ctd_locations): + print( + "WARN: some CTD casts were not planned along the route and have not been performed." ) - print("Cruise has ended. Please wait for drifters and/or Argo floats to finish.") - # simulate drifter deployments - drifter_deployments(config, drifter_time) + print("Simulating onboard salinity and temperature measurements.") + simulate_ship_underwater_st( + fieldset=config.ship_underwater_st_fieldset, + out_file_name=os.path.join("results", "ship_underwater_st.zarr"), + depth=-2, + sample_points=ship_underwater_sts, + ) + + print("Simulating onboard ADCP.") + simulate_adcp( + fieldset=config.adcp_fieldset, + out_file_name=os.path.join("results", "adcp.zarr"), + max_depth=config.ADCP_settings["max_depth"], + min_depth=-5, + bin_size=config.ADCP_settings["bin_size_m"], + sample_points=adcps, + ) - # simulate argo deployments + print("Simulating CTD casts.") + simulate_ctd( + ctds=ctds, + fieldset=config.ctd_fieldset, + out_file_name=os.path.join("results", "ctd.zarr"), + outputdt=timedelta(seconds=10), + ) + + print("Simulating drifters") + simulate_drifters( + drifters=drifters, + fieldset=config.drifter_fieldset, + out_file_name=os.path.join("results", "drifters.zarr"), + outputdt=timedelta(minutes=5), + ) + + print("Simulating argo floats") simulate_argo_floats( argo_floats=argo_floats, fieldset=config.argo_float_fieldset, @@ -279,90 +235,21 @@ def SampleT(particle, fieldset, time): ) # convert CTD data to CSV + print("Postprocessing..") 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, 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. - """ - filenames = { - "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 = { - "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" - 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 - if config.CTD_settings["max_depth"] == "max": - fieldset.add_constant("max_depth", -fieldset.U.depth[-1]) - else: - fieldset.add_constant("max_depth", config.CTD_settings["max_depth"]) - fieldset.add_constant("maxtime", fieldset.U.grid.time_full[-1]) - - # add bathymetry data to the fieldset for CTD cast - 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( - bathymetry_file, bathymetry_variables, bathymetry_dimensions - ) - fieldset.add_field(bathymetry_field) - # read in data already - fieldset.computeTimeChunk(0, 1) - - if fieldset.U.lon.min() > config.region_of_interest["West"]: - raise ValueError( - "FieldSet western boundary is outside region of interest. Please run download_data.py again." - ) - if fieldset.U.lon.max() < config.region_of_interest["East"]: - raise ValueError( - "FieldSet eastern boundary is outside region of interest. Please run download_data.py again." - ) - if fieldset.U.lat.min() > config.region_of_interest["South"]: - raise ValueError( - "FieldSet southern boundary is outside region of interest. Please run download_data.py again." - ) - if fieldset.U.lat.max() < config.region_of_interest["North"]: - raise ValueError( - "FieldSet northern boundary is outside region of interest. Please run download_data.py again." - ) - return fieldset + cost = costs(config, time_past) + print(f"This cruise took {time_past} and would have cost {cost:,.0f} euros.") -def shiproute(config): +def shiproute(config: VirtualShipConfiguration, dt: timedelta) -> list[Location]: """ Take in route coordinates and return lat and lon points within region of interest to sample. :param config: The cruise configuration. + :param dt: Sailing time between each discrete route point. :returns: lat and lon points within region of interest to sample. """ # Initialize lists to store intermediate points @@ -382,7 +269,7 @@ def shiproute(config): cruise_speed = 5.14 geod = pyproj.Geod(ellps="WGS84") azimuth1, azimuth2, distance = geod.inv(startlong, startlat, endlong, endlat) - if distance > (cruise_speed * 60 * 5): + if distance > (cruise_speed * dt.total_seconds()): r = geod.inv_intermediate( startlong, startlat, @@ -414,4 +301,8 @@ def shiproute(config): if poly.contains(Point(lons[i], lats[i])): sample_lons.append(lons[i]) sample_lats.append(lats[i]) - return sample_lons, sample_lats + points = [ + Location(latitude=lat, longitude=lon) + for lat, lon in zip(sample_lats, sample_lons, strict=True) + ] + return points diff --git a/virtual_ship/spacetime.py b/virtual_ship/spacetime.py new file mode 100644 index 000000000..c277e6426 --- /dev/null +++ b/virtual_ship/spacetime.py @@ -0,0 +1,14 @@ +"""Location class. See class description.""" + +from dataclasses import dataclass + +from .location import Location + + +@dataclass +# TODO I take suggestions for a better name +class Spacetime: + """A location and time.""" + + location: Location + time: float diff --git a/virtual_ship/virtual_ship_configuration.py b/virtual_ship/virtual_ship_configuration.py index f74251b7b..ca2b8547a 100644 --- a/virtual_ship/virtual_ship_configuration.py +++ b/virtual_ship/virtual_ship_configuration.py @@ -11,6 +11,10 @@ class VirtualShipConfiguration: """Configuration of the virtual ship, initialized from a json file.""" + adcp_fieldset: FieldSet + ship_underwater_st_fieldset: FieldSet + argo_float_fieldset: FieldSet + drifter_fieldset: FieldSet ctd_fieldset: FieldSet drifter_fieldset: FieldSet argo_float_fieldset: FieldSet @@ -18,6 +22,8 @@ class VirtualShipConfiguration: def __init__( self, json_file, + adcp_fieldset: FieldSet, + ship_underwater_st_fieldset: FieldSet, ctd_fieldset: FieldSet, drifter_fieldset: FieldSet, argo_float_fieldset: FieldSet, @@ -26,11 +32,15 @@ def __init__( Initialize this object. :param json_file: Path to the JSON file to init from. + :param adcp_fieldset: Fieldset for ADCP measurements. + :param ship_underwater_st_fieldset: Fieldset for ship salinity temperature measurements. :param ctd_fieldset: Fieldset for CTD measurements. :param drifter_fieldset: Fieldset for CTD measurements. - :param argo_float_fieldset: FIeldset for argo float measurements. + :param argo_float_fieldset: Fieldset for argo float measurements. :raises ValueError: If JSON file not valid. """ + self.adcp_fieldset = adcp_fieldset + self.ship_underwater_st_fieldset = ship_underwater_st_fieldset self.ctd_fieldset = ctd_fieldset self.drifter_fieldset = drifter_fieldset self.argo_float_fieldset = argo_float_fieldset