Parcels-code · erikvansebille · May 4, 2026 · Apr 17, 2026 · Apr 17, 2026 · Apr 17, 2026
diff --git a/docs/development/policies.md b/docs/development/policies.md
@@ -22,7 +22,7 @@ Parcels follows [Intended Effort Versioning (EffVer)](https://jacobtomlinson.dev
 
 When making backward incompatible changes, we will make sure these changes and instructions to upgrade are communicated to the user via change logs or migration guides, and (where applicable) informative error messaging.
 
-Note when conducting research we highly recommend documenting which version of Parcels (and other packages) you are using. This can be as easy as doing `conda env export > environment.yml` alongside your project code. The Parcels version used to generate an output file is also stored as metadata entry in the `.zarr` output file.
+Note when conducting research we highly recommend documenting which version of Parcels (and other packages) you are using. This can be as easy as doing `conda env export > environment.yml` alongside your project code. The Parcels version used to generate an output file is also stored as metadata entry in the `.parquet` output file.
 
 ## Changes in policies
 

diff --git a/docs/getting_started/explanation_concepts.md b/docs/getting_started/explanation_concepts.md
@@ -186,7 +186,7 @@ pset.execute(kernels=kernels, dt=dt, runtime=runtime)
 
 ### Output
 
-To analyse the particle data generated in the simulation, we need to define a `parcels.ParticleFile` and add it as an argument to `parcels.ParticleSet.execute()`. The output will be written in a [zarr format](https://zarr.readthedocs.io/en/stable/), which can be opened as an `xarray.Dataset`. The dataset will contain the particle data with at least `time`, `z`, `lat` and `lon`, for each particle at timesteps defined by the `outputdt` argument.
+To analyse the particle data generated in the simulation, we need to define a `parcels.ParticleFile` and add it as an argument to `parcels.ParticleSet.execute()`. The output will be written in a [parquet format](https://parquet.apache.org/), which can be opened as a `polars.DataFrame`. The dataset will contain the particle data with at least `time`, `z`, `lat` and `lon`, for each particle at timesteps defined by the `outputdt` argument.
 
 There are many ways to analyze particle output, and although we provide [a short tutorial to get started](./tutorial_output.ipynb), we recommend writing your own analysis code and checking out [related Lagrangian analysis projects in our community page](../community/index.md#analysis-code).
 

diff --git a/docs/getting_started/tutorial_output.ipynb b/docs/getting_started/tutorial_output.ipynb
diff --git a/docs/getting_started/tutorial_quickstart.md b/docs/getting_started/tutorial_quickstart.md
@@ -13,12 +13,12 @@ read more, we have a [concepts overview](./explanation_concepts.md) discussing t
 
 ## Imports
 
-Parcels depends on `xarray`, expecting inputs in the form of [`xarray.Dataset`](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html)
-and writing output files that can be read with xarray.
+Parcels depends on `xarray`, expecting inputs in the form of [`xarray.Dataset`](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html). Output files can be read with `polars`.
 
 ```{code-cell}
 import numpy as np
 import xarray as xr
+import polars as pl
 import parcels
 import parcels.tutorial
 ```
@@ -123,11 +123,11 @@ Before starting the simulation, we must define where and how frequent we want to
 We can define this in a {py:obj}`parcels.ParticleFile` object:
 
 ```{code-cell}
-output_file = parcels.ParticleFile("output-quickstart.zarr", outputdt=np.timedelta64(1, "h"))
+output_file = parcels.ParticleFile("output-quickstart.parquet", outputdt=np.timedelta64(1, "h"))
 ```
 
-The output files are in `.zarr` [format](https://zarr.readthedocs.io/en/stable/), which can be read by `xarray`.
-See the [Parcels output tutorial](./tutorial_output.ipynb) for more information on the zarr format. We want to choose
+The output files are in `.parquet` [format](https://parquet.apache.org/), which can be read by [Polars](https://pola.rs/).
+See the [Parcels output tutorial](./tutorial_output.ipynb) for more information on the parquet format. We want to choose
 the `outputdt` argument so that it captures the smallest timescales of our interest.
 
 ## Run Simulation: `ParticleSet.execute()`
@@ -152,26 +152,25 @@ pset.execute(
 
 ## Read output
 
-To start analyzing the trajectories computed by **Parcels**, we can open the `ParticleFile` using `xarray`:
+To start analyzing the trajectories computed by **Parcels**, we can open the `ParticleFile` using the `read_particlefile()` utility, which itself uses `polars`:
 
 ```{code-cell}
-ds_particles = xr.open_zarr("output-quickstart.zarr")
-ds_particles
+df = parcels.read_particlefile("output-quickstart.parquet")
+df
 ```
 
-The 10 particle trajectories are stored along the `trajectory` dimension, and each trajectory contains 25 observations
-(initial values + 24 hourly timesteps) along the `obs` dimension. The [working with Parcels output tutorial](./tutorial_output.ipynb)
-provides more detail about the dataset and how to analyse it.
+The file contains 250 rows: 25 observations for the 10 particle trajectories.
+The [working with Parcels output tutorial](./tutorial_output.ipynb) provides more detail about the dataset and how to analyse it.
 
 Let's verify that Parcels has computed the advection of the virtual particles!
 
 ```{code-cell}
 import matplotlib.pyplot as plt
 
-# plot positions and color particles by number of observation
-scatter = plt.scatter(ds_particles.lon.T, ds_particles.lat.T, c=np.repeat(ds_particles.obs.values,npart))
-plt.scatter(ds_particles.lon[:,0],ds_particles.lat[:,0],facecolors="none",edgecolors='r') # starting positions
-plt.scatter(lon,lat,facecolors="none",edgecolors='r') # starting positions
+# plot positions and color particles by time
+scatter = plt.scatter(df['lon'], df['lat'], c=df['time'])
+plt.scatter(df['lon'][:npart], df['lat'][:npart], facecolors="none", edgecolors='r') # starting positions
+plt.scatter(lon, lat, facecolors="none", edgecolors='r') # starting positions
 plt.xlim(31,33)
 plt.ylabel("Latitude [deg N]")
 plt.ylim(-33,-30)
@@ -196,7 +195,7 @@ location!
 ```{code-cell}
 :tags: [hide-output]
 # set up output file
-output_file = parcels.ParticleFile("output-backwards.zarr", outputdt=np.timedelta64(1, "h"))
+output_file = parcels.ParticleFile("output-backwards.parquet", outputdt=np.timedelta64(1, "h"))
 
 # execute simulation in backwards time
 pset.execute(
@@ -210,10 +209,11 @@ pset.execute(
 When we check the output, we can see that the particles have returned to their original position!
 
 ```{code-cell}
-ds_particles_back = xr.open_zarr("output-backwards.zarr")
+df_back = parcels.read_particlefile("output-backwards.parquet")
 
-scatter = plt.scatter(ds_particles_back.lon.T, ds_particles_back.lat.T, c=np.repeat(ds_particles_back.obs.values,npart))
-plt.scatter(ds_particles_back.lon[:,0],ds_particles_back.lat[:,0],facecolors="none",edgecolors='r') # starting positions
+scatter = plt.scatter(df_back['lon'], df_back['lat'], c=df_back['time'])
+particles_at_max_time = df_back.filter(pl.col("time") == df_back["time"].max())
+plt.scatter(particles_at_max_time['lon'], particles_at_max_time['lat'], facecolors="none", edgecolors='r') # starting positions
 plt.xlabel("Longitude [deg E]")
 plt.xlim(31,33)
 plt.ylabel("Latitude [deg N]")
@@ -226,6 +226,7 @@ Using Euler forward advection, the final positions are equal to the original pos
 
 ```{code-cell}
 # testing that final location == original location
-np.testing.assert_almost_equal(ds_particles_back['lat'].values[:,-1],ds_particles['lat'].values[:,0], 2)
-np.testing.assert_almost_equal(ds_particles_back['lon'].values[:,-1],ds_particles['lon'].values[:,0], 2)
+particles_at_min_time = df_back.filter(pl.col("time") == df_back["time"].min())
+np.testing.assert_almost_equal(particles_at_min_time["lat"], lat, 2)
+np.testing.assert_almost_equal(particles_at_min_time['lon'], lon, 2)
 ```
diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md
@@ -108,7 +108,7 @@ lats = np.linspace(-32.5, -30.5, npart)
 
 pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, lat=lats, lon=lons)
 output_file = parcels.ParticleFile(
-    store="advection_then_wind.zarr", outputdt=np.timedelta64(6,'h')
+    path="advection_then_wind.parquet", outputdt=np.timedelta64(6,'h')
 )
 pset.execute(
     [parcels.kernels.AdvectionRK2, wind_kernel],
@@ -126,7 +126,7 @@ pset_reverse = parcels.ParticleSet(
     fieldset, pclass=parcels.Particle, z=z, lat=lats,  lon=lons
 )
 output_file_reverse = parcels.ParticleFile(
-    store="wind_then_advection.zarr", outputdt=np.timedelta64(6,"h")
+    path="wind_then_advection.parquet", outputdt=np.timedelta64(6,"h")
 )
 pset_reverse.execute(
     [wind_kernel, parcels.kernels.AdvectionRK2],
@@ -140,10 +140,14 @@ Finally, plot the trajectories to show that they are identical in the two simula
 
 ```{code-cell}
 # Plot the resulting particle trajectories overlapped for both cases
-advection_then_wind = xr.open_zarr("advection_then_wind.zarr")
-wind_then_advection = xr.open_zarr("wind_then_advection.zarr")
-plt.plot(wind_then_advection.lon.T, wind_then_advection.lat.T, "-")
-plt.plot(advection_then_wind.lon.T, advection_then_wind.lat.T, "--", c="k", alpha=0.7)
+advection_then_wind = parcels.read_particlefile("advection_then_wind.parquet")
+wind_then_advection = parcels.read_particlefile("wind_then_advection.parquet")
+
+fig, ax = plt.subplots(figsize=(5, 3))
+for traj in wind_then_advection.partition_by("particle_id", maintain_order=True):
+    ax.plot(traj["lon"], traj["lat"], "-")
+for traj in advection_then_wind.partition_by("particle_id", maintain_order=True):
+    ax.plot(traj["lon"], traj["lat"], "--", c="k", alpha=0.7)
 plt.show()
 ```
 

diff --git a/docs/user_guide/examples/tutorial_Argofloats.ipynb b/docs/user_guide/examples/tutorial_Argofloats.ipynb
@@ -159,9 +159,8 @@
     "\n",
     "# Create a ParticleFile object to store the output\n",
     "output_file = parcels.ParticleFile(\n",
-    "    store=\"argo_float.zarr\",\n",
+    "    \"argo_float.parquet\",\n",
     "    outputdt=timedelta(minutes=15),\n",
-    "    chunks=(1, 500),  # setting to write in chunks of 500 observations\n",
     ")\n",
     "\n",
     "# Now execute the Kernels for 30 days, saving data every 30 minutes\n",
@@ -183,20 +182,6 @@
     "First plot the depth as a function of time, with the temperature as color (only on the upcast)."
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "ds_particles = xr.open_zarr(output_file.store)\n",
-    "x = ds_particles[\"lon\"][:].squeeze()\n",
-    "y = ds_particles[\"lat\"][:].squeeze()\n",
-    "z = ds_particles[\"z\"][:].squeeze()\n",
-    "time = ds_particles[\"time\"][:].squeeze()\n",
-    "temp = ds_particles[\"temp\"][:].squeeze()"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -205,10 +190,12 @@
    "source": [
     "import matplotlib.pyplot as plt\n",
     "\n",
+    "df = parcels.read_particlefile(\"argo_float.parquet\")\n",
+    "\n",
     "fig = plt.figure(figsize=(13, 6))\n",
     "ax = plt.axes()\n",
-    "ax.plot(time, z, color=\"gray\")\n",
-    "cb = ax.scatter(time, z, c=temp, s=20, marker=\"o\", zorder=2)\n",
+    "ax.plot(df[\"time\"], df[\"z\"], color=\"gray\")\n",
+    "cb = ax.scatter(df[\"time\"], df[\"z\"], c=df[\"temp\"], s=20, marker=\"o\", zorder=2)\n",
     "ax.set_xlabel(\"Time [days]\")\n",
     "ax.set_ylabel(\"Depth (m)\")\n",
     "ax.invert_yaxis()\n",
@@ -234,12 +221,12 @@
     "fig = plt.figure(figsize=(13, 8))\n",
     "ax = plt.axes(projection=\"3d\")\n",
     "ax.view_init(azim=-145)\n",
-    "ax.plot3D(x, y, z, color=\"gray\")\n",
-    "cb = ax.scatter(x, y, z, c=temp, s=20, marker=\"o\", zorder=2)\n",
+    "ax.plot3D(df[\"lon\"], df[\"lat\"], df[\"z\"], color=\"gray\")\n",
+    "cb = ax.scatter(df[\"lon\"], df[\"lat\"], df[\"z\"], c=df[\"temp\"], s=20, marker=\"o\", zorder=2)\n",
     "ax.set_xlabel(\"Longitude\")\n",
     "ax.set_ylabel(\"Latitude\")\n",
     "ax.set_zlabel(\"Depth (m)\")\n",
-    "ax.set_zlim(np.max(z), 0)\n",
+    "ax.set_zlim(df[\"z\"].max(), 0)\n",
     "fig.colorbar(cb, label=\"Temperature (°C)\")\n",
     "plt.show()"
    ]
@@ -261,7 +248,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.14.2"
+   "version": "3.14.4"
   }
  },
  "nbformat": 4,

diff --git a/docs/user_guide/examples/tutorial_croco_3D.ipynb b/docs/user_guide/examples/tutorial_croco_3D.ipynb
@@ -34,7 +34,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt\nimport numpy as np\nimport xarray as xr\n\nimport parcels\nimport parcels.tutorial\n\nds_fields = parcels.tutorial.open_dataset(\"CROCOidealized_data/data\")\n\nds_fields.load();  # Preload data to speed up access"
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import polars as pl\n",
+    "import xarray as xr\n",
+    "\n",
+    "import parcels\n",
+    "import parcels.tutorial\n",
+    "\n",
+    "ds_fields = parcels.tutorial.open_dataset(\"CROCOidealized_data/data\")\n",
+    "\n",
+    "ds_fields.load();  # Preload data to speed up access"
    ]
   },
   {
@@ -119,7 +129,7 @@
     ")\n",
     "\n",
     "outputfile = parcels.ParticleFile(\n",
-    "    store=\"croco_particles3D.zarr\",\n",
+    "    path=\"croco_particles3D.parquet\",\n",
     "    outputdt=np.timedelta64(5000, \"s\"),\n",
     ")\n",
     "\n",
@@ -149,10 +159,11 @@
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
-    "ds = xr.open_zarr(\"croco_particles3D.zarr\")\n",
+    "df = pl.read_parquet(\"croco_particles3D.parquet\")\n",
     "\n",
     "ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
-    "ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
+    "for traj in df.partition_by(\"particle_id\", maintain_order=True):\n",
+    "    ax.plot(traj[\"lon\"] / 1e3, traj[\"z\"], \".-\")\n",
     "\n",
     "for z in ds_fields.s_w.values:\n",
     "    ax.plot(\n",
@@ -208,7 +219,7 @@
     ")\n",
     "\n",
     "outputfile = parcels.ParticleFile(\n",
-    "    store=\"croco_particles_noW.zarr\", outputdt=np.timedelta64(5000, \"s\")\n",
+    "    path=\"croco_particles_noW.parquet\", outputdt=np.timedelta64(5000, \"s\")\n",
     ")\n",
     "\n",
     "pset_noW.execute(\n",
@@ -219,10 +230,11 @@
     ")\n",
     "\n",
     "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
-    "ds = xr.open_zarr(\"croco_particles_noW.zarr\")\n",
+    "df = pl.read_parquet(\"croco_particles_noW.parquet\")\n",
     "\n",
     "ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
-    "ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
+    "for traj in df.partition_by(\"particle_id\", maintain_order=True):\n",
+    "    ax.plot(traj[\"lon\"] / 1e3, traj[\"z\"], \".-\")\n",
     "\n",
     "for z in ds_fields.s_w.values:\n",
     "    ax.plot(\n",
@@ -306,7 +318,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.14.2"
+   "version": "3.14.4"
   }
  },
  "nbformat": 4,