Temperature sampling kernel returns all zeros in OceanParcels

[hyc006]

Question

Question

Hi everyone,
I’m trying to sample the temperature field along particle trajectories using OceanParcels, following the tutorial example here:
🔗 https://docs.oceanparcels.org/en/latest/examples/tutorial_sampling.html
I added a new variable called "TEMP" to my fieldset and wrote a kernel like this:

def SampleTEMP(particle, fieldset, time):
    particle.temperature = fieldset.TEMP[time, particle.depth, particle.lat, particle.lon]

Then I exported the results to a CSV file. The output file correctly contains a column named "temperature", but all of the values are 0, suggesting that no temperature data was actually sampled.
I’m wondering what might be going wrong here. Could it be:
that the "TEMP" field is not being read or interpolated correctly in the kernel?
Here’s how I defined my fieldset (simplified):

filenames = {'U': fname, 'V': fname, 'W': fname, 'TEMP': fname}
variables = {'U': 'U', 'V': 'V', 'W': 'W', 'TEMP': 'TEMP'}
dimensions = {'U': {'lat': 'lat', 'lon': 'lon', 'depth': 'depth', 'time': 'time'},
              'V': {'lat': 'lat', 'lon': 'lon', 'depth': 'depth', 'time': 'time'},
              'W': {'lat': 'lat', 'lon': 'lon', 'depth': 'depth', 'time': 'time'},
              'TEMP': {'lat': 'lat', 'lon': 'lon', 'depth': 'depth', 'time': 'time'}}

I’d really appreciate any suggestions or common pitfalls I might be missing — for instance, whether TEMP needs a specific interpolation method, or if sampling order (time vs depth vs lat/lon) matters.
Thank you very much for your time and help!

[erikvansebille]

Hi @hyc006, thanks for reaching out!

I'm not sure what is the issue here. I don't see anything wrong with the two code snippets you posted above. Could you also provide the rest of the code?

And are you sure that your fieldset.TEMP has values? Can you do a Field interpolation outside of the Kernel and/or pset.execute()? For example, what is the output of

print(pset[0].time, pset[0].depth, pset[0].lat, pset[0].lon)
print(fieldset.TEMP[pset[0].time, pset[0].depth, pset[0].lat, pset[0].lon]

[hyc006]

Thank you so much for your reply. I’ve attached my code for running the parcel backtracking while trying to sample temperature:
back_track_sample.py

When I check the field using fieldset, I can see the temperature field as follows:
<Field> name : 'TEMP' grid : RectilinearZGrid(lon=array([ 53.68, 53.71, 53.73, ..., 57.14, 57.16, 57.19], shape=(350,), dtype=float32), lat=array([-5.99, -5.96, -5.93, ..., -3.41, -3.39, -3.36], shape=(281,), dtype=float32), time=array([ 0.00, 3600.00, 7200.00, ..., 5641200.00, 5644800.00, 5648400.00], shape=(1570,)), time_origin=1900-01-01T00:00:00.000000000, mesh='spherical') extrapolate time: False time_periodic : False gridindexingtype: 'nemo' to_write : False
However, when I try the code you suggested:
print(fieldset.TEMP[pset[0].time, pset[0].depth, pset[0].lat, pset[0].lon]
I get the following error:
AttributeError: 'DeferredArray' object has no attribute 'shape'
I think this means I did something wrong in my setup, but I’m not quite sure where. I’m still new to OceanParcels, so I really appreciate your help and guidance!

[erikvansebille]

Hi @hyc006, I had a look at your python script and I think the solution is quite simple. You defined a SampleTEMP Kernel, but didn't include it in the list of Kernels in pset.execute(). Hence, the temperature was never sampled and the output contained the initial particle.temperature values of 0

So simply adding SampleTEMP to combo might fix the issue

[hyc006]

Hi Erik,
Thanks a lot for your help, and sorry for this kind of basic question!
One quick follow-up — I tried adding SampleTEMP in the combo kernel like this:
combo = adv_kernel + keep_kernel + leave_kernel + SampleTEMP # run advection first, then sample temperature
But I now get a new error that didn’t appear before:

INFO: Output files are stored in parcels_back_sample.zarr.
63%|████████████████████████▋              | 547200.0/864000.0 [02:04<01:12, 4344.06it/s]
Traceback (most recent call last):
  ...
parcels.tools.statuscodes.FieldOutOfBoundSurfaceError: Field sampled out-of-bound at the surface, at (depth=-2.8e-10, lat=-4.69, lon=55.84)

I’m wondering if this could be because the parcels reach the surface, but my temperature field doesn’t have values defined exactly at the surface or bottom? My ROMS grid stores temperature at the midpoints of the vertical cells, so there are no data points on the boundaries.
Would you suggest a good way to handle this (e.g., clipping depth or interpolating)?
Thanks again for your time and help!

[michaeldenes]

To chime in, @hyc006, a possible quick fix is to determine what the 'surface' of your temperature field is, and then sample only on or below it. For example, something like this:

def SampleTEMP(particle, fieldset, time):
    surface_depth = -2.8e-10 # Change this to your surface
    if particle.depth < surface_depth: # The sign here depends on if your +ve down, or +ve up, check this too.
        particle.temperature = fieldset.TEMP[time, surface_depth, particle.lat, particle.lon]
    else:
        particle.temperature = fieldset.TEMP[time, particle.depth, particle.lat, particle.lon]