From cb0366d7078bb48c8f2f854c526f40e6c7ad7654 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 1 May 2026 15:00:13 +0200 Subject: [PATCH 1/6] Add failing unit test to check writing of initial sampling values --- tests/test_particlefile.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 759ac5274..aed00832e 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -303,6 +303,24 @@ def IncreaseAge(particles, fieldset): # pragma: no cover np.testing.assert_equal(df_traj["age"].astype("timedelta64[s]").values, (df_traj["time"] - release_time).values) +def test_sampling_initial_value(fieldset, tmp_parquet): + # Test that inital value of a field gets sampled + + SampleParticle = get_default_particle(np.float64).add_variable(Variable("sample", initial=np.nan)) + + def SampleKernel(particles, fieldset): # pragma: no cover + particles.sample = fieldset.U[particles] + + pset = ParticleSet(fieldset, pclass=SampleParticle, lon=[0], lat=[0]) + ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(1, "s")) + + pset.execute(SampleKernel, runtime=np.timedelta64(2, "s"), dt=np.timedelta64(1, "s"), output_file=ofile) + + df = parcels.read_particlefile(tmp_parquet) + assert np.isfinite(df["sample"][1:]).all() # should not be nan + assert np.isfinite(df["sample"][0]) # should not be nan + + def test_reset_dt(fieldset, tmp_parquet): # Assert that p.dt gets reset when a write_time is not a multiple of dt # for p.dt=0.02 to reach outputdt=0.05 and endtime=0.1, the steps should be [0.2, 0.2, 0.1, 0.2, 0.2, 0.1], resulting in 6 kernel executions From b87eb89febe1fec9b3034e41021cdfd8637fd276 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 28 Jul 2026 16:01:54 +0200 Subject: [PATCH 2/6] Fixing unit test --- tests/test_particlefile.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 531eb4d08..c9790efb7 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -312,14 +312,14 @@ def test_sampling_initial_value(fieldset, tmp_parquet): def SampleKernel(particles, fieldset): # pragma: no cover particles.sample = fieldset.U[particles] - pset = ParticleSet(fieldset, pclass=SampleParticle, lon=[0], lat=[0]) - ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(1, "s")) + pset = ParticleSet(fieldset, pclass=SampleParticle, x=[0], y=[0]) + pset.sample = fieldset.U[0, pset.z, pset.y, pset.x] # Sample initial value + ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(1, "s")) pset.execute(SampleKernel, runtime=np.timedelta64(2, "s"), dt=np.timedelta64(1, "s"), output_file=ofile) df = parcels.read_particlefile(tmp_parquet) - assert np.isfinite(df["sample"][1:]).all() # should not be nan - assert np.isfinite(df["sample"][0]) # should not be nan + np.testing.assert_allclose(df["sample"].is_finite().all(), True) def test_reset_dt(fieldset, tmp_parquet): From 2e5c3e1c2e58c8987f4d558fb77a8b1da147c9e7 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 28 Jul 2026 16:13:52 +0200 Subject: [PATCH 3/6] Also support direct sampling of ParticleSet --- src/parcels/_core/field.py | 5 +++-- tests/test_particlefile.py | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index 56a6917a0..22f9bdc1f 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -8,6 +8,7 @@ import numpy as np from parcels._core.index_search import GRID_SEARCH_ERROR, LEFT_OUT_OF_BOUNDS, RIGHT_OUT_OF_BOUNDS, _search_time_index +from parcels._core.particleset import ParticleSet from parcels._core.particlesetview import ParticleSetView from parcels._core.statuscodes import ( AllParcelsErrorCodes, @@ -186,7 +187,7 @@ def eval(self, t: datetime, z, y, x, particles=None): def __getitem__(self, key): self._check_velocitysampling() try: - if isinstance(key, ParticleSetView): + if isinstance(key, (ParticleSetView, ParticleSet)): return self.eval(key.t, key.z, key.y, key.x, key) else: return self.eval(*key) @@ -295,7 +296,7 @@ def eval(self, t: datetime, z, y, x, particles=None): def __getitem__(self, key): try: - if isinstance(key, ParticleSetView): + if isinstance(key, (ParticleSetView, ParticleSet)): return self.eval(key.t, key.z, key.y, key.x, key) else: return self.eval(*key) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index c9790efb7..cad1d38b8 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -312,8 +312,8 @@ def test_sampling_initial_value(fieldset, tmp_parquet): def SampleKernel(particles, fieldset): # pragma: no cover particles.sample = fieldset.U[particles] - pset = ParticleSet(fieldset, pclass=SampleParticle, x=[0], y=[0]) - pset.sample = fieldset.U[0, pset.z, pset.y, pset.x] # Sample initial value + pset = ParticleSet(fieldset, pclass=SampleParticle, x=[0], y=[0], t=[fieldset.time_interval.left]) + pset.sample = fieldset.U[pset] # Sample initial value ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(1, "s")) pset.execute(SampleKernel, runtime=np.timedelta64(2, "s"), dt=np.timedelta64(1, "s"), output_file=ofile) From c8398f0e964eea54a29e1a6cd9138519fc0b83fd Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 28 Jul 2026 16:44:05 +0200 Subject: [PATCH 4/6] Also testing for multiple particles --- tests/test_particlefile.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index cad1d38b8..b8a24f14b 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -304,7 +304,8 @@ def IncreaseAge(particles, fieldset): # pragma: no cover assert (df_traj["age"] == traj_time).all() -def test_sampling_initial_value(fieldset, tmp_parquet): +@pytest.mark.parametrize("npart", [1, 10]) +def test_sampling_initial_value(fieldset, npart, tmp_parquet): # Test that inital value of a field gets sampled SampleParticle = get_default_particle(np.float64).add_variable(Variable("sample", initial=np.nan)) @@ -312,7 +313,11 @@ def test_sampling_initial_value(fieldset, tmp_parquet): def SampleKernel(particles, fieldset): # pragma: no cover particles.sample = fieldset.U[particles] - pset = ParticleSet(fieldset, pclass=SampleParticle, x=[0], y=[0], t=[fieldset.time_interval.left]) + x = np.zeros(npart) + y = np.zeros(npart) + t = np.zeros(npart, dtype="timedelta64[s]") + + pset = ParticleSet(fieldset, pclass=SampleParticle, x=x, y=y, t=t) pset.sample = fieldset.U[pset] # Sample initial value ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(1, "s")) From 40ad7b666b73b454b3f6e6ed30dc992d6b07efe1 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 28 Jul 2026 16:44:27 +0200 Subject: [PATCH 5/6] Add initial value sampling to tutorial --- .../examples/tutorial_sampling.ipynb | 76 ++++++++++++++++++- 1 file changed, 74 insertions(+), 2 deletions(-) diff --git a/docs/user_guide/examples/tutorial_sampling.ipynb b/docs/user_guide/examples/tutorial_sampling.ipynb index 98336db2e..b83271fed 100644 --- a/docs/user_guide/examples/tutorial_sampling.ipynb +++ b/docs/user_guide/examples/tutorial_sampling.ipynb @@ -11,8 +11,8 @@ "In this tutorial we will go through how particles can sample `Fields`, using temperature as an example. Along the way we will get to know the parcels class `Variable` (see [here](https://parcels.readthedocs.io/en/latest/reference/particles.html#parcels.particle.Variable) for the documentation) and some of its methods. This tutorial covers several applications of a sampling setup:\n", "\n", "- [**Basic along trajectory sampling**](#basic-sampling)\n", - "- [**Sampling velocity fields**](#sampling-velocity-fields)\n", - "- [**Sampling initial conditions**](#sampling-initial-values)" + "- [**Sampling initial Field values**](#sampling-initial-field-values)\n", + "- [**Sampling velocity fields**](#sampling-velocity-fields)" ] }, { @@ -200,6 +200,78 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Sampling initial Field values\n", + "\n", + "When creating a ParticleSet, the initial values of `Variables` are set as a constant (default `0`, but in the example above set to `np.nan`). Those will also be the values that appear as first entry in the ParticleFile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\n", + " [\n", + " traj[\"temperature\"][0]\n", + " for traj in df.partition_by(\"particle_id\", maintain_order=True)\n", + " ]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you want to sample the initial value of a field, you can do so by calling the field with the ParticleSet as an argument. For example, in, the code above you would add one line (`pset.temperature = fieldset.thetao[pset]`) to sample the initial temperature values, beofe calling `pset.execute(...)`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-output" + ] + }, + "outputs": [], + "source": [ + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset, pclass=SampleParticle, x=lon, y=lat, z=z, t=time\n", + ")\n", + "pset.temperature = fieldset.thetao[pset] # <-- sample initial temperature values\n", + "\n", + "output_file = parcels.ParticleFile(\n", + " \"sampletemp_withinit.parquet\", outputdt=timedelta(hours=1)\n", + ")\n", + "\n", + "pset.execute(\n", + " [parcels.kernels.AdvectionRK2, SampleT],\n", + " runtime=timedelta(hours=30),\n", + " dt=timedelta(minutes=5),\n", + " output_file=output_file,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = parcels.read_particlefile(\"sampletemp_withinit.parquet\")\n", + "print(\n", + " [\n", + " round(traj[\"temperature\"][0], 2)\n", + " for traj in df.partition_by(\"particle_id\", maintain_order=True)\n", + " ]\n", + ")" + ] + }, { "attachments": {}, "cell_type": "markdown", From a441f543baa03b5994243fcfe0de154476465d9d Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 28 Jul 2026 17:00:12 +0200 Subject: [PATCH 6/6] Adding initial sampling info to quickstart and concepts tutorials --- docs/user_guide/getting_started/explanation_concepts.md | 1 + docs/user_guide/getting_started/tutorial_quickstart.md | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/docs/user_guide/getting_started/explanation_concepts.md b/docs/user_guide/getting_started/explanation_concepts.md index 8f0f27114..107e583a7 100644 --- a/docs/user_guide/getting_started/explanation_concepts.md +++ b/docs/user_guide/getting_started/explanation_concepts.md @@ -94,6 +94,7 @@ pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, t=t, z=z, ```{admonition} 🖥️ Learn more about how to create ParticleSets :class: seealso - [Release particles at different times](../examples/tutorial_delaystart.ipynb) +- [Sampling initial Field values](../examples/tutorial_sampling.ipynb#sampling-initial-field-values) ``` ## 3. Kernels diff --git a/docs/user_guide/getting_started/tutorial_quickstart.md b/docs/user_guide/getting_started/tutorial_quickstart.md index 508ce6044..c647173a4 100644 --- a/docs/user_guide/getting_started/tutorial_quickstart.md +++ b/docs/user_guide/getting_started/tutorial_quickstart.md @@ -112,6 +112,14 @@ ax = temperature.axes ax.scatter(lon, lat, s=40, c='w', edgecolors='r'); ``` +If you also want to sample the initial value of a field, you can do so by calling the field with the `ParticleSet` as an argument. For example, in the code above you would add one line (assuming that the ParticleSet has a variable `temperature`) to sample the initial temperature values, before calling `pset.execute(...)`. + +```python +pset.temperature = fieldset.thetao[pset] +``` + +See the [sampling tutorial](../examples/tutorial_sampling.ipynb#sampling-initial-field-values) for more information on initial sampling of fields. + ## Compute: `Kernel` After setting up the input data and particle start locations and times, we need to specify what calculations to do with