Skip to content

Commit ea3789c

Browse files
iurytpre-commit-ci[bot]erikvansebille
authored
Add XBT instrument simulation based on CTD (#76)
* add new XBT instrument simulation based on CTD * add xbt to the export list * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add test for xbt casts * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove S from FieldSet.from_data * add XBT in the instruments list * Update src/virtualship/instruments/xbt.py Co-authored-by: Erik van Sebille <[email protected]> * add fall-rate equation * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * change raising error for shallow waters * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add missing : * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * change order of operations * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update src/virtualship/instruments/xbt.py Co-authored-by: Erik van Sebille <[email protected]> * keep only fall_speed * fix typo and add 2 multiplier to the fall rate equation * fix simple typo * set particle depth to max depth if it's too deep * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fixed the if statement * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Erik van Sebille <[email protected]> Co-authored-by: Erik van Sebille <[email protected]>
1 parent b081984 commit ea3789c

File tree

5 files changed

+278
-5
lines changed

5 files changed

+278
-5
lines changed

‎README.md‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,9 @@
3333

3434
VirtualShipParcels is a command line simulator allowing students to plan and conduct a virtual research expedition, receiving measurements as if they were coming from actual oceanographic instruments including:
3535

36-
- ADCP (for currents)
37-
- CTD (for conductivity, and temperature)
36+
- ADCP (currents)
37+
- CTD (conductivity and temperature)
38+
- XBT (temperature)
3839
- underwater measurements (salinity and temperature)
3940
- surface drifters
4041
- argo float deployments

‎docs/contributing.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,4 +51,4 @@ The running of these commands is useful for local development and quick iteratio
5151
When adding a dependency, make sure to modify the following files where relevant:
5252

5353
- `environment.yml` for core and development dependencies (important for the development environment, and CI)
54-
- `pyproject.toml` for core dependencies (important for the pypi package, this should propogate through automatically to `recipe/meta.yml` in the conda-forge feedstock)
54+
- `pyproject.toml` for core dependencies (important for the pypi package, this should propagate through automatically to `recipe/meta.yml` in the conda-forge feedstock)

‎src/virtualship/instruments/__init__.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
"""Measurement instrument that can be used with Parcels."""
22

3-
from . import adcp, argo_float, ctd, drifter, ship_underwater_st
3+
from . import adcp, argo_float, ctd, drifter, ship_underwater_st, xbt
44

5-
__all__ = ["adcp", "argo_float", "ctd", "drifter", "ship_underwater_st"]
5+
__all__ = ["adcp", "argo_float", "ctd", "drifter", "ship_underwater_st", "xbt"]

‎src/virtualship/instruments/xbt.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
"""XBT instrument."""
2+
3+
from dataclasses import dataclass
4+
from datetime import timedelta
5+
from pathlib import Path
6+
7+
import numpy as np
8+
from parcels import FieldSet, JITParticle, ParticleSet, Variable
9+
10+
from ..spacetime import Spacetime
11+
12+
13+
@dataclass
14+
class XBT:
15+
"""Configuration for a single XBT."""
16+
17+
spacetime: Spacetime
18+
min_depth: float
19+
max_depth: float
20+
fall_speed: float
21+
deceleration_coefficient: float
22+
23+
24+
_XBTParticle = JITParticle.add_variables(
25+
[
26+
Variable("temperature", dtype=np.float32, initial=np.nan),
27+
Variable("max_depth", dtype=np.float32),
28+
Variable("min_depth", dtype=np.float32),
29+
Variable("fall_speed", dtype=np.float32),
30+
Variable("deceleration_coefficient", dtype=np.float32),
31+
]
32+
)
33+
34+
35+
def _sample_temperature(particle, fieldset, time):
36+
particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
37+
38+
39+
def _xbt_cast(particle, fieldset, time):
40+
particle_ddepth = -particle.fall_speed * particle.dt
41+
42+
# update the fall speed from the quadractic fall-rate equation
43+
# check https://doi.org/10.5194/os-7-231-2011
44+
particle.fall_speed = (
45+
particle.fall_speed - 2 * particle.deceleration_coefficient * particle.dt
46+
)
47+
48+
# delete particle if depth is exactly max_depth
49+
if particle.depth == particle.max_depth:
50+
particle.delete()
51+
52+
# set particle depth to max depth if it's too deep
53+
if particle.depth + particle_ddepth < particle.max_depth:
54+
particle_ddepth = particle.max_depth - particle.depth
55+
56+
57+
def simulate_xbt(
58+
fieldset: FieldSet,
59+
out_path: str | Path,
60+
xbts: list[XBT],
61+
outputdt: timedelta,
62+
) -> None:
63+
"""
64+
Use Parcels to simulate a set of XBTs in a fieldset.
65+
66+
:param fieldset: The fieldset to simulate the XBTs in.
67+
:param out_path: The path to write the results to.
68+
:param xbts: A list of XBTs to simulate.
69+
:param outputdt: Interval which dictates the update frequency of file output during simulation
70+
:raises ValueError: Whenever provided XBTs, fieldset, are not compatible with this function.
71+
"""
72+
DT = 10.0 # dt of XBT simulation integrator
73+
74+
if len(xbts) == 0:
75+
print(
76+
"No XBTs provided. Parcels currently crashes when providing an empty particle set, so no XBT simulation will be done and no files will be created."
77+
)
78+
# TODO when Parcels supports it this check can be removed.
79+
return
80+
81+
fieldset_starttime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[0])
82+
fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1])
83+
84+
# deploy time for all xbts should be later than fieldset start time
85+
if not all(
86+
[np.datetime64(xbt.spacetime.time) >= fieldset_starttime for xbt in xbts]
87+
):
88+
raise ValueError("XBT deployed before fieldset starts.")
89+
90+
# depth the xbt will go to. shallowest between xbt max depth and bathymetry.
91+
max_depths = [
92+
max(
93+
xbt.max_depth,
94+
fieldset.bathymetry.eval(
95+
z=0, y=xbt.spacetime.location.lat, x=xbt.spacetime.location.lon, time=0
96+
),
97+
)
98+
for xbt in xbts
99+
]
100+
101+
# initial fall speeds
102+
initial_fall_speeds = [xbt.fall_speed for xbt in xbts]
103+
104+
# XBT depth can not be too shallow, because kernel would break.
105+
# This shallow is not useful anyway, no need to support.
106+
for max_depth, fall_speed in zip(max_depths, initial_fall_speeds, strict=False):
107+
if not max_depth <= -DT * fall_speed:
108+
raise ValueError(
109+
f"XBT max_depth or bathymetry shallower than maximum {-DT * fall_speed}"
110+
)
111+
112+
# define xbt particles
113+
xbt_particleset = ParticleSet(
114+
fieldset=fieldset,
115+
pclass=_XBTParticle,
116+
lon=[xbt.spacetime.location.lon for xbt in xbts],
117+
lat=[xbt.spacetime.location.lat for xbt in xbts],
118+
depth=[xbt.min_depth for xbt in xbts],
119+
time=[xbt.spacetime.time for xbt in xbts],
120+
max_depth=max_depths,
121+
min_depth=[xbt.min_depth for xbt in xbts],
122+
fall_speed=[xbt.fall_speed for xbt in xbts],
123+
)
124+
125+
# define output file for the simulation
126+
out_file = xbt_particleset.ParticleFile(name=out_path, outputdt=outputdt)
127+
128+
# execute simulation
129+
xbt_particleset.execute(
130+
[_sample_temperature, _xbt_cast],
131+
endtime=fieldset_endtime,
132+
dt=DT,
133+
verbose_progress=False,
134+
output_file=out_file,
135+
)
136+
137+
# there should be no particles left, as they delete themselves when they finish profiling
138+
if len(xbt_particleset.particledata) != 0:
139+
raise ValueError(
140+
"Simulation ended before XBT finished profiling. This most likely means the field time dimension did not match the simulation time span."
141+
)

