Trouble with particle times and variable extraction

[amargolskee]

Thanks in advance to anyone who helps with this.

I'm having an issue with reading variables in at specific particle times. I want to evaluate variables for each particle at their initial seeding location at the their initial seeding time. I tried to follow one of the tutorials to do that.

First I read in x,y,z,t particle locations from a csv file. The times for particle seeding iterate from 1 to 67 in intervals of 6. I've checked that the fieldset time (set as fieldset.U.grid.time) produces particle times that are (0...2592000...51840360... and so on), as I expected.

However, when I run my code, I get this error:

INFO: Compiled ArraySampleParticleInitZeroSampleJNO3 ==> /tmp/parcels-347546/libb16df9b1eb30e586550f130776c8a1e9_0.so
WARNING: dt or runtime are zero, or endtime is equal to Particle.time. The kernels will be executed once, without incrementing time
Traceback (most recent call last):
  File "/scratch/gpfs/am9753/Parcels/Parcels_mpi/sose_parcels_d20.py", line 16, in <module>
    mitparc.run_mitgcm("jit")
  File "/scratch/gpfs/am9753/Parcels/Parcels_mpi/mitgcm_parcels_d20.py", line 81, in run_mitgcm
    pset.execute(
  File "/home/am9753/.conda/envs/py3_parcels_mpi3/lib/python3.11/site-packages/parcels/particleset/baseparticleset.py", line 450, in execute
    self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
  File "/home/am9753/.conda/envs/py3_parcels_mpi3/lib/python3.11/site-packages/parcels/kernel/kernelsoa.py", line 227, in execute
    recovery_kernel(p, self.fieldset, p.time)
  File "/home/am9753/.conda/envs/py3_parcels_mpi3/lib/python3.11/site-packages/parcels/tools/statuscodes.py", line 221, in recovery_kernel_time_extrapolation
    raise OutOfTimeError(particle, fieldset)
parcels.tools.statuscodes.OutOfTimeError: 0
Particle P[311370](lon=3.333333, lat=-77.790306, depth=-2.100000, BLGBIONO3=0.000000, BLGMLD=0.000000, time=2592000.000000)
Time: 2013-02-04T00:00:00.000000000,	timestep dt: -0.000000
Field sampled outside time domain at time 2013-02-04T00:00:00.000000000. Try setting allow_time_extrapolation to True

2013-02-04 is definitely well within the runtime for the netcdf file that I'm reading from, so I'm very confused by this error. There are 511 5-day timesteps in the file. Here's a snippet of the netcdf structure. I would have thought that 2013-02-04 should just be 35 days into the model output. Is it necessary to set allow_time_extrapolation to True in this case? And if so, why?

It seems to run if I comment out the pset.execute() calls that have dt=0.. but it would be great to be able to get the initial variable reads to work at each particle's initial location.

netcdf bsose_Iter135_2013to2019_5day_BLGBIONO3 {
dimensions:
	time = 511 ;
	XC = 2160 ;
	YC = 588 ;
	Z = 52 ;
variables:
	int64 iter(time) ;
		iter:standard_name = "timestep" ;
		iter:long_name = "model timestep number" ;
	int64 time(time) ;
		time:standard_name = "time" ;
		time:long_name = "Time" ;
		time:axis = "T" ;
		time:units = "seconds since 2012-12-01" ;
		time:calendar = "proleptic_gregorian" ;

Here is the full file I'm using, in case that's helpful.

from datetime import timedelta as delta
from os import path
import sys
import numpy as np
import xarray as xr
import pandas as pd
import numpy.matlib
from parcels import (
    AdvectionRK4_3D,
    FieldSet,
    JITParticle,
    ParticleFile,
    ParticleSet,
    ScipyParticle,
    Variable,
    ErrorCode,
)

ptype = {"scipy": ScipyParticle, "jit": JITParticle}

class SampleParticleInitZero(JITParticle):         # Define a new particle class
    BLGBIONO3 = Variable('BLGBIONO3', initial=0)  # Variable 'JNO3' initialised by sampling the BLGBION
    BLGMLD = Variable('BLGMLD', initial=0)  # Variable 'MLD' initialised by sampling the BLGMLD


def SampleJNO3(particle, fieldset, time):
    particle.BLGBIONO3 = fieldset.BLGBIONO3[time, particle.depth, particle.lat, particle.lon]

def SampleMLD(particle, fieldset, time):
    particle.BLGMLD = fieldset.BLGMLD[time, 0, particle.lat, particle.lon]

def DeleteParticle(particle, fieldset, time):
    particle.delete()

def run_mitgcm(mode):
    """Function that shows how to load MITgcm data in a zonally periodic domain"""

    data_path = path.join(path.dirname(__file__), "/scratch/gpfs/am9753/SOSE/")
    filenames = {
        "U": data_path + "bsose_Iter135_2013to2019_5day_AllVel.nc",
        "V": data_path + "bsose_Iter135_2013to2019_5day_AllVel.nc",
        "W": data_path + "bsose_Iter135_2013to2019_5day_AllVel.nc",
        "BLGBIONO3": data_path + "bsose_Iter135_2013to2019_5day_BLGBIONO3.nc",
        "BLGMLD": data_path + "bsose_Iter135_2013to2019_5day_BLGMLD.nc"
    }
    variables = {"U": "u", "V": "v","W": "w", "BLGBIONO3": "BLGBIONO3","BLGMLD": "BLGMLD"}
    dimensions = {
        "U": {"lon": "XG", "lat": "YG", "time": "time","depth": "Z"},
        "V": {"lon": "XG", "lat": "YG", "time": "time","depth": "Z"},
        "W": {"lon": "XG", "lat": "YG", "time": "time","depth": "Z"},
        "BLGBIONO3": {"lon": "XC", "lat": "YC", "time": "time","depth": "Z"},
        "BLGMLD": {"lon": "XC", "lat": "YC", "time": "time"}
    }
    fieldset = FieldSet.from_mitgcm(filenames, variables,dimensions,
           allow_time_extrapolation=False)

#    data = pd.read_csv('init_loc.csv',header=None)
    data = pd.read_csv('init_loc_alld20xyt.csv',header=None)

    tmp = np.matrix(data)
    z = tmp[:,0]-1
    y = tmp[:,1]-1
    x = tmp[:,2]-1
    t = tmp[:,3]-1
    
    pset = ParticleSet(
         fieldset=fieldset,
         pclass=SampleParticleInitZero,
         lon = fieldset.U.grid.lon[x],
         lat = fieldset.U.grid.lat[y],
         depth = fieldset.U.grid.depth[z],
         time  = fieldset.U.grid.time[t],
    )

    pfile = ParticleFile(
        "SOSE_d20long_30" + ".zarr", pset, outputdt=delta(days=30)
    )
    sampleJNO3_kernel = pset.Kernel(SampleJNO3)
    sampleMLD_kernel = pset.Kernel(SampleMLD)
    kernels = AdvectionRK4_3D + sampleJNO3_kernel + sampleMLD_kernel
    pset.execute(
        sampleJNO3_kernel, dt=0,
        recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}
    )
    pset.execute(
        sampleMLD_kernel, dt=0,
        recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}
    )

    pset.execute(
        kernels, runtime=delta(days=2160), dt=delta(minutes=30), 
        output_file=pfile,
        recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle} 
    )

    pfile.close()

