From 16437aa30a23eea3b52baa9b108c6659a7ef79a0 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 23 Jul 2026 11:43:00 +0200
Subject: [PATCH 01/25] first changes; Parcels v4 API

---
 src/virtualship/instruments/adcp.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 26c8122d..92756f41 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -41,7 +41,6 @@ def _sample_velocity(particles, fieldset):
         particles.z,
         particles.x,
         particles.y,
-        applyConversion=False,
     )
 
 
@@ -99,15 +98,18 @@ def simulate(self, measurements, out_path) -> None:
             adcp_config.sensors, _ADCP_NONSENSOR_VARIABLES
         )
 
+        times = [np.datetime64(point.time) for point in measurements]
+        lons = [point.location.lon for point in measurements]
+        lats = [point.location.lat for point in measurements]
+        depth_full = ...  # noqa
+
         bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS)
-        num_particles = len(bins)
         particleset = ParticleSet(
             fieldset=fieldset,
             pclass=_ADCPParticle,
-            x=np.full(
-                num_particles, 0.0
-            ),  # initial lat/lon are irrelevant and will be overruled later
-            y=np.full(num_particles, 0.0),
+            t=times,
+            x=lons,
+            y=lats,
             z=bins,
         )
 

From d7a5d58c9ebcb2895b278bc4b365c45097838336 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 23 Jul 2026 14:50:55 +0200
Subject: [PATCH 02/25] move to using field.eval for UV sampling, plus call to
 custom write-to-parquet

---
 src/virtualship/instruments/adcp.py | 86 +++++++++++------------------
 1 file changed, 32 insertions(+), 54 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 92756f41..57d2e352 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -3,12 +3,11 @@
 from typing import ClassVar
 
 import numpy as np
-from parcels import ParticleFile, ParticleSet
 
 from virtualship.instruments.base import FetchSpec, Instrument
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
-from virtualship.utils import build_particle_class_from_sensors, register_instrument
+from virtualship.utils import _write_underway_to_parquet, register_instrument
 
 # =====================================================
 # SECTION: Dataclass
@@ -22,14 +21,6 @@ class ADCP:
     name: ClassVar[str] = "ADCP"
 
 
-# =====================================================
-# SECTION: non-sensor Particle Variables (non-sampling)
-# =====================================================
-
-# ADCP has no non-sensor variables, only sensor variables.
-_ADCP_NONSENSOR_VARIABLES: list = []
-
-
 # =====================================================
 # SECTION: Kernels
 # =====================================================
@@ -92,50 +83,37 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # build dynamic particle class from the active sensors
-        adcp_config = self.expedition.instruments_config.adcp_config
-        _ADCPParticle = build_particle_class_from_sensors(
-            adcp_config.sensors, _ADCP_NONSENSOR_VARIABLES
+        # use first active field for time reference
+        _time_ref_key = next(iter(self.variables))
+        _time_ref_field = getattr(fieldset, _time_ref_key)
+        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
+
+        # times in seconds since fieldset time origin, expanded across depth bins
+        times = np.array(
+            [
+                (np.datetime64(point.time) - fieldset_starttime)
+                / np.timedelta64(1, "s")
+                for point in measurements
+            ]
         )
-
-        times = [np.datetime64(point.time) for point in measurements]
-        lons = [point.location.lon for point in measurements]
-        lats = [point.location.lat for point in measurements]
-        depth_full = ...  # noqa
-
+        lons = np.array([point.location.lon for point in measurements])
+        lats = np.array([point.location.lat for point in measurements])
         bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS)
-        particleset = ParticleSet(
-            fieldset=fieldset,
-            pclass=_ADCPParticle,
-            t=times,
-            x=lons,
-            y=lats,
-            z=bins,
-        )
 
-        out_file = ParticleFile(path=out_path, outputdt=np.inf)
-
-        # build kernel list from active sensors only
-        sampling_kernels = [
-            self.sensor_kernels[sc.sensor_type]
-            for sc in adcp_config.sensors
-            if sc.enabled and sc.sensor_type in self.sensor_kernels
-        ]
-
-        # TODO: need to overhaul ADCP/underway instruments generally... don't think this Parcels API works anymore
-        # TODO: a good time to implement https://github.com/Parcels-code/virtualship/issues/231
-
-        for point in measurements:
-            particleset.lon_nextloop[:] = point.location.lon
-            particleset.lat_nextloop[:] = point.location.lat
-            particleset.time_nextloop[:] = fieldset.time_origin.reltime(
-                np.datetime64(point.time)
-            )
-
-            particleset.execute(
-                sampling_kernels,
-                dt=1,
-                runtime=1,
-                verbose_progress=self.verbose_progress,
-                output_file=out_file,
-            )
+        times_full = np.repeat(times, NUM_BINS)
+        lons_full = np.repeat(lons, NUM_BINS)
+        lats_full = np.repeat(lats, NUM_BINS)
+        depths_full = np.tile(bins, len(times))
+
+        u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
+
+        _write_underway_to_parquet(
+            dat_arrays=[u, v],
+            var_names=self.variables.keys(),
+            times_full=times_full,
+            lons_full=lons_full,
+            lats_full=lats_full,
+            depths_full=depths_full,
+            fieldset_time_origin=fieldset_starttime,
+            out_path=out_path,
+        )

From fa178c21e2a99590f82418449552eed257330cca Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 23 Jul 2026 14:51:22 +0200
Subject: [PATCH 03/25] custom write-to-parquet for underway instruments

---
 src/virtualship/utils.py | 86 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)

diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py
index 9cb028f3..97220e1c 100644
--- a/src/virtualship/utils.py
+++ b/src/virtualship/utils.py
@@ -14,6 +14,8 @@
 import copernicusmarine
 import numpy as np
 import parcels
+import pyarrow as pa
+import pyarrow.parquet as pq
 import pyproj
 import xarray as xr
 from parcels import FieldSet, Particle, Variable
