Issue writing a particle trajectory netCDF file using Parcels v2.1.1

[Lachlan00]

Hi,

I'm having issues writing a trajectory file in Parcels version 2.1.1. In previous version I would generate my file using the following code.

outputfn = 'trajectoryfile.nc'
sampledt = timedelta(days=1)
runtime = timedelta(days=4)

pset = ParticleSet(fieldset=fieldset, pclass=oceancc_particle, lon=lons, lat=lats, time=time)
# set ageing kernel
kernels = pset.Kernel(ageingParticle) + pset.Kernel(stuckParticle) + pset.Kernel(killSwitch) + pset.Kernel(AdvectionRK4)
# collect data
pset.execute(kernels,
             runtime=runtime,  # runtime controls the interval of the plots
             dt=sampledt,
             recovery={ErrorCode.ErrorOutOfBounds: deleteParticle},
             output_file=pset.ParticleFile(name=outputfn, outputdt=sampledt))

Now I would expect this to produce a trajectory file with 25 particles (there are 5 released every day) with 25 trajectories and 5 observations. This is what I got using the previous releases of Ocean Parcels. However what I have gotten is a trajectory file with only 5 trajectories and just 1 observation. I have also tried the following in line with the new tutorials but this produces the same issue...

pset = ParticleSet(fieldset=fieldset, pclass=oceancc_particle, lon=lons, lat=lats, time=time)
# set ageing kernel
kernels = pset.Kernel(ageingParticle) + pset.Kernel(stuckParticle) + pset.Kernel(killSwitch) + pset.Kernel(AdvectionRK4)

# collect data
output_file = pset.ParticleFile(name=outputfn, outputdt=sampledt)
pset.execute(kernels,
             runtime=runtime,
             dt=sampledt,
             recovery={ErrorCode.ErrorOutOfBounds: deleteParticle},
             output_file=output_file)
output_file.export()

Here's the netcdf output when using v2.1.1...

  dimensions:
    obs = 1;
    traj = 5;
  variables:
    int trajectory(traj=5, obs=1);
      :_FillValue = -2147483648; // int
      :long_name = "Unique identifier for each particle";
      :cf_role = "trajectory_id";

And the netCDF output with version 2.0.0...

  dimensions:
    obs = UNLIMITED;   // (5 currently)
    traj = UNLIMITED;   // (25 currently)
  variables:
    int trajectory(traj=25, obs=5);
      :cf_role = "trajectory_id";
      :long_name = "Unique identifier for each particle";
      :_ChunkSizes = 25U, 1U; // uint

Here's my custom kernels and particle classes.

####################
# Custom Particles #
####################
# add particle age
class oceancc_particle(JITParticle):
    # add age of particle
    age = Variable('age', dtype=np.int32, initial=0.)
    stuck = Variable('stuck', dtype=np.int32, initial=0.)
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.

####################
# Particle Kernels #
####################
# Delete particles out of bounds
def deleteParticle(particle, fieldset, time):
    particle.delete()

# Particle ageing kernel
def ageingParticle(particle, fieldset, time):
    particle.age += 1

# Stuck paerticle kernel
def stuckParticle(particle, fieldset, time):
    if (particle.prev_lon == particle.lon) and (particle.prev_lat == particle.lat):
        particle.stuck += 1
    # Set the stored values for next iteration.
    particle.prev_lon = particle.lon
    particle.prev_lat = particle.lat

# Delete old particles
def killSwitch(particle, fieldset, time):
    if particle.age >= 365:
        particle.delete()
    elif particle.stuck >= 20:
        particle.delete()

Thanks!

[erikvansebille]

Thanks for reporting, @Lachlan00. I'd be happy to have a look what's going wrong there. Could you provide (a subset of) the field that you've used to run this? And also the full code? That will help me to debug what's going on. You can email me e.g. a wetransfer link to e.vansebille@uu.nl

[Lachlan00]

Hey no worries, and thanks a lot for this Erik. I just forwarded a data subset to you via WeTransfer. Full code can be found here (I just did another push): https://github.com/Lachlan00/oceancc-lagrangian

Specifically you would run python control_script.py. In its current state it pulls the configuration data from config.py (should be straight forward to interpret). Then only one function is run from lagrangian.py called particle_positions_filterV. Just for context, the script is generating particles in a particular region in the model (specified in config.py) filtered by the vertically averaged southward velocity of the current (in this case only seeding when V < -0.4).

[erikvansebille]

I haven't been able to run your simulation, because I don't have the start locations of your particles.

Can you check how many .npy files there are in the output, at the end? Perhaps at convert_at_end = False when you create your ParticleFile object? And what is the output of print(output_file.__dict__)?

How different are the files between v2.0.0 and v2.1.1 actually? It might also be that the difference is not in the particlewriting, but in e.g. outofbounds particles being stuck?

[Lachlan00]

Hi Erik,

Sorry if the code wasn't clear. The particles start positions are generated within the function particle_positions_filterV() by another function called particle_generator_region_filterV() that's called within it. Basically it just first scans the fieldset in the netCDF and then only 'seeds' particles within the specified "generation region" on cells where V < -0.4.

The only thing you should have to change is the ROMS_dir(the directory where the file I sent you is - important: this should be the only netCDF in the directory, I originally wrote the script for a dataset which was divided into multiple netCDF files for each month so if you multiple in the directory it will attempt to stitch them all together...) and filenames (which is the filepath to the file) variable in the config.py script.

Just in case you have other errors I've emailed you the two output trajectory files. In these file I have not changed the code at all but only swapped between using parcels v2.0.0 and v2.1.1.

