Fieldset date does not sync with the netcdf downloaded

[sohinidudhat-dotcom]

Question

Question

Hi, I have downloaded netcdfs for currents, winds and stokes drifts for the year 2024. I imported them into parcels and proceeded to build a fieldset. However, the fieldset date-times are showing up as ' FieldSet time coverage:
Min: 1970-01-01 00:00:00+00:00
Max: 1970-01-02 00:00:00+00:00' and hence I am not able to advect my particle set which is released at varying dates in the year 2024 'Particle release times:
DatetimeIndex(['2024-01-10 12:00:00+00:00', '2024-05-10 12:00:00+00:00',
'2024-07-10 12:00:00+00:00', '2024-09-10 12:00:00+00:00',
'2024-12-10 12:00:00+00:00'],
dtype='datetime64[ns, UTC]', freq=None)'. I get an error when I use pset.execute ' TimeExtrapolationError: U sampled outside time domain at time 2078-12-10T12:00:00.000000000.'

Supporting code/error messages

# Paste your code within this block
```import math
from datetime import timedelta
from operator import attrgetter

import matplotlib.pyplot as plt
import numpy as np
import trajan as ta
import xarray as xr
from IPython.display import HTML
from matplotlib.animation import FuncAnimation

from parcels import (
    AdvectionRK4,
    FieldSet,
    JITParticle,
    ParticleSet,
    Variable,
    download_example_dataset,
)
import netCDF4 as nc
import numpy as np

# STEP 1: importing all the netcdfs and naming variables


import os

cur24 =  #my netcdf file path
if not os.path.isfile(cur24):
    raise FileNotFoundError(f"File not found: {cur24}")
# Define the mapping for the variables
variables = {'U': 'uo', 'V': 'vo'}

# Define the dimensions in the NetCDF file
dimensions = {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}

# Create the filenames dictionary
filenames = {'U': cur24, 'V': cur24}

wind24 = #my netcdf file path

stokes24= #my netcdf file path

# Variables and dimensions
wind_vars = {'U': 'u10', 'V': 'v10'}
stokes_vars = {'U': 'ust', 'V': 'vst'}

dims = {'lon': 'longitude', 'lat': 'latitude', 'time': 'valid_time'}

# Load FieldSets
wind_fs = FieldSet.from_netcdf(wind24, wind_vars, dims, mesh='flat')
stokes_fs = FieldSet.from_netcdf(stokes24, stokes_vars, dims, mesh='flat')

from parcels import FieldSet


# Wind
wind_vars = {'U': 'u10', 'V': 'v10'}
wind_dims = {'lon': 'longitude', 'lat': 'latitude', 'time': 'valid_time'}
wind_fieldset = FieldSet.from_netcdf(wind24, wind_vars, wind_dims, mesh='flat')

# Stokes
stokes_vars = {'U': 'ust', 'V': 'vst'}
stokes_dims = {'lon': 'longitude', 'lat': 'latitude', 'time': 'valid_time'}
stokes_fieldset = FieldSet.from_netcdf(stokes24, stokes_vars, stokes_dims, mesh='flat')


# STEP 2: in this step I tried to fix the date issue before proceeding to fieldset computing, but this did not really work. 


import xarray as xr

def fix_time_encoding(ds, timevar="time"):
    if timevar in ds:
        encoded = xr.coding.times.encode_cf_datetime(
            ds[timevar].values,
            units="seconds since 1970-01-01 00:00:00",
            calendar="gregorian"
        )
        # Some versions of xarray return (values, units, calendar)
        # Extract just the numeric values
        if isinstance(encoded, tuple):
            numeric_time = encoded[0]
        else:
            numeric_time = encoded

        # Replace coordinate values directly
        ds = ds.assign_coords({timevar: numeric_time.astype("int64")})
        
        # Add CF metadata
        ds[timevar].attrs["units"] = "seconds since 1970-01-01 00:00:00"
        ds[timevar].attrs["calendar"] = "gregorian"
    return ds

# Open and fix each dataset
ds_cur = xr.open_dataset(cur24)
ds_wind = xr.open_dataset(wind24)
ds_stokes = xr.open_dataset(stokes24)

ds_cur = fix_time_encoding(ds_cur, "time")
ds_wind = fix_time_encoding(ds_wind, "valid_time")
ds_stokes = fix_time_encoding(ds_stokes, "valid_time")

# Save back to temporary NetCDFs for Parcels to use
ds_cur.to_netcdf("cur24_fixed.nc")
ds_wind.to_netcdf("wind24_fixed.nc")
ds_stokes.to_netcdf("stokes24_fixed.nc")


# STEP 3: fieldset 

filenames = {'U': 'cur24_fixed.nc', 'V': 'cur24_fixed.nc'}
variables = {'U': 'uo', 'V': 'vo'}
dimensions = {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}

fieldset = FieldSet.from_netcdf(filenames=filenames, variables=variables, dimensions=dimensions, time_origin=pd.Timestamp("2024-01-01"))

fieldset.computeTimeChunk()


# checking if this worked 
import pandas as pd

print("FieldSet time coverage:")
print("Min:", pd.to_datetime(fieldset.U.grid.time[0], unit="s", utc=True))
print("Max:", pd.to_datetime(fieldset.U.grid.time[-1], unit="s", utc=True))

### OUTPUT:  FieldSet time coverage:
#Min: 1970-01-01 00:00:00+00:00
#Max: 1970-01-02 00:00:00+00:00



#STEP4 : setting particles set
# releasing five particles in different months

from parcels import JITParticle, ParticleSet, Variable
from datetime import datetime, timezone
import numpy as np
import pandas as pd


# Custom particle with an 'age' variable
class AgeParticle(JITParticle):
    age = Variable("age", initial=0)
    # Add wind and Stokes components as particle variables
    Uw = Variable("Uw", initial=0.0)
    Vw = Variable("Vw", initial=0.0)
    Us = Variable("Us", initial=0.0)
    Vs = Variable("Vs", initial=0.0)

# Initial positions (lon, lat)
initial_positions = [
    (58.939909, 20.488748),
    (58.822802, 20.382851),
    (58.779116, 20.265725),
    (58.623517, 20.246082),
    (58.702360, 20.432772)
]

# Assign a different release date to each particle
release_dates = [
    datetime(2024, 1, 10, 12, 0, tzinfo=timezone.utc),
    datetime(2024, 5, 10, 12, 0, tzinfo=timezone.utc),
    datetime(2024, 7, 10, 12, 0, tzinfo=timezone.utc),
    datetime(2024, 9, 10, 12, 0, tzinfo=timezone.utc),
    datetime(2024, 12, 10, 12, 0, tzinfo=timezone.utc),
]

# Convert to seconds since 1970 in a robust way
release_times = pd.to_datetime(release_dates, utc=True).astype("int64") // 10**9
release_times = release_times.astype(float).tolist()

print(release_times)
# Split lons and lats
lons, lats = zip(*initial_positions)

# Create ParticleSet with one particle per (lon, lat, time)
pset = ParticleSet(
    fieldset=fieldset,
    pclass=AgeParticle,
    lon=lons,
    lat=lats,
    time=release_times
)

# STEP 5:  I defined some kernels and proceeded to the pset execute

from datetime import timedelta

output_file = pset.ParticleFile(
    name="EddyParticles.zarr",  # the file name
    outputdt=timedelta(hours=1),  # the time step of the outputs
)

# Execute the simulation with your new kernel function (move_particles)
pset.execute(
    kernels,  
    runtime=timedelta(days=3),  # set to 3 days for the simulation
    dt=-timedelta(minutes=5),  # - TIMEDELTA FOR HINDCAST (backward advection)
    output_file=output_file,
)




```text
# Paste your error message within this block
```INFO: Output files are stored in EddyParticles.zarr.
  0%|                                                                                     | 0/259200.0 [00:00<?, ?it/s]
---------------------------------------------------------------------------
TimeExtrapolationError                    Traceback (most recent call last)
Cell In[49], line 9
      3 output_file = pset.ParticleFile(
      4     name="EddyParticles.zarr",  # the file name
      5     outputdt=timedelta(hours=1),  # the time step of the outputs
      6 )
      8 # Execute the simulation with your new kernel function (move_particles)
----> 9 pset.execute(
     10     kernels,  
     11     runtime=timedelta(days=3),  # set to 3 days for the simulation
     12     dt=-timedelta(minutes=5),  # - TIMEDELTA FOR HINDCAST (backward advection)
     13     output_file=output_file,
     14 )

File ~\AppData\Local\anaconda3\Lib\site-packages\parcels\particleset.py:1169, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt, delete_cfiles)
   1165     skip_kernel = True if max(self.time) < (time + dt) else False
   1167 time_at_startofloop = time
-> 1169 next_input = self.fieldset.computeTimeChunk(time, dt)
   1171 # Define next_time (the timestamp when the execution needs to be handed back to python)
   1172 if dt > 0:

File ~\AppData\Local\anaconda3\Lib\site-packages\parcels\fieldset.py:1553, in FieldSet.computeTimeChunk(self, time, dt)
   1551     continue
   1552 if f.grid._update_status == "not_updated":
-> 1553     nextTime_loc = f.grid._computeTimeChunk(f, time, signdt)
   1554     if time == nextTime_loc and signdt != 0:
   1555         raise TimeExtrapolationError(time, field=f)

File ~\AppData\Local\anaconda3\Lib\site-packages\parcels\grid.py:376, in Grid._computeTimeChunk(self, f, time, signdt)
    374 if self._ti == -1:
    375     self.time = self.time_full
--> 376     self._ti, _ = f._time_index(time)
    377     periods = self.periods.value if isinstance(self.periods, c_int) else self.periods
    378     if (
    379         signdt == -1
    380         and self._ti == 0
    381         and (time - periods * (self.time_full[-1] - self.time_full[0])) == self.time[0]
    382         and f.time_periodic
    383     ):

File ~\AppData\Local\anaconda3\Lib\site-packages\parcels\field.py:1102, in Field._time_index(self, time)
   1092 """Find the index in the time array associated with a given time.
   1093 
   1094 Note that we normalize to either the first or the last index
   1095 if the sampled value is outside the time value range.
   1096 """
   1097 if (
   1098     not self.time_periodic
   1099     and not self.allow_time_extrapolation
   1100     and (time < self.grid.time[0] or time > self.grid.time[-1])
   1101 ):
-> 1102     raise TimeExtrapolationError(time, field=self)
   1103 time_index = self.grid.time <= time
   1104 if self.time_periodic:

TimeExtrapolationError: U sampled outside time domain at time 2078-12-10T12:00:00.000000000. Try setting allow_time_extrapolation to True.

[michaeldenes]

Hi @sohinidudhat-dotcom , given your error throws 2078-12-10T12:00:00.000000000, which is approx. 54 years after 2024, and 2024 is 54 years after 1970, my speculation is that, because you fix your datasets to have units of "seconds since 1970-01-01 00:00:00", and then you set the time origin in parcels with time_origin=pd.Timestamp("2024-01-01"), it is reading the date units as "seconds since 2024-01-01 00:00:00".

What happens if you remove the time_origin=pd.Timestamp("2024-01-01") argument when creating your fieldset?

[sohinidudhat-dotcom]

Hi Michael, I removed the time_origin=pd.Timestamp("2024-01-01"), but it still shows the same error.

…

On Tue, Sep 16, 2025 at 1:22 PM Michael Denes ***@***.***> wrote: Hi @sohinidudhat-dotcom <https://github.com/sohinidudhat-dotcom> , given your error throws 2078-12-10T12:00:00.000000000, which is approx. 54 years after 2024, and 2024 is 54 years after 1970, my speculation is that, because you fix your datasets to have units of "seconds since 1970-01-01 00:00:00", and then you set the time origin in parcels with time_origin=pd.Timestamp("2024-01-01"), it is reading the date units as "seconds since 1970-01-01 00:00:00". What happens if you remove the time_origin=pd.Timestamp("2024-01-01") argument when creating your fieldset? — Reply to this email directly, view it on GitHub <#2219 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BXOKKJPBKAPGYVZO4I3VS3L3S76H3AVCNFSM6AAAAACGULHRQOVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINBRGQ3TQOI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Sohini Dudhat PhD student CRMG, Scottish Oceans Institute University of St Andrews, Scotland ***@***.*** https://crmg.st-andrews.ac.uk/

[michaeldenes]

Hi @sohinidudhat-dotcom, what was the purpose of the fix_time_encoding function? What does the data you downloaded look like, is it not in regular units of seconds since 1970-01-01? And does your function actually fix the time dimension (i.e. cur24_fixed.nc has the corrected time values)?

I think the issue is with your data, and that you should explore that avenue first. With that said, it's hard to determine the root cause of the issue without seeing the data.

[sohinidudhat-dotcom]

Hi Michael, Well, the issue is really how the fieldset is reading the netCDF. I do not think the way I tried before fixed it. I went through the TimeStamps and calendars — Parcels Documentation <https://docs.oceanparcels.org/en/stable/examples/tutorial_timestamps.html>, and added timestamps. While the particles are advected by doing this, more than one particle trajectory is not stored in the .zarr file and I am unable to plot it. Best Sohini

…

On Tue, Sep 16, 2025 at 1:47 PM Michael Denes ***@***.***> wrote: Hi @sohinidudhat-dotcom <https://github.com/sohinidudhat-dotcom>, what was the purpose of the fix_time_encoding function? What does the data you downloaded look like, is it not in regular units of seconds since 1970-01-01? And does your function actually fix the time dimension (i.e. cur24_fixed.nc has the corrected time values)? I think the issue is with your data, and that you should explore that avenue first. With that said, it's hard to determine the root cause of the issue without seeing the data. — Reply to this email directly, view it on GitHub <#2219 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BXOKKJLKKNY5DZAGV7WO6AT3TABGPAVCNFSM6AAAAACGULHRQOVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINBRGUZDEOI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Sohini Dudhat PhD student CRMG, Scottish Oceans Institute University of St Andrews, Scotland ***@***.*** https://crmg.st-andrews.ac.uk/

[michaeldenes]

Hi @sohinidudhat-dotcom, without the updated code, or data, or particlefile, it's hard to determine what is going on. For example, we don't know what kernels you are running, and perhaps you are using a deleteParticle style kernel to handle errors, which may be the cause for what may seem like an empty trajectory file (but would actually be expected behaviour).

If you can come up with a very simple example that reproduces the problem then we can help, otherwise we can only speculate as to what the problem is.

[sohinidudhat-dotcom]

I understand, here is the updates, clean code:

Loading libraries

import math
from datetime import timedelta
from operator import attrgetter

import matplotlib.pyplot as plt
import numpy as np
import trajan as ta
import xarray as xr
from IPython.display import HTML
from matplotlib.animation import FuncAnimation

from parcels import (
AdvectionRK4,
FieldSet,
JITParticle,
ParticleSet,
Variable,
download_example_dataset,
)
import netCDF4 as nc
import numpy as np

import cartopy.crs as ccrs
import cartopy.feature as cfeature

Currents daily means netCDF for 2024 imported as velocity field

filename = "cur.nc"

# Define the mapping for the variables
variables = {'U': 'uo', 'V': 'vo'}

Define the dimensions in the NetCDF file

dimensions = {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}

Create the filenames dictionary

filenames = {'U': filename, 'V': filename}

Creating fieldset

fieldset = FieldSet.from_netcdf(
filenames,
variables=variables,
dimensions=dimensions)

Three particles written for different locations and different start time in the year 2024

from parcels import JITParticle, ParticleSet
from datetime import datetime, timedelta

Example particle positions and release times

initial_positions = [(58.939909, 20.488748),
(58.822802, 20.382851),
(58.779116, 20.265725)]

release_dates = [
datetime(2024, 1, 10, 12, 0),
datetime(2024, 5, 10, 12, 0),
datetime(2024, 7, 10, 12, 0)
]

lons, lats = zip(*initial_positions)

pset = ParticleSet(
fieldset=fieldset,
pclass=JITParticle,
lon=lons, lat=lats,
time=[d.replace(tzinfo=None) for d in release_dates] # naive datetime
)

Advection of the particles

from parcels import AdvectionRK4, ParticleFile

output_file = pset.ParticleFile(
name="EddyParticles.zarr", # the file name
outputdt=timedelta(hours=1), # the time step of the outputs
)

pset.execute(
AdvectionRK4,
runtime=timedelta(days=3), # or whatever runtime you want
dt=-timedelta(minutes=5), # Hindcasting the particles for 3 days
output_file=output_file
)
INFO: Output files are stored in EddyParticles.zarr.
100%|███████████████████████████████████████████████████████████████████| 259200.0/259200.0 [00:05<00:00, 48250.54it/s]

the advection was successful but plotting was not

#checking the fieldset and the particles

how many particles are there in the zarr file

num_particles_in_file = ds.sizes['trajectory']
print(f"Particles in zarr file: {num_particles_in_file}")
Output: Particles in zarr file: 1 # output says 1, and not the 3 which are written.

#fieldset coverage
import pandas as pd
print("FieldSet time coverage:")
print("Min:", pd.to_datetime(fieldset.U.grid.time[0], unit="s", utc=True))
print("Max:", pd.to_datetime(fieldset.U.grid.time[-1], unit="s", utc=True))

Output: FieldSet time coverage:
Min: 1970-07-08 00:00:00+00:00
Max: 1970-07-09 00:00:00+00:00 # This is the issue! this coverage is not in 2024

checking the particle date:

import xarray as xr
import pandas as pd

Open the output file created during pset.execute()

ds = xr.open_zarr("EddyParticles.zarr")

The time variable is stored under the "time" dimension

particle_times = ds["time"].values # numpy array of times

If already decoded into datetime64

if "datetime64" in str(particle_times.dtype):
start = particle_times.min()
end = particle_times.max()
else:
# If stored as seconds since epoch
start = pd.to_datetime(particle_times.min(), unit="s")
end = pd.to_datetime(particle_times.max(), unit="s")

print(f"Particles simulated from {start} to {end}")
output: Particles simulated from 2024-07-07T13:00:00.000000000 to 2024-07-10T12:00:00.000000000

[michaeldenes]

As you say

Output: FieldSet time coverage:
Min: 1970-07-08 00:00:00+00:00
Max: 1970-07-09 00:00:00+00:00 # This is the issue! this coverage is not in 2024

then the issue is with your data. I tried your fix_time_encoding on a dataset of my own, but it did not work as expected on my dataset, so perhaps your data correction is not working on your dataset? Can you check this?

In an ideal world, I would have the data provider correct the data to fit the native model output. If you can't do that, then indeed, correctly editing the time stamps should work.

While debugging this, I recommend that you iteratively build up to the simulation you want to run so that it's easier to find where the bug occurs. For example, a simple forward advection at the initial model time (i.e. not providing the time argument in the ParticleSet) will help determine if a simulation runs with output that makes sense. Then adding in backwards advection, and then lastly adding in particular release times.