From 09ee8762102c28e39ab3cc1ac25063237db3ae7c Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:10:58 +0100 Subject: [PATCH 1/7] use bathymetry criterion to raise argos in shallow waters --- src/virtualship/instruments/argo_float.py | 40 +++++++++++++++++++---- src/virtualship/instruments/base.py | 2 +- 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index 204f0b3d..860383f3 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -4,6 +4,7 @@ from typing import ClassVar import numpy as np + from parcels import ( AdvectionRK4, JITParticle, @@ -11,7 +12,6 @@ StatusCode, Variable, ) - from virtualship.instruments.base import Instrument from virtualship.instruments.types import InstrumentType from virtualship.models.spacetime import Spacetime @@ -53,6 +53,7 @@ class ArgoFloat: Variable("vertical_speed", dtype=np.float32), Variable("cycle_days", dtype=np.int32), Variable("drift_days", dtype=np.int32), + Variable("grounded", dtype=np.int32, initial=0), ] ) @@ -64,10 +65,24 @@ class ArgoFloat: 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 Parcels defines particle_* variables, which code checkers cannot know. + particle_ddepth += ( # noqa particle.vertical_speed * particle.dt ) - if particle.depth + particle_ddepth <= particle.drift_depth: + + # bathymetry at particle location + loc_bathy = fieldset.bathymetry.eval( + time, particle.depth, particle.lat, particle.lon + ) + if particle.depth + particle_ddepth <= loc_bathy: + particle_ddepth = loc_bathy - particle.depth + 50.0 # 50m above bathy + particle.cycle_phase = 1 + particle.grounded = 1 + # TODO: print statments not working properly with JIT compiler... + # print( + # f"Argo float grounded at bathymetry depth: {loc_bathy}m during sinking to drift depth. Location: ({particle.lat}, {particle.lon}). Raising by 50m above bathymetry and continuing cycle." + # ) + + elif particle.depth + particle_ddepth <= particle.drift_depth: particle_ddepth = particle.drift_depth - particle.depth particle.cycle_phase = 1 @@ -81,7 +96,17 @@ def _argo_float_vertical_movement(particle, fieldset, time): 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: + loc_bathy = fieldset.bathymetry.eval( + time, particle.depth, particle.lat, particle.lon + ) + if particle.depth + particle_ddepth <= loc_bathy: + particle_ddepth = loc_bathy - particle.depth + 50.0 # 50m above bathy + particle.cycle_phase = 3 + particle.grounded = 1 + # print( + # f"Argo float grounded at bathymetry depth: {loc_bathy}m during sinking to max depth. Location: ({particle.lat}, {particle.lon}). Raising by 50m above bathymetry and continuing cycle." + # ) + elif particle.depth + particle_ddepth <= particle.max_depth: particle_ddepth = particle.max_depth - particle.depth particle.cycle_phase = 3 @@ -91,6 +116,7 @@ def _argo_float_vertical_movement(particle, fieldset, time): particle.cycle_age += ( particle.dt ) # solve issue of not updating cycle_age during ascent + particle.grounded = 0 if particle.depth + particle_ddepth >= particle.min_depth: particle_ddepth = particle.min_depth - particle.depth particle.temperature = ( @@ -148,7 +174,7 @@ def __init__(self, expedition, from_data): super().__init__( expedition, variables, - add_bathymetry=False, + add_bathymetry=True, allow_time_extrapolation=False, verbose_progress=True, spacetime_buffer_size=spacetime_buffer_size, @@ -162,6 +188,8 @@ def simulate(self, measurements, out_path) -> None: OUTPUT_DT = timedelta(minutes=5) ENDTIME = None + # TODO: insert checker that none of the waypoints are under 50m depth; if so, raise error. + if len(measurements) == 0: print( "No Argo floats provided. Parcels currently crashes when providing an empty particle set, so no argo floats simulation will be done and no files will be created." @@ -191,7 +219,7 @@ def simulate(self, measurements, out_path) -> None: out_file = argo_float_particleset.ParticleFile( name=out_path, outputdt=OUTPUT_DT, - chunks=[len(argo_float_particleset), 100], + chunks=[len(argo_float_particleset), 100], # TODO: is this necessary? ) # get earliest between fieldset end time and provide end time diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py index 22b0b54a..5faf9a76 100644 --- a/src/virtualship/instruments/base.py +++ b/src/virtualship/instruments/base.py @@ -8,9 +8,9 @@ import copernicusmarine import xarray as xr -from parcels import FieldSet from yaspin import yaspin +from parcels import FieldSet from virtualship.errors import CopernicusCatalogueError from virtualship.utils import ( COPERNICUSMARINE_PHYS_VARIABLES, From a953ccf67ece8b3f7589e0bf0c2342ccce5b829e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 8 Dec 2025 14:32:03 +0000 Subject: [PATCH 2/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/virtualship/instruments/argo_float.py | 2 +- src/virtualship/instruments/base.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index 860383f3..cfbf7c06 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -4,7 +4,6 @@ from typing import ClassVar import numpy as np - from parcels import ( AdvectionRK4, JITParticle, @@ -12,6 +11,7 @@ StatusCode, Variable, ) + from virtualship.instruments.base import Instrument from virtualship.instruments.types import InstrumentType from virtualship.models.spacetime import Spacetime diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py index 5faf9a76..22b0b54a 100644 --- a/src/virtualship/instruments/base.py +++ b/src/virtualship/instruments/base.py @@ -8,9 +8,9 @@ import copernicusmarine import xarray as xr +from parcels import FieldSet from yaspin import yaspin -from parcels import FieldSet from virtualship.errors import CopernicusCatalogueError from virtualship.utils import ( COPERNICUSMARINE_PHYS_VARIABLES, From 8c2dbcef0f7daa5cab28c0ec89cf8ac2f2c9cc48 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 8 Dec 2025 16:01:11 +0100 Subject: [PATCH 3/7] add logging messages when encountering shallow waters --- src/virtualship/instruments/argo_float.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index 860383f3..ecf3f269 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -77,10 +77,9 @@ def _argo_float_vertical_movement(particle, fieldset, time): particle_ddepth = loc_bathy - particle.depth + 50.0 # 50m above bathy particle.cycle_phase = 1 particle.grounded = 1 - # TODO: print statments not working properly with JIT compiler... - # print( - # f"Argo float grounded at bathymetry depth: {loc_bathy}m during sinking to drift depth. Location: ({particle.lat}, {particle.lon}). Raising by 50m above bathymetry and continuing cycle." - # ) + print( + "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to drift depth. Raising by 50m above bathymetry and continuing cycle." + ) elif particle.depth + particle_ddepth <= particle.drift_depth: particle_ddepth = particle.drift_depth - particle.depth @@ -103,9 +102,9 @@ def _argo_float_vertical_movement(particle, fieldset, time): particle_ddepth = loc_bathy - particle.depth + 50.0 # 50m above bathy particle.cycle_phase = 3 particle.grounded = 1 - # print( - # f"Argo float grounded at bathymetry depth: {loc_bathy}m during sinking to max depth. Location: ({particle.lat}, {particle.lon}). Raising by 50m above bathymetry and continuing cycle." - # ) + print( + "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to max depth. Raising by 50m above bathymetry and continuing cycle." + ) elif particle.depth + particle_ddepth <= particle.max_depth: particle_ddepth = particle.max_depth - particle.depth particle.cycle_phase = 3 From 6c7189c133ace96d365d44b65630ba4846f3abed Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 8 Dec 2025 16:32:01 +0100 Subject: [PATCH 4/7] throw error where waypoint is too shallow for Argo --- src/virtualship/instruments/argo_float.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index ecf3f269..64b5b68f 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -187,8 +187,6 @@ def simulate(self, measurements, out_path) -> None: OUTPUT_DT = timedelta(minutes=5) ENDTIME = None - # TODO: insert checker that none of the waypoints are under 50m depth; if so, raise error. - if len(measurements) == 0: print( "No Argo floats provided. Parcels currently crashes when providing an empty particle set, so no argo floats simulation will be done and no files will be created." @@ -198,6 +196,24 @@ def simulate(self, measurements, out_path) -> None: fieldset = self.load_input_data() + breakpoint() + + shallow_waypoints = {} + for i, m in enumerate(measurements): + loc_bathy = fieldset.bathymetry.eval( + time=0, + z=0, + y=m.spacetime.location.lat, + x=m.spacetime.location.lon, + ) + if abs(loc_bathy) < 50.0: + shallow_waypoints[f"Waypoint {i + 1}"] = f"{abs(loc_bathy):.2f}m depth" + + if len(shallow_waypoints) > 0: + raise ValueError( + f"{self.__class__.__name__} cannot be deployed in waters shallower than 50m. The following waypoints are too shallow: {shallow_waypoints}." + ) + # define parcel particles argo_float_particleset = ParticleSet( fieldset=fieldset, From 40d0d6d46b789ddbaafb0d91abb54ee117f29478 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Mon, 8 Dec 2025 16:53:43 +0100 Subject: [PATCH 5/7] add bathymetry to test argo fieldset --- src/virtualship/instruments/argo_float.py | 4 +--- tests/instruments/test_argo_float.py | 12 +++++++++++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index db4f3280..f3cf3dc8 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -4,6 +4,7 @@ from typing import ClassVar import numpy as np + from parcels import ( AdvectionRK4, JITParticle, @@ -11,7 +12,6 @@ StatusCode, Variable, ) - from virtualship.instruments.base import Instrument from virtualship.instruments.types import InstrumentType from virtualship.models.spacetime import Spacetime @@ -196,8 +196,6 @@ def simulate(self, measurements, out_path) -> None: fieldset = self.load_input_data() - breakpoint() - shallow_waypoints = {} for i, m in enumerate(measurements): loc_bathy = fieldset.bathymetry.eval( diff --git a/tests/instruments/test_argo_float.py b/tests/instruments/test_argo_float.py index cbe25d76..58bdd275 100644 --- a/tests/instruments/test_argo_float.py +++ b/tests/instruments/test_argo_float.py @@ -4,8 +4,8 @@ import numpy as np import xarray as xr -from parcels import FieldSet +from parcels import FieldSet from virtualship.instruments.argo_float import ArgoFloat, ArgoFloatInstrument from virtualship.models import Location, Spacetime @@ -27,6 +27,7 @@ def test_simulate_argo_floats(tmpdir) -> None: u = np.full((2, 2, 2), 1.0) t = np.full((2, 2, 2), CONST_TEMPERATURE) s = np.full((2, 2, 2), CONST_SALINITY) + bathy = np.full((2, 2), -5000.0) fieldset = FieldSet.from_data( {"V": v, "U": u, "T": t, "S": s}, @@ -39,6 +40,15 @@ def test_simulate_argo_floats(tmpdir) -> None: ], }, ) + fieldset.add_field( + FieldSet.from_data( + {"bathymetry": bathy}, + { + "lon": np.array([0.0, 10.0]), + "lat": np.array([0.0, 10.0]), + }, + ).bathymetry + ) # argo floats to deploy argo_floats = [ From 06d279febed0232012fb9ee5f5b53d5925f58251 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:54:01 +0000 Subject: [PATCH 6/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/virtualship/instruments/argo_float.py | 2 +- tests/instruments/test_argo_float.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index f3cf3dc8..63eac0db 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -4,7 +4,6 @@ from typing import ClassVar import numpy as np - from parcels import ( AdvectionRK4, JITParticle, @@ -12,6 +11,7 @@ StatusCode, Variable, ) + from virtualship.instruments.base import Instrument from virtualship.instruments.types import InstrumentType from virtualship.models.spacetime import Spacetime diff --git a/tests/instruments/test_argo_float.py b/tests/instruments/test_argo_float.py index 58bdd275..22c9727f 100644 --- a/tests/instruments/test_argo_float.py +++ b/tests/instruments/test_argo_float.py @@ -4,8 +4,8 @@ import numpy as np import xarray as xr - from parcels import FieldSet + from virtualship.instruments.argo_float import ArgoFloat, ArgoFloatInstrument from virtualship.models import Location, Spacetime From ab8879c029d6ffdc6c60aed58af43a1e2aaec960 Mon Sep 17 00:00:00 2001 From: j-atkins <106238905+j-atkins@users.noreply.github.com> Date: Tue, 9 Dec 2025 09:09:50 +0100 Subject: [PATCH 7/7] remove todo --- src/virtualship/instruments/argo_float.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index f3cf3dc8..7d8f5126 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -232,7 +232,7 @@ def simulate(self, measurements, out_path) -> None: out_file = argo_float_particleset.ParticleFile( name=out_path, outputdt=OUTPUT_DT, - chunks=[len(argo_float_particleset), 100], # TODO: is this necessary? + chunks=[len(argo_float_particleset), 100], ) # get earliest between fieldset end time and provide end time