I don't think it's an outofbounds error. Here's an animation of the particle simulation the code I sent you should run ( I have increased the run length from 4 to 29 for the animation). Also I have a kernel to delete out of bounds particles and particles that have been stuck for > 20 days. Animation was produced using v2.1.1 with the same setup as for the code producing the trajectory files (see the function particle_animation_filterV(). The very first time step of the particles is not shown in the animation (time of seeding). The black dotted box is the generation region.

[Lachlan00]

If it still doesn't work, here's some particle starting values you can use.

lons = [153.66922791, 153.72326675, 153.88181741, 153.95149217,
       153.76835555, 153.68243381, 153.74658414, 153.91465868,
       153.72307588, 153.56464782, 153.8139013 , 153.7138997 ,
       153.63494884, 153.6315857 , 153.73463102, 153.64318576,
       153.91497795, 153.92575691, 153.94976721, 153.76389011,
       153.78685182, 153.73275929, 153.83047259, 153.59427803,
       153.61823107]

lats = [-29.46804062, -28.91361114, -28.44006511, -27.98574743,
       -27.35424904, -28.94825326, -28.30141886, -27.73546792,
       -27.19637444, -27.04952806, -28.75192198, -27.76621771,
       -27.21560279, -26.88013201, -26.3402223 , -28.84043456,
       -28.35536709, -28.12050136, -28.12820825, -27.63917067,
       -29.02913672, -27.58140788, -27.46977463, -27.25020175,
       -27.2579518 ]

time = [     0.,      0.,      0.,      0.,      0.,  86400.,  86400.,
        86400.,  86400.,  86400., 172800., 172800., 172800., 172800.,
       172800., 259200., 259200., 259200., 259200., 259200., 345600.,
       345600., 345600., 345600., 345600.]

[erikvansebille]

Thanks for sending through, @Lachlan00 . I'm somewhat busy the coming days but maybe have time to have a look later this week. Can you in the meantime check what the contents of your out-*/* directory looks like?

[Lachlan00]

Thanks @erikvansebille,

Here's the out-*/* output and file tree:

out-KYYYLKCQ
.
└── 0
    ├── 0.npy
    └── pset_info.npy

out-KYYYLKCQ.zip

[erikvansebille]

OK thanks; that's just one time step. Could you, if you have time, also check what happens if you run the kernels independently? Is there a combination of kernels where the output file contains multiple timesteps again? Which combination? And is there one kernel that makes that only one timestep is outputted?

[Lachlan00]

HI @erikvansebille,

I ran the program with only the AdvectionRK4 kernel and the results were exactly the same. Curiously however if I run a longer simulation (the full run is 8309 days running with a dt at timedelta(days=1)) there are more numpy arrays produced but not the full time steps that would be expected. I keyboard interrupted the run at 10% (so ~ approximately 800 days - 67,478,400 of 717,811,200 advections). Here's the output and the file tree of that process.

out-SHFESALV
.
└── 0
    ├── 0.npy
    ├── 1.npy
    ├── 10.npy
    ├── 11.npy
    ├── 12.npy
    ├── 13.npy
    ├── 14.npy
    ├── 15.npy
    ├── 16.npy
    ├── 17.npy
    ├── 18.npy
    ├── 19.npy
    ├── 2.npy
    ├── 20.npy
    ├── 21.npy
    ├── 22.npy
    ├── 23.npy
    ├── 24.npy
    ├── 25.npy
    ├── 26.npy
    ├── 27.npy
    ├── 28.npy
    ├── 29.npy
    ├── 3.npy
    ├── 30.npy
    ├── 31.npy
    ├── 32.npy
    ├── 33.npy
    ├── 34.npy
    ├── 35.npy
    ├── 36.npy
    ├── 37.npy
    ├── 38.npy
    ├── 39.npy
    ├── 4.npy
    ├── 40.npy
    ├── 41.npy
    ├── 42.npy
    ├── 43.npy
    ├── 44.npy
    ├── 45.npy
    ├── 46.npy
    ├── 47.npy
    ├── 48.npy
    ├── 49.npy
    ├── 5.npy
    ├── 50.npy
    ├── 51.npy
    ├── 52.npy
    ├── 53.npy
    ├── 54.npy
    ├── 55.npy
    ├── 56.npy
    ├── 57.npy
    ├── 6.npy
    ├── 7.npy
    ├── 8.npy
    ├── 9.npy
    └── pset_info.npy

out-SHFESALV.zip

Let me know if there's any other outputs I can produce that might help the debugging process.

[erikvansebille]

OK, thanks for letting me know. I now heard that @daanreijnders also had a similar bug, so this seems like a serious bug. I will investigate asap

[erikvansebille]

Hi @Lachlan00 and @daanreijnders , I may have found the bug. Could you check whether the change in PR #672 fixes the writing issue for you?

[daanreijnders]

Yes! It seems to work for me. The .nc files contain the correct amount of particles and timesteps and they seem to behave correctly. Thanks!

[Lachlan00]

Hey @erikvansebille,

That seemed to do it! Working as expected on my end now. Thanks very much!

[erikvansebille]

I've updated the code to make it more robust to other use cases too. Could you check out the latest version of #672 and test again? If it's still ok, we will merge into master and probably release a new version very soon, as this is a serious bug that needs to be fixed in conda asap

[daanreijnders]

Still working for me (correct .nc dimensions and correct particle behavior when animated).

[Lachlan00]

@erikvansebille Your new patch appears to be all good on my end too. Particle trajectory file has the correct number of observations and trajectories.