Interpolation at coastlines

[Abobie]

I just started looking into the problem with boundary values and I don't really remember/understand what the problem is... I tried setting the top right quarter of the peninsula grid to NaN. I used a grid size of (20,20) and let it run for a bit longer (36 hours) and got the results displayed in the picture below (the background is the U field). Now, could you please explain to me again what the problem is here? I was under the impression that if some values in the grid were NaN it would mess up the entire grid, but I think the results look pretty good. It even manages to do some nice interpolation when the particles are in a grid cell between NaN and non-NaN, like the cyan particle.

[mlange05]

Yes, we currently remove the NaN values and replace them with 0s (see parcels/field.py:49). This, however, is just a temporary hack that avoids all kinds of problems but is not sustainable long term, since we do need detect boundaries properly. This should be done by either setting them to NaN or using masks to skip the unwanted values.

[erikvansebille]

Michael is right, the replacing of NaNs with zeros clearly is a hack. The issue of boundary conditions is closely related to the issue of spatial interpolation. Perhaps we'll have to think about more advanced interpolation schemes. I've heard colleagues talk about how Wachpress interpolation can assure that particles never hit coastlines, but must say I'm not completely sure what Wachpress interpolation is and whether it is what we are looking for. Joel, could you investigate?

[mlange05]

For completeness, the reason we do this right now is SciPy issue #2360; which comes down to interpolator splines in FITPACK being global and thus choking as soon as there is a single NaN in the data.

[Abobie]

Ok, so are you saying we might have to write our own interpolation script? I was hoping I could just post a question to stackoverflow for a solution, maybe just to use another method from scipy.interpolate, but it sounds like it's not quite that simple. Anyway, I'll see if I can find something useful on Wachspress interpolation.

[erikvansebille]

Searching around is always a good idea. And it doesn't hurt to ask a question to stackoverflow, of course. Maybe someone has found a solution already

[mlange05]

Ok, after some further digging I think we should attempt to create a workaround for this, which I believe requires us to to two things:

• Python-mode: Make sure the auto-inserted values do not affect the accuracy of the generated SciPy interpolator object. This can either be done by performing interpolation for the missing value, similar to the way Pandas gets around this, or adding 0-weights to the spline creation to mask the missing values. The latter might require us to change the type of interpolator we are using (scipy.interpolation.LSQBivariateSpline allows you to define weights).
• JIT and Python: Keep a mask/bitmask to signal NaN/missing values and test this before evaluating fields. This will then allow us to skip interpolation and do whatever we deem appropriate in this case.

[erikvansebille]

Closing this old Issue, as most of it is redundant now.

With the new interpolation schemes in release v2.0.0beta, particles should not get stuck anymore if the flow fields is non-divergent near a coastline. note that this is very often not true, especially when model outputs have been reinterpolated before being shared. But in these cases, there is nothing much that Parcels can do about stuck particles. One option is to then gently push particles back to the ocean is they're on land