Output time column truncated to seconds, causing duplicate timestamps at milli-second output frequency

[ChitYanToe]

Parcels version

Parcels v4

Description

When running a Parcels simulation with an output frequency (outputdt) on the order of milliseconds, the time variable written to the output file does not preserve the millisecond precision. Instead, time values appear to be truncated/rounded to the nearest second.
As a result, multiple consecutive output records for the same trajectory (Particle_ID) end up sharing an identical time value, even though they correspond to distinct simulation steps. This creates duplicate (trajectory, time) pairs in the output dataset, which breaks downstream analysis that assumes each (Particle_ID, time) pair is unique (e.g., indexing, grouping, or merging on time).

Code sample

output_file = parcels.ParticleFile(resultDir+"output-quickstart.parquet", outputdt=timedelta(seconds=0.4), mode="w") #

[erikvansebille]

Hmm, I looked into this put can't immediately reproduce the issue. See the unit test below, which does pass as expected

@pytest.mark.parametrize("dt", [100, 200])
def test_subsecond_outputdt(fieldset, dt, tmp_parquet):
    """Test that outputdt can be subsecond and that the output times are correct."""

    def Update_lon(particles, fieldset):  # pragma: no cover
        particles.dx += dt / 1000.0  # Move at a rate of 1 unit per second

    pset = ParticleSet(fieldset, x=[0], y=[0])
    ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(dt, "ms"))
    pset.execute(Update_lon, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(dt, "ms"), output_file=ofile)

    df = pd.read_parquet(tmp_parquet)
    np.testing.assert_allclose(df["x"].values, np.arange(0, 1+1e-6, dt / 1000.), atol=1e-6)
    np.testing.assert_allclose(df["t"].values, np.arange(0, 1+1e-6, dt / 1000.), atol=1e-6)

@ChitYanToe can you check what is fundamentally different between this unit test and your code? For example, what happens if you use outputdt=np.timedelta64(400, "ms")?

[erikvansebille]

OK, a follow-up. I think (part of) the issue is with parcels.read_particlefile(). If I change the test above to use that method to read the output, it fails

def test_subsecond_outputdt(fieldset, dt, tmp_parquet):
    """Test that outputdt can be subsecond and that the output times are correct."""

    def Update_lon(particles, fieldset):  # pragma: no cover
        particles.dx += dt / 1000.0  # Move at a rate of 1 unit per second

    pset = ParticleSet(fieldset, x=[0], y=[0])
    ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(dt, "ms"))
    pset.execute(Update_lon, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(dt, "ms"), output_file=ofile)

    df = parcels.read_particlefile(tmp_parquet)

    np.testing.assert_allclose(df["x"], np.arange(0, 1+1e-6, dt / 1000.), atol=1e-6)
    expected_t = np.arange(0, 1001, dt).astype("timedelta64[ms]")
    elapsed_t = (df["t"] - df["t"].min()).to_numpy().astype("timedelta64[ms]")
    np.testing.assert_allclose(elapsed_t, expected_t, atol=1)

If I however make one change to readparticlefile, it passes

--- a/src/parcels/_core/particlefile.py
+++ b/src/parcels/_core/particlefile.py
@@ -275,7 +275,7 @@ def read_particlefile(path: PathLike, decode_times: bool = True) -> pd.DataFrame
 
     values = table.column("t").to_numpy()
     var = xr.Variable(("t",), values, attrs)
-    values = xr.coders.CFDatetimeCoder(time_unit="s").decode(var).values
+    values = xr.coders.CFDatetimeCoder(time_unit="ns").decode(var).values
     if "since" in attrs["units"]:
         values = values.astype("datetime64[ns]")
         df = df.with_columns(pl.Series("t", values, dtype=pl.Datetime("ns")))

[ChitYanToe]

Thanks a lot.

I tried to use pd.read_parquet from pandas, and I re-formatted the time t column as mm:ss.sss for instance. Then, further analysis worked.

Hmm, I looked into this put can't immediately reproduce the issue. See the unit test below, which does pass as expected

@pytest.mark.parametrize("dt", [100, 200])
def test_subsecond_outputdt(fieldset, dt, tmp_parquet):
"""Test that outputdt can be subsecond and that the output times are correct."""

def Update_lon(particles, fieldset):  # pragma: no cover
    particles.dx += dt / 1000.0  # Move at a rate of 1 unit per second

pset = ParticleSet(fieldset, x=[0], y=[0])
ofile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(dt, "ms"))
pset.execute(Update_lon, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(dt, "ms"), output_file=ofile)

df = pd.read_parquet(tmp_parquet)
np.testing.assert_allclose(df["x"].values, np.arange(0, 1+1e-6, dt / 1000.), atol=1e-6)
np.testing.assert_allclose(df["t"].values, np.arange(0, 1+1e-6, dt / 1000.), atol=1e-6)

@ChitYanToe can you check what is fundamentally different between this unit test and your code? For example, what happens if you use outputdt=np.timedelta64(400, "ms")?