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 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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],