@@ -563,6 +565,90 @@ def _find_files_in_timerange(
     return [fname for _, fname in files_with_dates]
 
 
+def _write_underway_to_parquet(
+    dat_arrays: list[np.ndarray],
+    var_names: list[str],
+    times_full: np.ndarray,
+    lons_full: np.ndarray,
+    lats_full: np.ndarray,
+    fieldset_time_origin: np.datetime64,
+    out_path: Path | str,
+    depths_full: np.ndarray | None = None,
+    depth_fill: float = np.nan,
+    compression: Literal["zstd", "gzip", "snappy", "brotli", None] = "zstd",
+) -> None:
+    """
+    Write underway instrument data to a Parquet file mirroring the Parcels v4 ParticleFile schema.
+
+    Designed so that output files can be re-read back in with Parcels.read_particlefile for downstream workflows.
+    """
+    breakpoint()
+    assert len(dat_arrays) == len(var_names), (
+        "dat_arrays and var_names must have the same length"
+    )
+
+    n = len(times_full)
+
+    origin_str = str(fieldset_time_origin).replace("T", " ")
+    t_metadata = {"units": f"seconds since {origin_str}", "calendar": "standard"}
+
+    # base schema mirroring Parcels ParticleFile schema, not yet with sampled variables
+    base_schema = pa.schema(
+        [
+            pa.field("t", pa.float64(), metadata=t_metadata),
+            pa.field("z", pa.float32()),
+            pa.field("y", pa.float32()),
+            pa.field("x", pa.float32()),
+            pa.field("particle_id", pa.int64()),
+            pa.field("dt", pa.float64()),
+            pa.field("state", pa.int32()),
+        ],
+        metadata={
+            "feature_type": "trajectory",
+            "Conventions": "CF-1.6/CF-1.7",
+            "ncei_template_version": "NCEI_NetCDF_Trajectory_Template_v2.0",
+            "parcels_version": parcels.__version__,
+            "parcels_grid_mesh": "spherical",  # TODO: is the case as long as using Copernicus Marine data...
+        },
+    )
+
+    for var in var_names:
+        base_schema = base_schema.append(
+            pa.field(var, pa.float32())
+        )  # add sampled variable to schema
+
+    out_path = Path(out_path)
+    if out_path.suffix != ".parquet":
+        raise ValueError(f"out_path must end in '.parquet', got {out_path.suffix!r}")
+
+    if depths_full is None:
+        depths_full = np.full(n, depth_fill, dtype=np.float32)
+
+    # build table with all data, including sampled variables
+    table = pa.table(
+        {
+            "t": pa.array(times_full.astype(np.float64)),
+            "z": pa.array(depths_full.astype(np.float32))
+            if depths_full is not None
+            else pa.array(np.full(n, depth_fill, dtype=np.float32)),
+            "y": pa.array(lats_full.astype(np.float32)),
+            "x": pa.array(lons_full.astype(np.float32)),
+            "particle_id": pa.array(
+                np.arange(n, dtype=np.int64)
+            ),  # doesn't necessarily correspond to an actual particle_id in the simulation, but required for Parcels ParticleFile schema
+            "dt": pa.array(np.full(n, np.nan, dtype=np.float64)),
+            "state": pa.array(np.zeros(n, dtype=np.int32)),
+            **{
+                var: pa.array(dat.astype(np.float32))
+                for var, dat in zip(var_names, dat_arrays, strict=True)
+            },
+        },
+        schema=base_schema,
+    )
+
+    pq.write_table(table, out_path, compression=compression)
+
+
 def _compute_max_depths(measurements, fieldset) -> list[float]:
     """Compute the effective max depth for each measurement, capped by bathymetry."""
     return [

From a14a4f2541f2df405935d4a7d3bd3ad9e6914724 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Mon, 27 Jul 2026 10:32:13 +0200
Subject: [PATCH 04/25] instrument type property, conditional windowed arrays

---
 src/virtualship/instruments/base.py | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index bbb32889..968ab708 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -15,8 +15,10 @@
 from yaspin import yaspin
 
 from virtualship.errors import CopernicusCatalogueError
+from virtualship.instruments.types import InstrumentType
 from virtualship.utils import (
     COPERNICUSMARINE_PHYS_VARIABLES,
+    INSTRUMENT_CLASS_MAP,
     _find_files_in_timerange,
     _find_nc_file_with_variable,
     _get_bathy_data,
@@ -130,7 +132,7 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = True  # TODO: just for dev; remove before merging
+        TMP = False  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
@@ -239,7 +241,11 @@ def _generate_fieldset(self) -> parcels.FieldSet:
             ds_fset = self._via_tmp_ds(ds_fset)
 
             fs = parcels.FieldSet.from_sgrid_conventions(ds_fset)
-            fs.to_windowed_arrays()  # always to windowed arrays, just in case any ds is Dask backed
+
+            # non-underway instruments to windowed arrays, just in case any ds is Dask backed
+            # underway instruments not to windowed arrays, they use fieldset.eval() which could cause a big memory usage if the fieldset is windowed
+            if not self.instrument_type.is_underway:
+                fs.to_windowed_arrays()
 
             fieldsets_list.append(fs)
 
@@ -268,3 +274,8 @@ def _via_tmp_ds(ds) -> xr.Dataset:
         ds.to_netcdf(tmp_fpath)
         del ds
         return xr.open_dataset(tmp_fpath)
+
+    @property
+    def instrument_type(self) -> InstrumentType:
+        """Return the InstrumentType for this instrument instance."""
+        return next(k for k, v in INSTRUMENT_CLASS_MAP.items() if type(self) is v)

From 9251cca6dc371b7108f54c27f26651c83bf00a84 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Mon, 27 Jul 2026 11:05:20 +0200
Subject: [PATCH 05/25] tidy up

---
 src/virtualship/utils.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py
index 97220e1c..26563bb3 100644
--- a/src/virtualship/utils.py
+++ b/src/virtualship/utils.py
@@ -582,7 +582,6 @@ def _write_underway_to_parquet(
 
     Designed so that output files can be re-read back in with Parcels.read_particlefile for downstream workflows.
     """
-    breakpoint()
     assert len(dat_arrays) == len(var_names), (
         "dat_arrays and var_names must have the same length"
     )

From 610da5f33a8a43fa6e0dd73cd832dff136041c88 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Tue, 28 Jul 2026 11:33:25 +0200
Subject: [PATCH 06/25] sample individual U and V fields for m s-1 units

---
 src/virtualship/instruments/adcp.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 57d2e352..21d4be8d 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -96,6 +96,7 @@ def simulate(self, measurements, out_path) -> None:
                 for point in measurements
             ]
         )
+
         lons = np.array([point.location.lon for point in measurements])
         lats = np.array([point.location.lat for point in measurements])
         bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS)
@@ -105,7 +106,10 @@ def simulate(self, measurements, out_path) -> None:
         lats_full = np.repeat(lats, NUM_BINS)
         depths_full = np.tile(bins, len(times))
 
-        u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
+        u, v = (
+            fieldset.U.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full),
+            fieldset.V.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full),
+        )
 
         _write_underway_to_parquet(
             dat_arrays=[u, v],

From 436d4f3094f5c855ec39c2d7cc85a79eef40e301 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 10:39:24 +0200
Subject: [PATCH 07/25] fix u, v sampling to use correct sampling and
 conversions

---
 src/virtualship/instruments/adcp.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 21d4be8d..65788c90 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -106,10 +106,9 @@ def simulate(self, measurements, out_path) -> None:
         lats_full = np.repeat(lats, NUM_BINS)
         depths_full = np.tile(bins, len(times))
 
-        u, v = (
-            fieldset.U.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full),
-            fieldset.V.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full),
-        )
+        u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
+        u = u * 1852 * 60 * np.cos(np.deg2rad(lats_full))
+        v = v * 1852 * 60
 
         _write_underway_to_parquet(
             dat_arrays=[u, v],

From b7119810cde4a38085c2757bb050d5763a202047 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 11:59:55 +0200
Subject: [PATCH 08/25] new intermediate UnderwayInstrument base class (parquet
 writing moved over from utils), change kernels to special-case underway
 kernels

---
 src/virtualship/instruments/adcp.py |  73 ++++++++--------
 src/virtualship/instruments/base.py | 131 +++++++++++++++++++++++++++-
 src/virtualship/utils.py            |  85 ------------------
 3 files changed, 164 insertions(+), 125 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 65788c90..07f196be 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -1,38 +1,35 @@
 from collections.abc import Callable
-from dataclasses import dataclass
 from typing import ClassVar
 
 import numpy as np
 
-from virtualship.instruments.base import FetchSpec, Instrument
+from virtualship.instruments.base import (
+    FetchSpec,
+    UnderwayCoordinates,
+    UnderwayInstrument,
+)
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
-from virtualship.utils import _write_underway_to_parquet, register_instrument
+from virtualship.utils import register_instrument
 
 # =====================================================
-# SECTION: Dataclass
+# SECTION: Kernels
 # =====================================================
 
 
-@dataclass
-class ADCP:
-    """ADCP configuration."""
-
-    name: ClassVar[str] = "ADCP"
+# N.B. underway 'kernels' are special cases, where the particleset is not needed, and the kernel is not passed to `pset.execute()` as would be done for a typical Parcels workflow.
+# Instead, the 'kernel' function is used only once to evaluate the fieldset at given times, depths, lats, lons.
 
 
-# =====================================================
-# SECTION: Kernels
-# =====================================================
+def _sample_underway_velocity(fieldset, times_full, depths_full, lats_full, lons_full):
+    # eval
+    u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
 
+    # convert from degrees s-1 to metres s-1
+    u = u * 1852 * 60 * np.cos(np.deg2rad(lats_full))
+    v = v * 1852 * 60
 
-def _sample_velocity(particles, fieldset):
-    particles.U, particles.V = fieldset.UV.eval(
-        particles.t,
-        particles.z,
-        particles.x,
-        particles.y,
-    )
+    return u, v
 
 
 # =====================================================
@@ -41,11 +38,11 @@ def _sample_velocity(particles, fieldset):
 
 
 @register_instrument(InstrumentType.ADCP)
-class ADCPInstrument(Instrument):
+class ADCPInstrument(UnderwayInstrument):
     """ADCP instrument class."""
 
     sensor_kernels: ClassVar[dict[SensorType, Callable]] = {
-        SensorType.VELOCITY: _sample_velocity,
+        SensorType.VELOCITY: _sample_underway_velocity,
     }
 
     def __init__(self, expedition, from_data):
@@ -64,9 +61,9 @@ def __init__(self, expedition, from_data):
 
     def simulate(self, measurements, out_path) -> None:
         """Simulate ADCP measurements."""
-        config_max_depth = (
-            self.expedition.instruments_config.adcp_config.max_depth_meter
-        )
+        adcp_config = self.expedition.instruments_config.adcp_config
+
+        config_max_depth = adcp_config.max_depth_meter
 
         if config_max_depth < -1600.0:
             print(
@@ -77,7 +74,7 @@ def simulate(self, measurements, out_path) -> None:
 
         MAX_DEPTH = config_max_depth
         MIN_DEPTH = -5.0
-        NUM_BINS = self.expedition.instruments_config.adcp_config.num_bins
+        NUM_BINS = adcp_config.num_bins
 
         measurements.sort(key=lambda p: p.time)
 
@@ -101,22 +98,22 @@ def simulate(self, measurements, out_path) -> None:
         lats = np.array([point.location.lat for point in measurements])
         bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS)
 
-        times_full = np.repeat(times, NUM_BINS)
-        lons_full = np.repeat(lons, NUM_BINS)
-        lats_full = np.repeat(lats, NUM_BINS)
-        depths_full = np.tile(bins, len(times))
+        # full sampling coordinates
+        coords = UnderwayCoordinates(
+            times=np.repeat(times, NUM_BINS),
+            lons=np.repeat(lons, NUM_BINS),
+            lats=np.repeat(lats, NUM_BINS),
+            depths=np.tile(bins, len(times)),
+        )
 
-        u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
-        u = u * 1852 * 60 * np.cos(np.deg2rad(lats_full))
-        v = v * 1852 * 60
+        sampled = self._sample_underway(
+            config_sensors=adcp_config.sensors, fieldset=fieldset, coords=coords
+        )
 
-        _write_underway_to_parquet(
-            dat_arrays=[u, v],
+        self._write_underway_to_parquet(
+            dat_arrays=sampled,
             var_names=self.variables.keys(),
-            times_full=times_full,
-            lons_full=lons_full,
-            lats_full=lats_full,
-            depths_full=depths_full,
             fieldset_time_origin=fieldset_starttime,
             out_path=out_path,
+            coords=coords,
         )
diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 6df3b513..7c4473f7 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -7,10 +7,13 @@
 from datetime import timedelta
 from itertools import pairwise
 from pathlib import Path
-from typing import TYPE_CHECKING, ClassVar
+from typing import TYPE_CHECKING, ClassVar, Literal
 
 import copernicusmarine
+import numpy as np
 import parcels
+import pyarrow as pa
+import pyarrow.parquet as pq
 import xarray as xr
 from yaspin import yaspin
 
@@ -132,7 +135,7 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = False  # TODO: just for dev; remove before merging
+        TMP = True  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
@@ -279,3 +282,127 @@ def _via_tmp_ds(ds) -> xr.Dataset:
     def instrument_type(self) -> InstrumentType:
         """Return the InstrumentType for this instrument instance."""
         return next(k for k, v in INSTRUMENT_CLASS_MAP.items() if type(self) is v)
+
+
+@dataclass()
+class UnderwayCoordinates(frozen=True):
+    """1D evaluation grid arrays for underway instruments."""
+
+    times: np.ndarray  # seconds since origin
+    lons: np.ndarray
+    lats: np.ndarray
+    depths: np.ndarray | None = None
+
+
+class UnderwayInstrument(Instrument):
+    """Intermediate base class for underway instruments, which perform variable sampling without ParticleSets."""
+
+    def _sample_underway(
+        self,
+        config_sensors: list,
+        fieldset: parcels.FieldSet,
+        coords: UnderwayCoordinates,
+    ):
+        """Perform variable sampling for underway instruments and their active sensors."""
+        sampling_kernels = [
+            self.sensor_kernels[sc.sensor_type]
+            for sc in config_sensors
+            if sc.enabled and sc.sensor_type in self.sensor_kernels
+        ]  # active sensors only
+
+        sampled = [
+            kernel(fieldset, coords) for kernel in sampling_kernels
+        ]  # perform sampling
+
+        # ensure that sampled is a flat list of arrays, even if some kernels return tuples/lists of arrays
+        sampled_flat = [
+            arr
+            for item in sampled
+            for arr in (item if isinstance(item, (tuple, list)) else (item,))
+        ]
+
+        return sampled_flat
+
+    @staticmethod
+    def _write_underway_to_parquet(
+        dat_arrays: list[np.ndarray],
+        var_names: list[str],
+        fieldset_time_origin: np.datetime64,
+        out_path: Path | str,
+        coords: UnderwayCoordinates,
+        depth_fill: float = np.nan,
+        compression: Literal["zstd", "gzip", "snappy", "brotli", None] = "zstd",
+    ) -> None:
+        """
+        Write underway instrument data to a Parquet file mirroring the Parcels v4 ParticleFile schema.
+
+        Designed so that output files can be re-read back in with Parcels.read_particlefile for downstream workflows.
+        """
+        assert len(dat_arrays) == len(var_names), (
+            "dat_arrays and var_names must have the same length"
+        )
+
+        n = len(coords.times)
+
+        origin_str = str(fieldset_time_origin).replace("T", " ")
+        t_metadata = {"units": f"seconds since {origin_str}", "calendar": "standard"}
+
+        # base schema mirroring Parcels ParticleFile schema, not yet with sampled variables
+        base_schema = pa.schema(
+            [
+                pa.field("t", pa.float64(), metadata=t_metadata),
+                pa.field("z", pa.float32()),
+                pa.field("y", pa.float32()),
+                pa.field("x", pa.float32()),
+                pa.field("particle_id", pa.int64()),
+                pa.field("dt", pa.float64()),
+                pa.field("state", pa.int32()),
+            ],
+            metadata={
+                "feature_type": "trajectory",
+                "Conventions": "CF-1.6/CF-1.7",
+                "ncei_template_version": "NCEI_NetCDF_Trajectory_Template_v2.0",
+                "parcels_version": parcels.__version__,
+                "parcels_grid_mesh": "spherical",  # TODO: is the case as long as using Copernicus Marine data...
+            },
+        )
+
+        for var in var_names:
+            base_schema = base_schema.append(
+                pa.field(var, pa.float32())
+            )  # add sampled variable to schema
+
+        out_path = Path(out_path)
+        if out_path.suffix != ".parquet":
+            raise ValueError(
+                f"out_path must end in '.parquet', got {out_path.suffix!r}"
+            )
+
+        if coords.depths is None:
+            depths_full = np.full(n, depth_fill, dtype=np.float32)
+        else:
+            depths_full = coords.depths
+
+        # build table with all data, including sampled variables
+        table = pa.table(
+            {
+                "t": pa.array(coords.times.astype(np.float64)),
+                "z": pa.array(depths_full.astype(np.float32))
+                if depths_full is not None
+                else pa.array(np.full(n, depth_fill, dtype=np.float32)),
+                "y": pa.array(coords.lats.astype(np.float32)),
+                "x": pa.array(coords.lons.astype(np.float32)),
+                "particle_id": pa.array(
+                    np.arange(n, dtype=np.int64)
+                ),  # doesn't necessarily correspond to an actual particle_id in the simulation, but required for Parcels ParticleFile schema
+                "dt": pa.array(np.full(n, np.nan, dtype=np.float64)),
+                "state": pa.array(np.zeros(n, dtype=np.int32)),
+                **{
+                    var: pa.array(dat.astype(np.float32))
+                    for var, dat in zip(var_names, dat_arrays, strict=True)
+                },
+            },
+            schema=base_schema,
+        )
+
+        pq.write_table(table, out_path, compression=compression)
diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py
index 26563bb3..9cb028f3 100644
--- a/src/virtualship/utils.py
+++ b/src/virtualship/utils.py
@@ -14,8 +14,6 @@
 import copernicusmarine
 import numpy as np
 import parcels
-import pyarrow as pa
-import pyarrow.parquet as pq
 import pyproj
 import xarray as xr
 from parcels import FieldSet, Particle, Variable
@@ -565,89 +563,6 @@ def _find_files_in_timerange(
     return [fname for _, fname in files_with_dates]
 
 
-def _write_underway_to_parquet(
-    dat_arrays: list[np.ndarray],
-    var_names: list[str],
-    times_full: np.ndarray,
-    lons_full: np.ndarray,
-    lats_full: np.ndarray,
-    fieldset_time_origin: np.datetime64,
-    out_path: Path | str,
-    depths_full: np.ndarray | None = None,
-    depth_fill: float = np.nan,
-    compression: Literal["zstd", "gzip", "snappy", "brotli", None] = "zstd",
-) -> None:
-    """
-    Write underway instrument data to a Parquet file mirroring the Parcels v4 ParticleFile schema.
-
-    Designed so that output files can be re-read back in with Parcels.read_particlefile for downstream workflows.
-    """
-    assert len(dat_arrays) == len(var_names), (
-        "dat_arrays and var_names must have the same length"
-    )
-
-    n = len(times_full)
-
-    origin_str = str(fieldset_time_origin).replace("T", " ")
-    t_metadata = {"units": f"seconds since {origin_str}", "calendar": "standard"}
-
-    # base schema mirroring Parcels ParticleFile schema, not yet with sampled variables
-    base_schema = pa.schema(
-        [
-            pa.field("t", pa.float64(), metadata=t_metadata),
-            pa.field("z", pa.float32()),
-            pa.field("y", pa.float32()),
-            pa.field("x", pa.float32()),
-            pa.field("particle_id", pa.int64()),
-            pa.field("dt", pa.float64()),
-            pa.field("state", pa.int32()),
-        ],
-        metadata={
-            "feature_type": "trajectory",
-            "Conventions": "CF-1.6/CF-1.7",
-            "ncei_template_version": "NCEI_NetCDF_Trajectory_Template_v2.0",
-            "parcels_version": parcels.__version__,
-            "parcels_grid_mesh": "spherical",  # TODO: is the case as long as using Copernicus Marine data...
-        },
-    )
-
-    for var in var_names:
-        base_schema = base_schema.append(
-            pa.field(var, pa.float32())
-        )  # add sampled variable to schema
-
-    out_path = Path(out_path)
-    if out_path.suffix != ".parquet":
-        raise ValueError(f"out_path must end in '.parquet', got {out_path.suffix!r}")
-
-    if depths_full is None:
-        depths_full = np.full(n, depth_fill, dtype=np.float32)
-
-    # build table with all data, including sampled variables
-    table = pa.table(
-        {
-            "t": pa.array(times_full.astype(np.float64)),
-            "z": pa.array(depths_full.astype(np.float32))
-            if depths_full is not None
-            else pa.array(np.full(n, depth_fill, dtype=np.float32)),
-            "y": pa.array(lats_full.astype(np.float32)),
-            "x": pa.array(lons_full.astype(np.float32)),
-            "particle_id": pa.array(
-                np.arange(n, dtype=np.int64)
-            ),  # doesn't necessarily correspond to an actual particle_id in the simulation, but required for Parcels ParticleFile schema
-            "dt": pa.array(np.full(n, np.nan, dtype=np.float64)),
-            "state": pa.array(np.zeros(n, dtype=np.int32)),
-            **{
-                var: pa.array(dat.astype(np.float32))
-                for var, dat in zip(var_names, dat_arrays, strict=True)
-            },
-        },
-        schema=base_schema,
-    )
-
-    pq.write_table(table, out_path, compression=compression)
-
-
 def _compute_max_depths(measurements, fieldset) -> list[float]:
     """Compute the effective max depth for each measurement, capped by bathymetry."""
     return [

From a35628273fddb55a0f87dd233044d589e7581e19 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 12:22:28 +0200
Subject: [PATCH 09/25] fix intermediate class logic so that sensor_kernels
 check is only triggered for final instrument child classes

---
 src/virtualship/instruments/base.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 7c4473f7..779d3192 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -2,6 +2,7 @@
 
 import abc
 import collections
+import inspect
 import tempfile
 from dataclasses import dataclass
 from datetime import timedelta
@@ -53,8 +54,11 @@ class Instrument(abc.ABC):
     sensor_kernels: ClassVar[dict[SensorType, collections.abc.Callable]]
 
     def __init_subclass__(cls, **kwargs: object) -> None:
-        """Ensure subclasses define sensor_kernels as class attribute."""
+        """Ensure non-abstract subclasses (i.e. final instrument classes) define sensor_kernels as a class attribute."""
         super().__init_subclass__(**kwargs)
+        if inspect.isabstract(cls):
+            return
+
         if "sensor_kernels" not in cls.__dict__:
             raise TypeError(
                 f"Instrument subclass '{cls.__name__}' must define 'sensor_kernels' as a class attribute."
@@ -135,7 +139,7 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = True  # TODO: just for dev; remove before merging
+        TMP = False  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
@@ -284,8 +288,8 @@ def instrument_type(self) -> InstrumentType:
         return next(k for k, v in INSTRUMENT_CLASS_MAP.items() if type(self) is v)
 
 
-@dataclass()
-class UnderwayCoordinates(frozen=True):
+@dataclass(frozen=True)
+class UnderwayCoordinates:
     """1D evaluation grid arrays for underway instruments."""
 
     times: np.ndarray  # seconds since origin

From 1614a7c8f29ce2b3f3bbbf2318984cccb3d50053 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 13:33:25 +0200
Subject: [PATCH 10/25] remove support for depth as None

---
 src/virtualship/instruments/base.py | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 779d3192..3b0af330 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -334,7 +334,6 @@ def _write_underway_to_parquet(
         fieldset_time_origin: np.datetime64,
         out_path: Path | str,
         coords: UnderwayCoordinates,
-        depth_fill: float = np.nan,
         compression: Literal["zstd", "gzip", "snappy", "brotli", None] = "zstd",
     ) -> None:
         """
@@ -382,18 +381,13 @@ def _write_underway_to_parquet(
                 f"out_path must end in '.parquet', got {out_path.suffix!r}"
             )
 
-        if coords.depths is None:
-            depths_full = np.full(n, depth_fill, dtype=np.float32)
-        else:
-            depths_full = coords.depths
-
         # build table with all data, including sampled variables
         table = pa.table(
             {
                 "t": pa.array(coords.times.astype(np.float64)),
-                "z": pa.array(depths_full.astype(np.float32))
-                if depths_full is not None
-                else pa.array(np.full(n, depth_fill, dtype=np.float32)),
+                "z": pa.array(coords.depths.astype(np.float32))
+                if coords.depths is not None
+                else pa.array(np.full(n, np.nan, dtype=np.float32)),
                 "y": pa.array(coords.lats.astype(np.float32)),
                 "x": pa.array(coords.lons.astype(np.float32)),
                 "particle_id": pa.array(

From 40784c8e440b2562c551baec64fa22015363398c Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 14:36:58 +0200
Subject: [PATCH 11/25] use coords for adcp 'kernel'

---
 src/virtualship/instruments/adcp.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 07f196be..23a21c76 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -2,6 +2,7 @@
 from typing import ClassVar
 
 import numpy as np
+import parcels
 
 from virtualship.instruments.base import (
     FetchSpec,
@@ -21,12 +22,14 @@
 # Instead, the 'kernel' function is used only once to evaluate the fieldset at given times, depths, lats, lons.
 
 
-def _sample_underway_velocity(fieldset, times_full, depths_full, lats_full, lons_full):
+def _sample_underway_velocity(fieldset: parcels.FieldSet, coords: UnderwayCoordinates):
     # eval
-    u, v = fieldset.UV.eval(t=times_full, z=depths_full, x=lons_full, y=lats_full)
+    u, v = fieldset.UV.eval(
+        t=coords.times, z=coords.depths, x=coords.lons, y=coords.lats
+    )
 
     # convert from degrees s-1 to metres s-1
-    u = u * 1852 * 60 * np.cos(np.deg2rad(lats_full))
+    u = u * 1852 * 60 * np.cos(np.deg2rad(coords.lats))
     v = v * 1852 * 60
 
     return u, v

From 0b7c0defa34ccfbe62482a4c3e6ed52dfe9c22b1 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 14:37:44 +0200
Subject: [PATCH 12/25] add validtion to UnderwayCoordinates class, rename func

---
 src/virtualship/instruments/base.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 3b0af330..08bdf5d3 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -295,7 +295,28 @@ class UnderwayCoordinates:
     times: np.ndarray  # seconds since origin
     lons: np.ndarray
     lats: np.ndarray
-    depths: np.ndarray | None = None
+    depths: np.ndarray
+
+    def __post_init__(self):
+        """Validate that all arrays are 1D and have the same length."""
+        shapes = {
+            "times": self.times.shape,
+            "lons": self.lons.shape,
+            "lats": self.lats.shape,
+            "depths": self.depths.shape,
+        }
+
+        for name, shape in shapes.items():
+            if len(shape) != 1:
+                raise ValueError(f"Array '{name}' must be 1D, but got shape {shape}.")
+
+        n = len(self.times)
+        if not (len(self.lons) == len(self.lats) == len(self.depths) == n):
+            raise ValueError(
+                f"Array length mismatch in UnderwayCoordinates: "
+                f"times={len(self.times)}, lons={len(self.lons)}, "
+                f"lats={len(self.lats)}, depths={len(self.depths)}"
+            )
 
 
 class UnderwayInstrument(Instrument):
@@ -328,7 +349,7 @@ def _sample_underway(
         return sampled_flat
 
     @staticmethod
-    def _write_underway_to_parquet(
+    def _to_parquet(
         dat_arrays: list[np.ndarray],
         var_names: list[str],
         fieldset_time_origin: np.datetime64,

From 2c32b9c2c029b555040f561a41c85bda5d519c8a Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 14:38:19 +0200
Subject: [PATCH 13/25] migrate UnderwaterST instrument to new underway
 instrument logic

---
 .../instruments/ship_underwater_st.py         | 119 ++++++++----------
 1 file changed, 50 insertions(+), 69 deletions(-)

diff --git a/src/virtualship/instruments/ship_underwater_st.py b/src/virtualship/instruments/ship_underwater_st.py
index 1dc7522a..3519c2c8 100644
--- a/src/virtualship/instruments/ship_underwater_st.py
+++ b/src/virtualship/instruments/ship_underwater_st.py
@@ -1,51 +1,40 @@
 from collections.abc import Callable
-from dataclasses import dataclass
 from typing import ClassVar
 
 import numpy as np
-from parcels import ParticleFile, ParticleSet
+import parcels
 
-from virtualship.instruments.base import FetchSpec, Instrument
+from virtualship.instruments.base import (
+    FetchSpec,
+    UnderwayCoordinates,
+    UnderwayInstrument,
+)
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
 from virtualship.utils import (
-    build_particle_class_from_sensors,
     register_instrument,
 )
 
-# =====================================================
-# SECTION: Dataclass
-# =====================================================
-
-
-@dataclass
-class Underwater_ST:
-    """Underwater_ST configuration."""
-
-    name: ClassVar[str] = "Underwater_ST"
-
-
-# =====================================================
-# SECTION: non-sensor Particle Variables (non-sampling)
-# =====================================================
-
-# Underwater ST has no non-sensor variables, only sensor variables.
-_ST_NONSENSOR_VARIABLES: list = []
-
-
 # =====================================================
 # SECTION: Kernels
 # =====================================================
 
+# N.B. underway 'kernels' are special cases, where the particleset is not needed, and the kernel is not passed to `pset.execute()` as would be done for a typical Parcels workflow.
+# Instead, the 'kernel' function is used only once to evaluate the fieldset at given times, depths, lats, lons.
 
-# define function sampling Salinity
-def _sample_salinity(particles, fieldset):
-    particles.S = fieldset.S[particles.t, particles.z, particles.y, particles.x]
 
+def _sample_underway_salinity(fieldset: parcels.FieldSet, coords: UnderwayCoordinates):
+    return fieldset.S.eval(
+        t=coords.times, z=coords.depths, x=coords.lons, y=coords.lats
+    )
 
-# define function sampling Temperature
-def _sample_temperature(particles, fieldset):
-    particles.T = fieldset.T[particles.t, particles.z, particles.y, particles.x]
+
+def _sample_underway_temperature(
+    fieldset: parcels.FieldSet, coords: UnderwayCoordinates
+):
+    return fieldset.T.eval(
+        t=coords.times, z=coords.depths, x=coords.lons, y=coords.lats
+    )
 
 
 # =====================================================
@@ -54,12 +43,12 @@ def _sample_temperature(particles, fieldset):
 
 
 @register_instrument(InstrumentType.UNDERWATER_ST)
-class Underwater_STInstrument(Instrument):
+class Underwater_STInstrument(UnderwayInstrument):
     """Underwater_ST instrument class."""
 
     sensor_kernels: ClassVar[dict[SensorType, Callable]] = {
-        SensorType.TEMPERATURE: _sample_temperature,
-        SensorType.SALINITY: _sample_salinity,
+        SensorType.TEMPERATURE: _sample_underway_temperature,
+        SensorType.SALINITY: _sample_underway_salinity,
     }
 
     def __init__(self, expedition, from_data):
@@ -80,49 +69,41 @@ def __init__(self, expedition, from_data):
 
     def simulate(self, measurements, out_path) -> None:
         """Simulate underway salinity and temperature measurements."""
+        st_config = self.expedition.instruments_config.ship_underwater_st_config
+
         DEPTH = -2.0
 
         measurements.sort(key=lambda p: p.time)
 
         fieldset = self.load_input_data()
 
-        # build dynamic particle class from the active sensors
-        st_config = self.expedition.instruments_config.ship_underwater_st_config
-        _ShipSTParticle = build_particle_class_from_sensors(
-            st_config.sensors, _ST_NONSENSOR_VARIABLES
+        # use first active field for time reference
+        _time_ref_key = next(iter(self.variables))
+        _time_ref_field = getattr(fieldset, _time_ref_key)
+        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
+
+        # sampling times and locations
+        times = np.array(
+            [
+                (np.datetime64(point.time) - fieldset_starttime)
+                / np.timedelta64(1, "s")
+                for point in measurements
+            ]
+        )
+        lons = np.array([point.location.lon for point in measurements])
+        lats = np.array([point.location.lat for point in measurements])
+        coords = UnderwayCoordinates(
+            times, lons, lats, depths=np.full_like(times, DEPTH)
         )
 
-        particleset = ParticleSet(
-            fieldset=fieldset,
-            pclass=_ShipSTParticle,
-            x=0.0,
-            y=0.0,
-            z=DEPTH,
+        sampled = self._sample_underway(
+            config_sensors=st_config.sensors, fieldset=fieldset, coords=coords
         )
 
-        out_file = ParticleFile(path=out_path, outputdt=np.inf)
-
-        # build kernel list from active sensors only
-        sampling_kernels = [
-            self.sensor_kernels[sc.sensor_type]
-            for sc in st_config.sensors
-            if sc.enabled and sc.sensor_type in self.sensor_kernels
-        ]
-
-        # TODO: need to overhaul UNDERWATER_ST/underway instruments generally... don't think this Parcels API works anymore
-        # TODO: a good time to implement https://github.com/Parcels-code/virtualship/issues/231
-
-        for point in measurements:
-            particleset.lon_nextloop[:] = point.location.lon
-            particleset.lat_nextloop[:] = point.location.lat
-            particleset.time_nextloop[:] = fieldset.time_origin.reltime(
-                np.datetime64(point.time)
-            )
-
-            particleset.execute(
-                sampling_kernels,
-                dt=1,
-                runtime=1,
-                verbose_progress=self.verbose_progress,
-                output_file=out_file,
-            )
+        self._to_parquet(
+            dat_arrays=sampled,
+            var_names=self.variables.keys(),
+            fieldset_time_origin=fieldset_starttime,
+            out_path=out_path,
+            coords=coords,
+        )

From 511d06b2bd896b51319bc4bc1425ba01d4cb7286 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 14:44:59 +0200
Subject: [PATCH 14/25] fix func name

---
 src/virtualship/instruments/adcp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 23a21c76..4fe0459f 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -113,7 +113,7 @@ def simulate(self, measurements, out_path) -> None:
             config_sensors=adcp_config.sensors, fieldset=fieldset, coords=coords
         )
 
-        self._write_underway_to_parquet(
+        self._to_parquet(
             dat_arrays=sampled,
             var_names=self.variables.keys(),
             fieldset_time_origin=fieldset_starttime,

From c86c78ad30dab4a3cc55b149e6238a0d2ebb32dd Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Wed, 29 Jul 2026 15:41:23 +0200
Subject: [PATCH 15/25] add new tests and refine for intermediate
 UnderwayInstrument class

---
 tests/instruments/test_base.py | 190 ++++++++++++++++++++++++++++++++-
 1 file changed, 189 insertions(+), 1 deletion(-)

diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py
index 93a38e90..1d03e1ec 100644
--- a/tests/instruments/test_base.py
+++ b/tests/instruments/test_base.py
@@ -1,12 +1,27 @@
+from dataclasses import dataclass
+from typing import ClassVar
 from unittest.mock import MagicMock, patch
 
+import numpy as np
+import parcels
+import pyarrow.parquet as pq
 import pytest
 import xarray as xr
 
-from virtualship.instruments.base import FetchSpec, Instrument
+from virtualship.instruments.base import (
+    FetchSpec,
+    Instrument,
+    UnderwayCoordinates,
+    UnderwayInstrument,
+)
+from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
 from virtualship.utils import get_instrument_class
 
+# =============================================================================
+# Instrument base class testing
+# =============================================================================
+
 
 def test_FetchSpec():
     fetch_spec = FetchSpec()
@@ -193,3 +208,176 @@ def test_instrument_subclass_without_sensor_kernels_error():
         class ErrorInstrument(Instrument):
             def simulate(self, data_dir, measurements, out_path):
                 pass
+
+
+# =============================================================================
+# UnderwayInstrument intermediate class testing
+# =============================================================================
+
+
+@dataclass
+class DummySensorConfig:
+    """Mock sensor configuration."""
+
+    sensor_type: SensorType
+    enabled: bool = True
+
+
+class ConcreteUnderwayInstrument(UnderwayInstrument):
+    """Concrete subclass of UnderwayInstrument for testing."""
+
+    sensor_kernels: ClassVar = {
+        SensorType.TEMPERATURE: lambda fieldset, coords: np.array(
+            [15.0, 16.0], dtype=np.float32
+        ),
+        SensorType.SALINITY: lambda fieldset, coords: np.array(
+            [35.0, 35.1], dtype=np.float32
+        ),
+        SensorType.VELOCITY: lambda fieldset, coords: (
+            np.array([0.5, 0.6], dtype=np.float32),  # U vector component
+            np.array([-0.1, -0.2], dtype=np.float32),  # V vector component
+        ),
+    }
+
+    def simulate(self, measurements, out_path) -> None:  # noqa
+        pass
+
+
+@pytest.fixture
+def sample_underway_coords():
+    """Fixture providing valid 1D UnderwayCoordinates."""
+    return UnderwayCoordinates(
+        times=np.array([0.0, 3600.0]),
+        lons=np.array([-5.0, -5.1]),
+        lats=np.array([50.0, 50.1]),
+        depths=np.array([-2.0, -2.0]),
+    )
+
+
+@pytest.fixture
+def dummy_underway_inst():
+    """Bypass __init__ and requirements for expedition object etc. for testing."""
+    return ConcreteUnderwayInstrument.__new__(ConcreteUnderwayInstrument)
+
+
+def test_underway_coordinates_validation():
+    """UnderwayCoordinates validates array lengths upon instantiation."""
+    # valid coordinates work cleanly
+    coords = UnderwayCoordinates(
+        times=np.array([0.0, 1.0]),
+        lons=np.array([10.0, 11.0]),
+        lats=np.array([20.0, 21.0]),
+        depths=np.array([-1.0, -1.0]),
+    )
+    assert len(coords.times) == 2
+
+    # mismatched array lengths raise ValueError
+    with pytest.raises(ValueError, match="Array length mismatch"):
+        UnderwayCoordinates(
+            times=np.array([0.0, 1.0]),
+            lons=np.array([10.0]),  # length 1 vs 2
+            lats=np.array([20.0, 21.0]),
+            depths=np.array([-1.0, -1.0]),
+        )
+
+
+def test_sample_underway_filters_and_flattens(
+    dummy_underway_inst, sample_underway_coords
+):
+    """_sample_underway evaluates active sensors and flattens multi-output tuple kernels."""
+    configs = [
+        DummySensorConfig(SensorType.VELOCITY, enabled=True),
+        DummySensorConfig(SensorType.TEMPERATURE, enabled=True),
+        DummySensorConfig(SensorType.SALINITY, enabled=True),
+    ]
+
+    sampled = dummy_underway_inst._sample_underway(
+        config_sensors=configs,
+        fieldset=None,
+        coords=sample_underway_coords,
+    )
+
+    assert len(sampled) == 4  # total flattened arrays (u, v, temp, sal)
+    np.testing.assert_array_equal(sampled[0], np.array([0.5, 0.6], dtype=np.float32))
+    np.testing.assert_array_equal(sampled[1], np.array([-0.1, -0.2], dtype=np.float32))
+    np.testing.assert_array_equal(sampled[2], np.array([15.0, 16.0], dtype=np.float32))
+    np.testing.assert_array_equal(sampled[3], np.array([35.0, 35.1], dtype=np.float32))
+
+
+def test_to_parquet_writes_valid_file(
+    tmp_path, dummy_underway_inst, sample_underway_coords
+):
+    """_to_parquet writes a valid Parquet table with expected schema metadata and data values."""
+    out_path = tmp_path / "output.parquet"
+    dat_arrays = [
+        np.array([15.0, 16.0], dtype=np.float32),
+        np.array([35.0, 35.1], dtype=np.float32),
+    ]
+    var_names = ["temp", "sal"]
+    origin = np.datetime64("2026-01-01T00:00:00")
+
+    dummy_underway_inst._to_parquet(
+        dat_arrays=dat_arrays,
+        var_names=var_names,
+        fieldset_time_origin=origin,
+        out_path=out_path,
+        coords=sample_underway_coords,
+    )
+
+    assert out_path.exists()
+
+    # verify parquet table, metadata, and columns
+    table = pq.read_table(out_path)
+    schema = table.schema
+
+    assert table.column_names == [
+        "t",
+        "z",
+        "y",
+        "x",
+        "particle_id",
+        "dt",
+        "state",
+        "temp",
+        "sal",
+    ]
+    assert schema.metadata[b"feature_type"] == b"trajectory"
+    assert b"units" in schema.field("t").metadata
+
+    np.testing.assert_array_equal(
+        table["x"].to_numpy(), np.array(sample_underway_coords.lons, dtype=np.float32)
+    )
+    np.testing.assert_array_equal(table["temp"].to_numpy(), dat_arrays[0])
+
+
+def test_parquet_openable_by_read_particles(tmp_path):
+    """Test that a parquet file written by _to_parquet can be read back by parcels.read_particlefile."""
+    coords = UnderwayCoordinates(
+        times=np.array([0.0, 3600.0]),
+        lons=np.array([-5.0, -5.1]),
+        lats=np.array([50.0, 50.1]),
+        depths=np.array([-2.0, -2.0]),
+    )
+    dat_arrays = [
+        np.array([15.0, 16.0], dtype=np.float32),
+        np.array([35.0, 35.1], dtype=np.float32),
+    ]
+    var_names = ["temp", "sal"]
+    origin = np.datetime64("2026-01-01T00:00:00")
+
+    UnderwayInstrument._to_parquet(
+        dat_arrays=dat_arrays,
+        var_names=var_names,
+        fieldset_time_origin=origin,
+        out_path=tmp_path / "test_particles.parquet",
+        coords=coords,
+    )
+
+    # read it back
+    results = parcels.read_particlefile(tmp_path / "test_particles.parquet")
+
+    assert len(results) == 2
+
+    # sampling results may differ slightly due to float32 precision
+    assert np.isclose(results["temp"][0], 15.0)
+    assert np.isclose(results["sal"][1], 35.1)

From 701ad044a31d9a4ce53c263aa019ee25e4d95a9a Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 30 Jul 2026 10:09:55 +0200
Subject: [PATCH 16/25] update ADCP tests

---
 tests/instruments/test_adcp.py | 69 ++++++++++++++++++----------------
 1 file changed, 37 insertions(+), 32 deletions(-)

diff --git a/tests/instruments/test_adcp.py b/tests/instruments/test_adcp.py
index 16604b67..ab13f99c 100644
--- a/tests/instruments/test_adcp.py
+++ b/tests/instruments/test_adcp.py
@@ -4,10 +4,11 @@
 from typing import ClassVar
 
 import numpy as np
+import parcels
+import polars as pl
 import pydantic
 import pytest
 import xarray as xr
-from parcels import FieldSet
 
 from virtualship.instruments.adcp import ADCPInstrument
 from virtualship.instruments.sensors import SensorType
@@ -22,6 +23,7 @@
 BASE_TIME = datetime.datetime.strptime(
     "1950-01-01", "%Y-%m-%d"
 )  # arbitrary time offset for the dummy fieldset
+MIN_DEPTH = -5
 MAX_DEPTH = -1000
 NUM_BINS = 40
 
@@ -53,8 +55,6 @@ class schedule:
 
 
 def test_simulate_adcp(tmpdir, adcp_expedition) -> None:
-    MIN_DEPTH = -5
-
     # where to sample
     sample_points = [
         Spacetime(Location(1, 2), BASE_TIME + datetime.timedelta(seconds=0)),
@@ -93,53 +93,58 @@ def test_simulate_adcp(tmpdir, adcp_expedition) -> None:
     u[1, 0, 1, 1] = expected_obs[1]["U"]["max_depth"]
     u[1, 1, 1, 1] = expected_obs[1]["U"]["surface"]
 
-    fieldset = FieldSet.from_data(
-        {
-            "V": v,
-            "U": u,
+    # make ds
+    times = np.array([expected_obs[0]["time"], expected_obs[1]["time"]])
+    lats = np.array([expected_obs[0]["lat"], expected_obs[1]["lat"]])
+    lons = np.array([expected_obs[0]["lon"], expected_obs[1]["lon"]])
+
+    ds_fields = xr.Dataset(
+        data_vars={
+            "U": (["time", "depth", "lat", "lon"], u, {"units": "m s-1"}),
+            "V": (["time", "depth", "lat", "lon"], v, {"units": "m s-1"}),
         },
-        {
-            "lat": np.array([expected_obs[0]["lat"], expected_obs[1]["lat"]]),
-            "lon": np.array([expected_obs[0]["lon"], expected_obs[1]["lon"]]),
-            "depth": np.array([MAX_DEPTH, MIN_DEPTH]),
-            "time": np.array(
-                [
-                    np.datetime64(expected_obs[0]["time"]),
-                    np.datetime64(expected_obs[1]["time"]),
-                ]
-            ),
+        coords={
+            "time": ("time", times, {"axis": "T"}),
+            "depth": ("depth", [MAX_DEPTH, MIN_DEPTH], {"units": "m", "axis": "Z"}),
+            "lat": ("lat", lats, {"units": "degrees_north"}),
+            "lon": ("lon", lons, {"units": "degrees_east"}),
         },
     )
 
+    # to fieldset
+    fields = {"U": ds_fields["U"], "V": ds_fields["V"]}
+    ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
+    fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
+
     adcp_instrument = ADCPInstrument(adcp_expedition, from_data=None)
-    out_path = tmpdir.join("out.zarr")
+    out_path = tmpdir.join("out.parquet")
 
     adcp_instrument.load_input_data = lambda: fieldset
     adcp_instrument.simulate(sample_points, out_path)
 
-    results = xr.open_zarr(out_path)
+    results = parcels.read_particlefile(out_path)
 
-    # test if output is as expected
-    assert len(results.trajectory) == NUM_BINS
+    assert np.unique(results["z"].to_numpy()).size == NUM_BINS
 
     # for every obs, check if the variables match the expected observations
     # we only verify at the surface and max depth of the adcp, because in between is tricky
-    for traj, vert_loc in [
-        (results.trajectory[0], "max_depth"),
-        (results.trajectory[-1], "surface"),
+    for df_depth, vert_loc in [
+        (results.filter(pl.col("z") == MAX_DEPTH), "max_depth"),
+        (results.filter(pl.col("z") == MIN_DEPTH), "surface"),
     ]:
-        obs_all = results.sel(trajectory=traj).obs
-        assert len(obs_all) == len(sample_points)
-        for i, (obs_i, exp) in enumerate(zip(obs_all, expected_obs, strict=True)):
-            obs = results.sel(trajectory=traj, obs=obs_i)
-            for var in ["lat", "lon"]:
-                obs_value = obs[var].values.item()
-                exp_value = exp[var]
+        assert len(df_depth) == len(sample_points)
+
+        for i, (obs_i, exp) in enumerate(
+            zip(df_depth.iter_rows(named=True), expected_obs, strict=True)
+        ):
+            for var in [("y", "lat"), ("x", "lon")]:
+                obs_value = obs_i[var[0]]
+                exp_value = exp[var[1]]
                 assert np.isclose(obs_value, exp_value), (
                     f"Observation incorrect {vert_loc=} {obs_i=} {var=} {obs_value=} {exp_value=}."
                 )
             for var in ["V", "U"]:
-                obs_value = obs[var].values.item()
+                obs_value = obs_i[var]
                 exp_value = exp[var][vert_loc]
                 assert np.isclose(obs_value, exp_value), (
                     f"Observation incorrect {vert_loc=} {i=} {var=} {obs_value=} {exp_value=}."

From 763636abfdaa736902ed066117ba08e4872a8748 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 30 Jul 2026 10:23:23 +0200
Subject: [PATCH 17/25] tidy up comments

---
 src/virtualship/instruments/base.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 08bdf5d3..e87d1b82 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -54,7 +54,7 @@ class Instrument(abc.ABC):
     sensor_kernels: ClassVar[dict[SensorType, collections.abc.Callable]]
 
     def __init_subclass__(cls, **kwargs: object) -> None:
-        """Ensure non-abstract subclasses (i.e. final instrument classes) define sensor_kernels as a class attribute."""
+        """Ensure non-abstract subclasses (i.e. final/concrete instrument classes) define sensor_kernels as a class attribute."""
         super().__init_subclass__(**kwargs)
         if inspect.isabstract(cls):
             return
@@ -139,7 +139,7 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = False  # TODO: just for dev; remove before merging
+        TMP = True  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
@@ -340,6 +340,7 @@ def _sample_underway(
         ]  # perform sampling
 
         # ensure that sampled is a flat list of arrays, even if some kernels return tuples/lists of arrays
+        # e.g. ADCP kernel returns (u, v) tuple of arrays, whilst UnderwaterST returns single array of temperature/salinity
         sampled_flat = [
             arr
             for item in sampled
@@ -360,7 +361,7 @@ def _to_parquet(
         """
         Write underway instrument data to a Parquet file mirroring the Parcels v4 ParticleFile schema.
 
-        Designed so that output files can be re-read back in with Parcels.read_particlefile for downstream workflows.
+        Designed so that output files can be re-read back in with Parcels.read_particlefile for consistent downstream workflows with non-underway instruments.
         """
         assert len(dat_arrays) == len(var_names), (
             "dat_arrays and var_names must have the same length"
@@ -387,7 +388,7 @@ def _to_parquet(
                 "Conventions": "CF-1.6/CF-1.7",
                 "ncei_template_version": "NCEI_NetCDF_Trajectory_Template_v2.0",
                 "parcels_version": parcels.__version__,
-                "parcels_grid_mesh": "spherical",  # TODO: is the case as long as using Copernicus Marine data...
+                "parcels_grid_mesh": "spherical",
             },
         )
 

From 96d474d4b5120787917046bfe943e1a1df9f440e Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 30 Jul 2026 10:39:51 +0200
Subject: [PATCH 18/25] update UnderwaterST tests

---
 tests/instruments/test_ship_underwater_st.py | 95 ++++++++++----------
 1 file changed, 48 insertions(+), 47 deletions(-)

diff --git a/tests/instruments/test_ship_underwater_st.py b/tests/instruments/test_ship_underwater_st.py
index 016734a1..1bf1689c 100644
--- a/tests/instruments/test_ship_underwater_st.py
+++ b/tests/instruments/test_ship_underwater_st.py
@@ -4,10 +4,10 @@
 from typing import ClassVar
 
 import numpy as np
+import parcels
 import pydantic
 import pytest
 import xarray as xr
-from parcels import FieldSet
 
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.ship_underwater_st import Underwater_STInstrument
@@ -54,87 +54,88 @@ class schedule:
 
 
 def test_simulate_ship_underwater_st(tmpdir, underwater_st_expedition) -> None:
-    # arbitrary time offset for the dummy fieldset
-    base_time = datetime.datetime.strptime("1950-01-01", "%Y-%m-%d")
-
     # where to sample
     sample_points = [
-        Spacetime(Location(1, 2), base_time + datetime.timedelta(seconds=0)),
-        Spacetime(Location(3, 4), base_time + datetime.timedelta(seconds=1)),
+        Spacetime(Location(1, 2), BASE_TIME + datetime.timedelta(seconds=0)),
+        Spacetime(Location(3, 4), BASE_TIME + datetime.timedelta(seconds=1)),
     ]
 
     # expected observations at sample points
     expected_obs = [
         {
-            "salinity": 5,
-            "temperature": 6,
+            "S": 5,
+            "T": 6,
             "lat": sample_points[0].location.lat,
             "lon": sample_points[0].location.lon,
-            "time": base_time + datetime.timedelta(seconds=0),
+            "time": BASE_TIME + datetime.timedelta(seconds=0),
         },
         {
-            "salinity": 7,
-            "temperature": 8,
+            "S": 7,
+            "T": 8,
             "lat": sample_points[1].location.lat,
             "lon": sample_points[1].location.lon,
-            "time": base_time + datetime.timedelta(seconds=1),
+            "time": BASE_TIME + datetime.timedelta(seconds=1),
         },
     ]
 
     # create fieldset based on the expected observations
     # indices are time, latitude, longitude
     salinity = np.zeros((2, 2, 2))
-    salinity[0, 0, 0] = expected_obs[0]["salinity"]
-    salinity[1, 1, 1] = expected_obs[1]["salinity"]
+    salinity[0, 0, 0] = expected_obs[0]["S"]
+    salinity[1, 1, 1] = expected_obs[1]["S"]
 
     temperature = np.zeros((2, 2, 2))
-    temperature[0, 0, 0] = expected_obs[0]["temperature"]
-    temperature[1, 1, 1] = expected_obs[1]["temperature"]
-
-    fieldset = FieldSet.from_data(
-        {
-            "V": np.zeros((2, 2, 2)),
-            "U": np.zeros((2, 2, 2)),
-            "S": salinity,
-            "T": temperature,
+    temperature[0, 0, 0] = expected_obs[0]["T"]
+    temperature[1, 1, 1] = expected_obs[1]["T"]
+
+    # make ds
+    times = np.array([expected_obs[0]["time"], expected_obs[1]["time"]])
+    lats = np.array([expected_obs[0]["lat"], expected_obs[1]["lat"]])
+    lons = np.array([expected_obs[0]["lon"], expected_obs[1]["lon"]])
+
+    ds_fields = xr.Dataset(
+        data_vars={
+            "T": (["time", "lat", "lon"], temperature, {"units": "degC"}),
+            "S": (["time", "lat", "lon"], salinity, {"units": "psu"}),
         },
-        {
-            "lat": np.array([expected_obs[0]["lat"], expected_obs[1]["lat"]]),
-            "lon": np.array([expected_obs[0]["lon"], expected_obs[1]["lon"]]),
-            "time": np.array(
-                [
-                    np.datetime64(expected_obs[0]["time"]),
-                    np.datetime64(expected_obs[1]["time"]),
-                ]
-            ),
+        coords={
+            "time": ("time", times, {"axis": "T"}),
+            "lat": ("lat", lats, {"units": "degrees_north"}),
+            "lon": ("lon", lons, {"units": "degrees_east"}),
         },
     )
 
-    from_data = None
+    # to fieldset
+    fields = {"T": ds_fields["T"], "S": ds_fields["S"]}
+    ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
+    fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
 
-    st_instrument = Underwater_STInstrument(underwater_st_expedition, from_data)
-    out_path = tmpdir.join("out.zarr")
+    st_instrument = Underwater_STInstrument(underwater_st_expedition, from_data=None)
+    out_path = tmpdir.join("out.parquet")
 
     st_instrument.load_input_data = lambda: fieldset
-    # The instrument expects measurements as sample_points
     st_instrument.simulate(sample_points, out_path)
 
-    # test if output is as expected
-    results = xr.open_zarr(out_path)
+    results = parcels.read_particlefile(out_path)
 
-    assert len(results.trajectory) == 1  # expect a single trajectory
-    traj = results.trajectory.item()
-    assert len(results.sel(trajectory=traj).obs) == len(
-        sample_points
-    )  # expect as many obs as sample points
+    # expect a single depth level
+    assert np.unique(results["z"].to_numpy()).size == 1
+
+    # expect as many obs as sample points (given the period is 5 minutes and the sample points are 1 second apart)
+    assert len(results) == len(sample_points)
 
     # for every obs, check if the variables match the expected observations
     for i, (obs_i, exp) in enumerate(
-        zip(results.sel(trajectory=traj).obs, expected_obs, strict=True)
+        zip(results.iter_rows(named=True), expected_obs, strict=True)
     ):
-        obs = results.sel(trajectory=traj, obs=obs_i)
-        for var in ["salinity", "temperature", "lat", "lon"]:
-            obs_value = obs[var].values.item()
+        for var in [("y", "lat"), ("x", "lon")]:
+            obs_value = obs_i[var[0]]
+            exp_value = exp[var[1]]
+            assert np.isclose(obs_value, exp_value), (
+                f"Observation incorrect {obs_i=} {var=} {obs_value=} {exp_value=}."
+            )
+        for var in ["T", "S"]:
+            obs_value = obs_i[var]
             exp_value = exp[var]
             assert np.isclose(obs_value, exp_value), (
                 f"Observation incorrect {i=} {var=} {obs_value=} {exp_value=}."

From 2149ef1f0c45448901d91af02b99bb6814f5d805 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 30 Jul 2026 15:15:12 +0200
Subject: [PATCH 19/25] use fieldset.time_interval.left / right for fieldset
 start / end times

---
 src/virtualship/instruments/adcp.py               | 6 +-----
 src/virtualship/instruments/ctd.py                | 8 ++------
 src/virtualship/instruments/ship_underwater_st.py | 6 +-----
 src/virtualship/instruments/xbt.py                | 8 ++------
 4 files changed, 6 insertions(+), 22 deletions(-)

diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index 4fe0459f..bd253230 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -83,12 +83,8 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # use first active field for time reference
-        _time_ref_key = next(iter(self.variables))
-        _time_ref_field = getattr(fieldset, _time_ref_key)
-        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
-
         # times in seconds since fieldset time origin, expanded across depth bins
+        fieldset_starttime = fieldset.time_interval.left
         times = np.array(
             [
                 (np.datetime64(point.time) - fieldset_starttime)
diff --git a/src/virtualship/instruments/ctd.py b/src/virtualship/instruments/ctd.py
index d6764130..5b8cba61 100644
--- a/src/virtualship/instruments/ctd.py
+++ b/src/virtualship/instruments/ctd.py
@@ -168,12 +168,8 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # use first active field for time reference
-        _time_ref_key = next(iter(self.variables))
-        _time_ref_field = getattr(fieldset, _time_ref_key)
-
-        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
-        fieldset_endtime = _time_ref_field.data.time.isel(time=-1).values
+        fieldset_starttime = fieldset.time_interval.left
+        fieldset_endtime = fieldset.time_interval.right
 
         # deploy time for all ctds should be later than fieldset start time
         if not all(
diff --git a/src/virtualship/instruments/ship_underwater_st.py b/src/virtualship/instruments/ship_underwater_st.py
index 3519c2c8..4aff0af5 100644
--- a/src/virtualship/instruments/ship_underwater_st.py
+++ b/src/virtualship/instruments/ship_underwater_st.py
@@ -77,12 +77,8 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # use first active field for time reference
-        _time_ref_key = next(iter(self.variables))
-        _time_ref_field = getattr(fieldset, _time_ref_key)
-        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
-
         # sampling times and locations
+        fieldset_starttime = fieldset.time_interval.left
         times = np.array(
             [
                 (np.datetime64(point.time) - fieldset_starttime)
diff --git a/src/virtualship/instruments/xbt.py b/src/virtualship/instruments/xbt.py
index 7046fde1..9de23092 100644
--- a/src/virtualship/instruments/xbt.py
+++ b/src/virtualship/instruments/xbt.py
@@ -120,12 +120,8 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # use first active field for time reference
-        _time_ref_key = next(iter(self.variables))
-        _time_ref_field = getattr(fieldset, _time_ref_key)
-
-        fieldset_starttime = _time_ref_field.data.time.isel(time=0).values
-        fieldset_endtime = _time_ref_field.data.time.isel(time=-1).values
+        fieldset_starttime = fieldset.time_interval.left
+        fieldset_endtime = fieldset.time_interval.right
 
         # deploy time for all xbts should be later than fieldset start time
         if not all(

From 2e5b5b0bf502ebf391465b35e70580839c8d6fc8 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Thu, 30 Jul 2026 15:22:04 +0200
Subject: [PATCH 20/25] particle_id is constant for underway instruments

---
 src/virtualship/instruments/base.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index e87d1b82..d5140817 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -139,7 +139,7 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = True  # TODO: just for dev; remove before merging
+        TMP = False  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
@@ -290,7 +290,7 @@ def instrument_type(self) -> InstrumentType:
 
 @dataclass(frozen=True)
 class UnderwayCoordinates:
-    """1D evaluation grid arrays for underway instruments."""
+    """1D sampling location arrays for underway instruments."""
 
     times: np.ndarray  # seconds since origin
     lons: np.ndarray
@@ -413,8 +413,8 @@ def _to_parquet(
                 "y": pa.array(coords.lats.astype(np.float32)),
                 "x": pa.array(coords.lons.astype(np.float32)),
                 "particle_id": pa.array(
-                    np.arange(n, dtype=np.int64)
-                ),  # doesn't necessarily correspond to an actual particle_id in the simulation, but required for Parcels ParticleFile schema
+                    np.zeros(n, dtype=np.int64)
+                ),  # ship is a single 'particle' (here represented by a constant particle_id of 0)
                 "dt": pa.array(np.full(n, np.nan, dtype=np.float64)),
                 "state": pa.array(np.zeros(n, dtype=np.int32)),
                 **{

From 79f0a0c749eb0a2bd67dfdbb1b68c709cd46c2db Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Fri, 31 Jul 2026 10:45:56 +0200
Subject: [PATCH 21/25] dt and state are unnecessary for public facing output

---
 src/virtualship/instruments/base.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index d5140817..4ac7dec7 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -380,8 +380,6 @@ def _to_parquet(
                 pa.field("y", pa.float32()),
                 pa.field("x", pa.float32()),
                 pa.field("particle_id", pa.int64()),
-                pa.field("dt", pa.float64()),
-                pa.field("state", pa.int32()),
             ],
             metadata={
                 "feature_type": "trajectory",

From ccd0ae3c17e847cc43eb9fcc282da5d5fd3c4474 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Fri, 31 Jul 2026 11:13:34 +0200
Subject: [PATCH 22/25] add test for monitoring for schema drift vs. parcels,
 plus some refactoring

---
 tests/instruments/test_base.py | 99 ++++++++++++++++++++++++++++------
 1 file changed, 84 insertions(+), 15 deletions(-)

diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py
index 1d03e1ec..3868f068 100644
--- a/tests/instruments/test_base.py
+++ b/tests/instruments/test_base.py
@@ -336,8 +336,6 @@ def test_to_parquet_writes_valid_file(
         "y",
         "x",
         "particle_id",
-        "dt",
-        "state",
         "temp",
         "sal",
     ]
@@ -350,34 +348,105 @@ def test_to_parquet_writes_valid_file(
     np.testing.assert_array_equal(table["temp"].to_numpy(), dat_arrays[0])
 
 
-def test_parquet_openable_by_read_particles(tmp_path):
-    """Test that a parquet file written by _to_parquet can be read back by parcels.read_particlefile."""
+def _create_underway_parquet(
+    out_path,
+    var_names,
+    dat_arrays=None,
+    origin=np.datetime64("2026-01-01T00:00:00"),  # noqa
+):
+    """Helper to generate an UnderwayInstrument parquet output file."""
     coords = UnderwayCoordinates(
         times=np.array([0.0, 3600.0]),
         lons=np.array([-5.0, -5.1]),
         lats=np.array([50.0, 50.1]),
         depths=np.array([-2.0, -2.0]),
     )
-    dat_arrays = [
-        np.array([15.0, 16.0], dtype=np.float32),
-        np.array([35.0, 35.1], dtype=np.float32),
-    ]
-    var_names = ["temp", "sal"]
-    origin = np.datetime64("2026-01-01T00:00:00")
+
+    if dat_arrays is None:
+        dat_arrays = [
+            np.array([15.0, 16.0], dtype=np.float32),
+            np.array([35.0, 35.1], dtype=np.float32),
+        ]
 
     UnderwayInstrument._to_parquet(
         dat_arrays=dat_arrays,
         var_names=var_names,
         fieldset_time_origin=origin,
-        out_path=tmp_path / "test_particles.parquet",
+        out_path=out_path,
         coords=coords,
     )
 
-    # read it back
-    results = parcels.read_particlefile(tmp_path / "test_particles.parquet")
 
-    assert len(results) == 2
+def dummy_sample_temperature(particles, fieldset):
+    particles.T = fieldset.T[particles.t, particles.z, particles.y, particles.x]
 
-    # sampling results may differ slightly due to float32 precision
+
+def test_parquet_openable_by_read_particles(tmp_path):
+    """Test that a parquet file written by _to_parquet can be read back by parcels.read_particlefile."""
+    parquet_path = tmp_path / "test_particles.parquet"
+    _create_underway_parquet(
+        out_path=parquet_path,
+        var_names=["temp", "sal"],
+        dat_arrays=[
+            np.array([15.0, 16.0], dtype=np.float32),
+            np.array([35.0, 35.1], dtype=np.float32),
+        ],
+    )
+
+    # read back and assert values
+    results = parcels.read_particlefile(parquet_path)
+    assert len(results) == 2
     assert np.isclose(results["temp"][0], 15.0)
     assert np.isclose(results["sal"][1], 35.1)
+
+
+def test_underway_schema_matches_parcels(tmp_path):
+    """Verify that underway instrument parquet output base schema matches Parcels' ParticleFile."""
+    # minimal Parcels FieldSet
+    T = np.zeros((2, 1, 1))
+    T[0, 0, 0], T[1, 0, 0] = 15.0, 16.0
+
+    t1 = np.datetime64("2024-01-01T00:00:00")
+    t2 = np.datetime64("2024-01-02T00:00:00")
+    ds_fields = xr.Dataset(
+        data_vars={"temperature": (["time", "lat", "lon"], T, {"units": "degC"})},
+        coords={
+            "time": ("time", [t1, t2], {"axis": "T"}),
+            "lat": ("lat", [0.0], {"units": "degrees_north"}),
+            "lon": ("lon", [0.0], {"units": "degrees_east"}),
+        },
+    )
+
+    fields = {"T": ds_fields["temperature"]}
+    ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
+    fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
+
+    # parcels simualtion
+    SampleParticle = parcels.Particle.add_variable(parcels.Variable("T"))
+
+    pset = parcels.ParticleSet(
+        fieldset=fieldset, pclass=SampleParticle, t=t1, y=[0.0], x=[0.0]
+    )
+
+    parcels_path = tmp_path / "parcels_particles.parquet"
+    parcels_output = parcels.ParticleFile(parcels_path, outputdt=3600.0)
+    pset.execute(
+        [dummy_sample_temperature],
+        runtime=np.timedelta64(1, "D"),
+        dt=np.timedelta64(60, "m"),
+        output_file=parcels_output,
+    )
+    parcels_df = parcels.read_particlefile(parcels_path)
+
+    # UnderwayInstrument output
+    underway_path = tmp_path / "underway_particles.parquet"
+    _create_underway_parquet(
+        out_path=underway_path,
+        var_names=["T"],
+        dat_arrays=[np.array([15.0, 16.0], dtype=np.float32)],
+        origin=np.datetime64("2024-01-01T00:00:00"),
+    )
+    underway_df = parcels.read_particlefile(underway_path)
+
+    # assert schemas match
+    assert parcels_df.schema == underway_df.schema

From 7af935d63ef752c5f6a1ab619ef3202f927c02d8 Mon Sep 17 00:00:00 2001
From: Jamie Atkins <106238905+j-atkins@users.noreply.github.com>
Date: Fri, 31 Jul 2026 11:21:11 +0200
Subject: [PATCH 23/25] Update tests/instruments/test_base.py

Co-authored-by: Erik van Sebille <erikvansebille@gmail.com>
---
 tests/instruments/test_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py
index 3868f068..44aa5614 100644
--- a/tests/instruments/test_base.py
+++ b/tests/instruments/test_base.py
@@ -381,7 +381,7 @@ def dummy_sample_temperature(particles, fieldset):
     particles.T = fieldset.T[particles.t, particles.z, particles.y, particles.x]
 
 
-def test_parquet_openable_by_read_particles(tmp_path):
+def test_parquet_openable_by_parcels_read_particlefile(tmp_path):
     """Test that a parquet file written by _to_parquet can be read back by parcels.read_particlefile."""
     parquet_path = tmp_path / "test_particles.parquet"
     _create_underway_parquet(

From 9675963f763e71f2cd5567c7aec6334941c35317 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Fri, 31 Jul 2026 11:24:38 +0200
Subject: [PATCH 24/25] remove dev spinner bypass option

---
 src/virtualship/instruments/base.py | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index 4ac7dec7..ce0f8324 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -139,20 +139,17 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = False  # TODO: just for dev; remove before merging
         instrument_name = self.__class__.__name__.split("Instrument")[0]
 
         if not self.verbose_progress:
-            if TMP:
-                with yaspin(
-                    text=f"Simulating {instrument_name} measurements... ",
-                    side="right",
-                    spinner=ship_spinner,
-                ) as spinner:
-                    self.simulate(measurements, out_path)
-                    spinner.ok("✅\n")
-            else:
+            with yaspin(
+                text=f"Simulating {instrument_name} measurements... ",
+                side="right",
+                spinner=ship_spinner,
+            ) as spinner:
                 self.simulate(measurements, out_path)
+                spinner.ok("✅\n")
+
         else:
             print(f"Simulating {instrument_name} measurements... ")
             self.simulate(measurements, out_path)

From eb276a2d2c51d17d80de64fc855ed086b53d2701 Mon Sep 17 00:00:00 2001
From: j-atkins <106238905+j-atkins@users.noreply.github.com>
Date: Fri, 31 Jul 2026 11:27:40 +0200
Subject: [PATCH 25/25] parcels simulation only needs one write step

---
 tests/instruments/test_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py
index 3868f068..268bc031 100644
--- a/tests/instruments/test_base.py
+++ b/tests/instruments/test_base.py
@@ -432,7 +432,7 @@ def test_underway_schema_matches_parcels(tmp_path):
     parcels_output = parcels.ParticleFile(parcels_path, outputdt=3600.0)
     pset.execute(
         [dummy_sample_temperature],
-        runtime=np.timedelta64(1, "D"),
+        runtime=np.timedelta64(60, "m"),
         dt=np.timedelta64(60, "m"),
         output_file=parcels_output,
     )