[erikvansebille]

Thanks for reporting, @amargolskee. It could be that the calendar isn't parsed well when you create the fieldset, do you get any warnings about that? And what if you print fieldset.U.grid.time_full? What dates are then in the fieldset?

[amargolskee]

Thanks for getting back to me. The fieldset.U.grid.time_full and fieldset.U.grid.time dates appear identical.

Both are in units of seconds (0... 432000... up to 432000*511) and seemingly correspond to the amount of time since "2012-12-01" as in the netcdf file. Should these instead have been converted to datetime somewhere?

[erikvansebille]

Yes, they should indeed have been converted to a date time object; that apparently hasn't happened. It's probably beacuse Parcels doesn't know how to deal with the time:calendar = "proleptic_gregorian" ; in your netcdf file...

What is the value of fieldset.U.grid.time_origin?

[amargolskee]

I see. I didn't prepare these files myself: I am using B-SOSE output. The fielsdet.U.grid.time_origin is 2013-01-05T00:00:00.000000000

My python is a little shaky. Is there a quick fix I can use to convert the fieldset.U.grid.time* fields?

[erikvansebille]

I'm not sure that the difference between setting dt=delta(seconds=0.1) and using the default (which should be the same as using dt=delta(seconds=1.0), can you check?) is not related to this issue. Could it be that this value of 45885600.0 is the end of your time series? I.e. is your input data 45885600.0 seconds long?

[hadarberman1]

Hi,

I'm not sure how to check the default dt value? Nevertheless, I'm sure that the default value differs from the specified one, since changing dt changes the runtime.

I have found that for certain dt values there is no failure, i.e. dt=30s, dt=0.09s, and for other dt values the code fails at different time steps with the same time extrapolation error. It seems to me that this issue is related to the ratio of dt and 1 hr (which is the time interval between every two input files).

Regarding your question, 45885600 is not the end of my input, and I am able to access the data with different dt values.

While I can make slight adjustments to dt, I'm still concerned that it could encounter issues at some point due to the 'dt' to one-hour ratio. Does this make sense/ do you have any suggestions?

Thanks for your assistance,
Hadar.

[erikvansebille]

I'm not sure how to check the default dt value? Nevertheless, I'm sure that the default value differs from the specified one, since changing dt changes the runtime.

Well, the default in parcels is dt=1. seconds, see below
https://github.com/OceanParcels/parcels/blob/0265b5e73d475dd2f2bed64b84393872d56a591f/parcels/particleset.py#L803-L804

I'm not sure then what causes your error. Could you try updating to the new Parcels v3.0.0; in which the kernel loop has been rewritten. Perhaps that fixes the issue?

[hadarberman1]

Hi Erik,

Sure I can give it a try, how can I download the 3.0.0 version (using conda would be the best option for me)? Is it compatible with multi processing?

Thanks again,
Hadar.

[erikvansebille]

We released parcels v3 via condo-forge last week. However, some people are experiencing problems with installing (see conda-forge/parcels-feedstock#100), so best way to install it is using (note that it can take a long time for the repo data to be resolved)

conda create -n parcels3 -c conda-forge parcels=3

And then parcels activate parcels3