Raise a Warning or Error when lon/lat in grids are not monotonically increasing

[erikvansebille]

As shown in https://github.com/OceanParcels/parcels/issues/1360#issuecomment-1738535942, Parcels can throw a confusing OutOfBoundsError when there are negative gradients in the longitude or latitude grids.

It would be good to capture these on FieldSet creation, so that users know that they can expect erroneous behaviour.

[Anna7651]

Hi, I would like to contribute to this issue if it’s still open. I have tried to add such a warning message in my fork in the FieldSet class, but I’m not sure if I understood the task correctly. So I’d be happy if you let me know how you want it to be implemented, and then I can try to work on it. (https://github.com/Anna7651/ParcelsF/tree/lonLatMonotone)

[VeckoTheGecko]

Hi @Anna7651 , thanks for your interest in this issue. I've taken a look at your fork, and I think we can make the diff smaller. I think something along the lines of the below would greatly simplify things. Thoughts?

import numpy as np

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)

X, Y = np.meshgrid(x, y)

def is_monotonic(array: np.ndarray, increasing: bool = True):
    """Checks that a provided array is all monotonically increasing/decreasing along all axes."""
    if increasing:
        def f(arr): return arr >= 0 
    else:
        def f(arr): return arr <= 0
    
    for axis in range(array.ndim):
        monotonic_along_axis = f(np.diff(array, axis=axis)).all()
        if not monotonic_along_axis:
            return False
    return True


assert is_monotonic(x)
assert is_monotonic(X)
assert not is_monotonic(-x)
assert not is_monotonic(-y)
assert not is_monotonic(-X)
assert not is_monotonic(-Y)

Z = X ** 2 + Y ** 2
assert not is_monotonic(Z)

Regarding the "monotonically increasing" argument (cc @erikvansebille ), I'm not 100% sure if this issue is correctly scoped. If we have an arbitrary curvilinear grid, do the values necessarily need to be monotonically increasing? See the image below from the original discussion, where the latitude isn't monotonically increasing

@Anna7651 I think a minimal example here where the issue is highlighted would help with the scoping of this issue, if that's something that you're interested in documenting :)

cc @michaeldenes since you were also running into this in your work (ie. "index couldn't be found after 10000 iterations")

[JamiePringle]

I do not know what Parcels requires. But in general, there is no requirement that a ROMS grid have monotonic increasing coordinates on a grid. Long ago I configured a grid as torrid... Fun fact -- all that ROMS needs to run is the grid metrics, (dx, dy, and the curvature matrics dx/d_eta and dy/d_xi)

@VeckoTheGecko Having said that, the grid above is monotonic in the sense (and there are several senses you can define monotonicity) that for the model coordinates (xi, eta), on a constant xi the coordinate longitude always increases, and on a constant eta the latitude always increases. (I may have xi and eta flipped).

Most grids are relatively straightforward... but check out this one from a very good discussion of ROMS grids at https://www.myroms.org/forum/viewtopic.php?t=295 . The east/west coordinate on the right-hand side of the grid is clearly non-monotonic, and this is allowed.

It would be helpful to precisely define what kind of monotonicity Parcels needs.

[erikvansebille]

Thanks for this input, @JamiePringle, @VeckoTheGecko and @Anna7651! I agree that perhaps the issue is not that 'non-monotonic' grids are not parsed well, but rather that parcels sometimes has difficulties with 'complicated' grids (however defined). I remember trying to create a simple minimal example that breaks; and couldn't create a grid that would cause the errors in https://github.com/OceanParcels/Parcels/issues/1360#issuecomment-1738535942.

I propose to pause this development fore now; and wait until we either encounter a simple minimal example where this breaks; or get support for unstructured grids into (hopefully!) Parcels v4.