‎tests/instruments/test_xbt.py‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
"""
2+
Test the simulation of XBT instruments.
3+
4+
Fields are kept static over time and time component of XBT measurements is not tested tested because it's tricky to provide expected measurements.
5+
"""
6+
7+
import datetime
8+
from datetime import timedelta
9+
10+
import numpy as np
11+
import xarray as xr
12+
from parcels import Field, FieldSet
13+
14+
from virtualship import Location, Spacetime
15+
from virtualship.instruments.xbt import XBT, simulate_xbt
16+
17+
18+
def test_simulate_xbts(tmpdir) -> None:
19+
# arbitrary time offset for the dummy fieldset
20+
base_time = datetime.datetime.strptime("1950-01-01", "%Y-%m-%d")
21+
22+
# where to cast XBTs
23+
xbts = [
24+
XBT(
25+
spacetime=Spacetime(
26+
location=Location(latitude=0, longitude=1),
27+
time=base_time + datetime.timedelta(hours=0),
28+
),
29+
min_depth=0,
30+
max_depth=float("-inf"),
31+
fall_speed=6.553,
32+
deceleration_coefficient=0.00242,
33+
),
34+
XBT(
35+
spacetime=Spacetime(
36+
location=Location(latitude=1, longitude=0),
37+
time=base_time,
38+
),
39+
min_depth=0,
40+
max_depth=float("-inf"),
41+
fall_speed=6.553,
42+
deceleration_coefficient=0.00242,
43+
),
44+
]
45+
46+
# expected observations for xbts at surface and at maximum depth
47+
xbt_exp = [
48+
{
49+
"surface": {
50+
"temperature": 6,
51+
"lat": xbts[0].spacetime.location.lat,
52+
"lon": xbts[0].spacetime.location.lon,
53+
},
54+
"maxdepth": {
55+
"temperature": 8,
56+
"lat": xbts[0].spacetime.location.lat,
57+
"lon": xbts[0].spacetime.location.lon,
58+
},
59+
},
60+
{
61+
"surface": {
62+
"temperature": 6,
63+
"lat": xbts[1].spacetime.location.lat,
64+
"lon": xbts[1].spacetime.location.lon,
65+
},
66+
"maxdepth": {
67+
"temperature": 8,
68+
"lat": xbts[1].spacetime.location.lat,
69+
"lon": xbts[1].spacetime.location.lon,
70+
},
71+
},
72+
]
73+
74+
# create fieldset based on the expected observations
75+
# indices are time, depth, latitude, longitude
76+
u = np.zeros((2, 2, 2, 2))
77+
v = np.zeros((2, 2, 2, 2))
78+
t = np.zeros((2, 2, 2, 2))
79+
80+
t[:, 1, 0, 1] = xbt_exp[0]["surface"]["temperature"]
81+
t[:, 0, 0, 1] = xbt_exp[0]["maxdepth"]["temperature"]
82+
t[:, 1, 1, 0] = xbt_exp[1]["surface"]["temperature"]
83+
t[:, 0, 1, 0] = xbt_exp[1]["maxdepth"]["temperature"]
84+
85+
fieldset = FieldSet.from_data(
86+
{"V": v, "U": u, "T": t},
87+
{
88+
"time": [
89+
np.datetime64(base_time + datetime.timedelta(hours=0)),
90+
np.datetime64(base_time + datetime.timedelta(hours=1)),
91+
],
92+
"depth": [-1000, 0],
93+
"lat": [0, 1],
94+
"lon": [0, 1],
95+
},
96+
)
97+
fieldset.add_field(Field("bathymetry", [-1000], lon=0, lat=0))
98+
99+
# perform simulation
100+
out_path = tmpdir.join("out.zarr")
101+
102+
simulate_xbt(
103+
xbts=xbts,
104+
fieldset=fieldset,
105+
out_path=out_path,
106+
outputdt=timedelta(seconds=10),
107+
)
108+
109+
# test if output is as expected
110+
results = xr.open_zarr(out_path)
111+
112+
assert len(results.trajectory) == len(xbts)
113+
114+
for xbt_i, (traj, exp_bothloc) in enumerate(
115+
zip(results.trajectory, xbt_exp, strict=True)
116+
):
117+
obs_surface = results.sel(trajectory=traj, obs=0)
118+
min_index = np.argmin(results.sel(trajectory=traj)["z"].data)
119+
obs_maxdepth = results.sel(trajectory=traj, obs=min_index)
120+
121+
for obs, loc in [
122+
(obs_surface, "surface"),
123+
(obs_maxdepth, "maxdepth"),
124+
]:
125+
exp = exp_bothloc[loc]
126+
for var in ["temperature", "lat", "lon"]:
127+
obs_value = obs[var].values.item()
128+
exp_value = exp[var]
129+
assert np.isclose(
130+
obs_value, exp_value
131+
), f"Observation incorrect {xbt_i=} {loc=} {var=} {obs_value=} {exp_value=}."

0 commit comments

Comments
 (0)