Update underway

src/virtualship/instruments/adcp.py

@@ -1,48 +1,38 @@
 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
+from virtualship.utils import register_instrument
 
 # =====================================================
-# SECTION: Dataclass
+# SECTION: Kernels
 # =====================================================
 
 
-@dataclass
-class ADCP:
-    """ADCP configuration."""
-
-    name: ClassVar[str] = "ADCP"
-
-
-# =====================================================
-# SECTION: non-sensor Particle Variables (non-sampling)
-# =====================================================
-
-# ADCP has no non-sensor variables, only sensor variables.
-_ADCP_NONSENSOR_VARIABLES: list = []
+# 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: parcels.FieldSet, coords: UnderwayCoordinates):
+    # eval
+    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(coords.lats))
+    v = v * 1852 * 60
 
-def _sample_velocity(particles, fieldset):
-    particles.U, particles.V = fieldset.UV.eval(
-        particles.t,
-        particles.z,
-        particles.x,
-        particles.y,
-        applyConversion=False,
-    )
+    return u, v
 
 
 # =====================================================
@@ -51,11 +41,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):
@@ -74,9 +64,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(
@@ -87,53 +77,42 @@ 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)
 
         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
+        # 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)
+                / 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])
         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),
-            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
+        # 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)),
+        )
 
-        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)
-            )
+        sampled = self._sample_underway(
+            config_sensors=adcp_config.sensors, fieldset=fieldset, coords=coords
+        )
 
-            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,
+        )

src/virtualship/instruments/base.py

@@ -2,15 +2,19 @@
 
 import abc
 import collections
+import inspect
 import tempfile
 from dataclasses import dataclass
 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
 
@@ -50,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/concrete 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."
@@ -132,20 +139,17 @@ def simulate(
 
     def execute(self, measurements: list, out_path: str | Path) -> None:
         """Run instrument simulation."""
-        TMP = True  # 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)
@@ -279,3 +283,141 @@ 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(frozen=True)
+class UnderwayCoordinates:
+    """1D sampling location arrays for underway instruments."""
+
+    times: np.ndarray  # seconds since origin
+    lons: np.ndarray
+    lats: np.ndarray
+    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):
+    """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
+        # 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
+            for arr in (item if isinstance(item, (tuple, list)) else (item,))
+        ]
+
+        return sampled_flat
+
+    @staticmethod
+    def _to_parquet(
+        dat_arrays: list[np.ndarray],
+        var_names: list[str],
+        fieldset_time_origin: np.datetime64,
+        out_path: Path | str,
+        coords: UnderwayCoordinates,
+        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 consistent downstream workflows with non-underway instruments.
+        """
+        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()),
+            ],
+            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",
+            },
+        )
+
+        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}"
+            )
+
+        # build table with all data, including sampled variables
+        table = pa.table(
+            {
+                "t": pa.array(coords.times.astype(np.float64)),
+                "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(
+                    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)),
+                **{
+                    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)

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(