ErrorOutOfBounds for tracer field sampling

[brogalla]

Hi Parcels Community,

I'm running into an ErrorOutOfBounds problem when sampling a tracer field based on particle tracks and was hoping that one of you might have an idea what the issue is. When I release particles at a particular location + time without field sampling, the particles run fine. However, running the same particles with my tracer sampling kernel, I get ErrorOutOfBounds for some, but not all particles. I've checked that the latitude and longitudes in the fieldsets for the tracer and velocity fields cover the sampled locations.

For my code, I tried to follow a method similar to this example notebook. You can see my notebook with the code that produces the problem here. I'm running with parcels version 2.4.0.

Thanks in advance!

[erikvansebille]

Thanks for reporting, @brogalla. I think the issue is with the depth of the particles. When you don't specify a depth in your ParticleSet, Parcels takes the shallowest depth of the U-Fields. But the Tracer field only starts deeper, as you see in your last cell:

print(np.amin(fieldset.U.depth), np.amax(fieldset.U.depth))
print(np.amin(fieldset.Pb.depth), np.amax(fieldset.Pb.depth))
0.0 5500.0015
0.49402538 5727.9165

So the solution could be to create your particles with e.g. depth=0.5?

[erikvansebille]

Hi @brogalla, I had a look at your Issue this afternoon, but I think this is a more fundamental problem for Parcels: it can't easily handle non-monotonous grids....

Your Grid is extremely curvilinear, which you can see if you issue

fieldset.Pb.show()

Now, try to fill your entire domain with Particles, by creating a (larger) regular grid of these and seeing which ones 'stick' when you delete all the ones out of bounds:

x = np.arange(-180, 180, 1)
y = np.arange(55, 89, 1)
plon, plat = np.meshgrid(x, y)
pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=plon, lat=plat)

pset.execute(AdvectionRK4_3D, runtime=0, dt=0, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

plt.plot(pset.lon, pset.lat, 'k.')

You see that all the ones west of approx 120W get deleted...

So my hunch is that this is not a problem with the Pb tracer field, but rather a problem that your grid is so extremely curved. The way that Parcels tries to find in which grid cell a particle is located (based on its lat/lon) just doesn't work for this Grid.

It will be a lot of under-the-hood development to fix this. So I'm very sorry, but I don't think that Parcels will work on your particular Grid...

[erikvansebille]

Actually, what I just wrote about difficult to fix this may not be entirely true. It occurred to me that you could help parcels by trying to find the indices of the grid cells yourself. Once these indices are correctly found, Parcels does a local search to update a grid cell when a particle moves out, but this is much simpler and more robust. So only the initial search is the problem

If you add the code below after you create the ParticleSet (but before the execute), things may work?

for i in range(len(pset)):
    d = np.hypot((pset.lon[i] - fieldset.U.grid.lon), (pset.lat[i] - fieldset.U.grid.lat))
    ind = np.unravel_index(np.argmin(d, axis=None), d.shape)
    pset.xi[i, :] = ind[1]
    pset.yi[i, :] = ind[0]

See also the image below, which now shows all particles in your curved domain

And it also works with sampling now:

import matplotlib.colors as cols
plt.scatter(pset.lon,pset.lat, c=pset.dissolpb, s=5, norm=cols.LogNorm())

[brogalla]

Thanks for looking into this more, Erik! Much appreciated. I did some more digging and it's odd... Although the curvature of the grid certainly could be part of the problem, I suspect it's not the main one. The subdomain and full domain have the exact same underlying grid (one is a subset of the other), so their curvature should be identical. However, the issue arises only with sampling the tracer field, so when I run particles with the velocities from the full domain and don't sample the tracer field, everything runs fine and OP doesn't have trouble finding the correct cell locations. You can see my debugging in this notebook.

I tried your nice suggestion of hinting at the coordinates for the start locations (I'll call it "nudging"). It definitely helps it start! When I run for a few timesteps, some of the parcels remain while the rest are removed because the correct cell location is not found within 100000 iterations:

When I run without indices specified in the Fieldset and without sampling the tracer, I don't need the nudging in order to run the particles in the entire region, while it was necessary when the indices were specified. For the next part, I've removed the indices in the Fieldset. If I run particles for multiple timesteps without sampling the tracer (so just with velocities from the full domain) and without nudging, it runs fine:

For the identical setup, if I add the tracer sampling kernel the subdomain is selected but excludes an area from the western end of the subdomain and it can't find a bunch of the particle locations (top figure), while if I nudge the start location the whole subdomain gets shifted which is odd (bottom figure):

So... I'm adding more confusion rather than reducing it 😅. I'll keep playing around with it; let me know if you have any further suggestions, but I also understand if this is just an odd use case.

[erikvansebille]

Hi @brogalla, I had another lok into your issue, and I think I may have found the problem. Your OutOfBounds kernel kills the particles also if they go through the ocean surface. This problem was also recently discussed in e.g. #1266, with solution here: #1266 (comment). Not sure why this is, but it may be related to the slight differences in the uppermost depth level of your U and Pb fields. Can you see if the solution above helps?

(by the way, I also found that the code for 'nudging' the starting locations was not perfect; as it only used the U and V grids; not the Pb grids. Especially when you don't use indexing, these can be very different. The proper code is below, where the correct index is found for each of the grids in the fieldset)

for i in range(len(pset)):
    for gi, g in enumerate(fieldset.gridset.grids):
        d = np.hypot((pset.lon[i] - g.lon), (pset.lat[i] - g.lat))
        ind = np.unravel_index(np.argmin(d, axis=None), d.shape)
        pset.xi[i, gi] = ind[1]
        pset.yi[i, gi] = ind[0]

[brogalla]

Woohoo! The solution from the through surface problem worked (and should be fine for my purposes) when combined with the corrected nudging! Thank you so much for the help, I really appreciate it and am happy to be able to use OP to look at tracer fields in another way.