Fixing OutOfBounds errors on field.eval (attempt 2) (#2681)

* Adding unit test for OutOfBounds field.eval

* Fixing OutOfBounds errors on field.eval

By masking out all the interpolated values to zero for particles out of bound

* Removing debug print statements

* Incorporating reviewer feedback

* Adding warning when field.eval is out of bounds

And filtering this warning in kernel execution

* Implementing reviewer feedback from @VeckoTheGecko

Author: erikvansebille

src/parcels/__init__.py

@@ -61,6 +61,7 @@
     "OutsideTimeInterval",
     "StatusCode",
     # Warnings
+    "FieldEvalWarning",
     "FieldSetWarning",
     "FileWarning",
     "KernelWarning",

src/parcels/_core/field.py

@@ -15,6 +15,7 @@
 )
 from parcels._core.utils.string import _assert_str_and_python_varname
 from parcels._core.uxgrid import UxGrid
+from parcels._core.warnings import FieldEvalWarning
 from parcels._core.xgrid import XGrid
 from parcels._typing import VectorType
 from parcels.interpolators._base import ScalarInterpolator, VectorInterpolator
@@ -175,6 +176,7 @@ def eval(self, time: datetime, z, y, x, particles=None):
         value = self.interp_method.interp(particle_positions, grid_positions, self)
 
         _update_particle_states_interp_value(particles, value)
+        _mask_outofbounds_values(grid_positions, value)
 
         return value
 
@@ -281,6 +283,7 @@ def eval(self, time: datetime, z, y, x, particles=None):
 
         for vel in (u, v, w):
             _update_particle_states_interp_value(particles, vel)
+            _mask_outofbounds_values(grid_positions, vel)
 
         if "3D" in self.vector_type:
             return (u, v, w)
@@ -349,6 +352,20 @@ def _update_particle_states_position(particles, grid_positions: dict):
             )
 
 
+def _mask_outofbounds_values(grid_positions: dict, value):
+    mask = np.zeros(value.shape, dtype=bool)
+    for dim in ["X", "Y", "Z", "FACE"]:
+        if dim in grid_positions:
+            mask[grid_positions[dim]["index"] < 0] = True
+    if np.any(mask):
+        warnings.warn(
+            "Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.",
+            FieldEvalWarning,
+            stacklevel=2,
+        )
+        value[mask] = 0.0
+
+
 def _update_particle_states_interp_value(particles, value):
     """Update the particle states based on the interpolated value, but only if state is not an Error already."""
     if particles:

src/parcels/_core/kernel.py

@@ -16,7 +16,7 @@
     _raise_grid_searching_error,
     _raise_outside_time_interval_error,
 )
-from parcels._core.warnings import KernelWarning
+from parcels._core.warnings import FieldEvalWarning, KernelWarning
 from parcels._python import assert_same_function_signature
 from parcels.kernels import (
     AdvectionAnalytical,
@@ -206,13 +206,16 @@ def execute(self, pset, endtime, dt):
 
             # run kernels for all particles that need to be evaluated
             for f in self._kernels:
-                f(pset[evaluate_particles], self._fieldset)
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore", FieldEvalWarning)
 
-                # check for particles that have to be repeated
-                repeat_particles = pset.state == StatusCode.Repeat
-                while np.any(repeat_particles):
-                    f(pset[repeat_particles], self._fieldset)
+                    f(pset[evaluate_particles], self._fieldset)
+
+                    # check for particles that have to be repeated
                     repeat_particles = pset.state == StatusCode.Repeat
+                    while np.any(repeat_particles):
+                        f(pset[repeat_particles], self._fieldset)
+                        repeat_particles = pset.state == StatusCode.Repeat
 
             # apply position/time update only to particles still in a normal state
             # (particles that signalled Stop*/Delete/errors should not have time/position advanced)

src/parcels/_core/warnings.py

@@ -28,6 +28,16 @@ class FileWarning(UserWarning):
     pass
 
 
+class FieldEvalWarning(UserWarning):
+    """Warning that is raised when there are issues during the evaluation of a Field.
+
+    These warnings can be related to out-of-bounds indices during interpolation,
+    or other issues that arise during the evaluation of a Field at particle positions.
+    """
+
+    pass
+
+
 class KernelWarning(RuntimeWarning):
     """Warning that is raised when there are issues with the Kernel.
 

src/parcels/interpolators/_uxinterpolators.py

@@ -177,5 +177,5 @@ def interp(
         if "3D" in vectorfield.vector_type:
             w = vectorfield.W.interp_method.interp(particle_positions, grid_positions, vectorfield.W)
         else:
-            w = 0.0
+            w = np.zeros_like(u)
         return u, v, w

src/parcels/interpolators/_xinterpolators.py

@@ -132,7 +132,7 @@ def interp(
         field: Field,
     ):
         """Returning the single value of a Constant Field (with a size=(1,1,1,1) array)"""
-        return field.data[0, 0, 0, 0].values
+        return field.data[0, 0, 0, 0].values * np.ones_like(particle_positions["lon"])
 
 
 class XLinear_Velocity(VectorInterpolator):  # noqa:  N801
@@ -155,7 +155,7 @@ def interp(
         if vectorfield.W:
             w = _xlinear.interp(particle_positions, grid_positions, vectorfield.W)
         else:
-            w = 0.0
+            w = np.zeros_like(u)
         return u, v, w
 
 
@@ -502,7 +502,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
 
         w *= f_w
     else:
-        w = None
+        w = np.zeros_like(u)
     return u, v, w
 
 

src/parcels/kernels/_sigmagrids.py

@@ -42,44 +42,44 @@ def AdvectionRK4_3D_CROCO(particles, fieldset):  # pragma: no cover
     It also uses linear interpolation of the W field, which gives much better results than the default C-grid interpolation.
     """
     dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt)
-    sigma = particles.z / fieldset.h[particles.time, 0, particles.lat, particles.lon]
+    sigma = particles.z / fieldset.h[particles.time, np.zeros_like(particles.z), particles.lat, particles.lon]
 
     sig = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles)
     (u1, v1) = fieldset.UV[particles.time, sig, particles.lat, particles.lon, particles]
     w1 = fieldset.W[particles.time, sig, particles.lat, particles.lon, particles]
-    w1 *= sigma / fieldset.h[particles.time, 0, particles.lat, particles.lon]
+    w1 *= sigma / fieldset.h[particles.time, np.zeros_like(particles.z), particles.lat, particles.lon]
     lon1 = particles.lon + u1 * 0.5 * dt
     lat1 = particles.lat + v1 * 0.5 * dt
     sig_dep1 = sigma + w1 * 0.5 * dt
-    dep1 = sig_dep1 * fieldset.h[particles.time, 0, lat1, lon1]
+    dep1 = sig_dep1 * fieldset.h[particles.time, np.zeros_like(particles.z), lat1, lon1]
 
     sig1 = convert_z_to_sigma_croco(fieldset, particles.time + 0.5 * dt, dep1, lat1, lon1, particles)
     (u2, v2) = fieldset.UV[particles.time + 0.5 * dt, sig1, lat1, lon1, particles]
     w2 = fieldset.W[particles.time + 0.5 * dt, sig1, lat1, lon1, particles]
-    w2 *= sig_dep1 / fieldset.h[particles.time, 0, lat1, lon1]
+    w2 *= sig_dep1 / fieldset.h[particles.time, np.zeros_like(particles.z), lat1, lon1]
     lon2 = particles.lon + u2 * 0.5 * dt
     lat2 = particles.lat + v2 * 0.5 * dt
     sig_dep2 = sigma + w2 * 0.5 * dt
-    dep2 = sig_dep2 * fieldset.h[particles.time, 0, lat2, lon2]
+    dep2 = sig_dep2 * fieldset.h[particles.time, np.zeros_like(particles.z), lat2, lon2]
 
     sig2 = convert_z_to_sigma_croco(fieldset, particles.time + 0.5 * dt, dep2, lat2, lon2, particles)
     (u3, v3) = fieldset.UV[particles.time + 0.5 * dt, sig2, lat2, lon2, particles]
     w3 = fieldset.W[particles.time + 0.5 * dt, sig2, lat2, lon2, particles]
-    w3 *= sig_dep2 / fieldset.h[particles.time, 0, lat2, lon2]
+    w3 *= sig_dep2 / fieldset.h[particles.time, np.zeros_like(particles.z), lat2, lon2]
     lon3 = particles.lon + u3 * dt
     lat3 = particles.lat + v3 * dt
     sig_dep3 = sigma + w3 * dt
-    dep3 = sig_dep3 * fieldset.h[particles.time, 0, lat3, lon3]
+    dep3 = sig_dep3 * fieldset.h[particles.time, np.zeros_like(particles.z), lat3, lon3]
 
     sig3 = convert_z_to_sigma_croco(fieldset, particles.time + dt, dep3, lat3, lon3, particles)
     (u4, v4) = fieldset.UV[particles.time + dt, sig3, lat3, lon3, particles]
     w4 = fieldset.W[particles.time + dt, sig3, lat3, lon3, particles]
-    w4 *= sig_dep3 / fieldset.h[particles.time, 0, lat3, lon3]
+    w4 *= sig_dep3 / fieldset.h[particles.time, np.zeros_like(particles.z), lat3, lon3]
     lon4 = particles.lon + u4 * dt
     lat4 = particles.lat + v4 * dt
     sig_dep4 = sigma + w4 * dt
 
-    dep4 = sig_dep4 * fieldset.h[particles.time, 0, lat4, lon4]
+    dep4 = sig_dep4 * fieldset.h[particles.time, np.zeros_like(particles.z), lat4, lon4]
     particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
     particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
     particles.dz += (

tests/test_advection.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pandas as pd
 import pytest
@@ -15,6 +17,7 @@
     convert,
 )
 from parcels._core.utils.time import timedelta_to_float
+from parcels._core.warnings import FieldEvalWarning
 from parcels._datasets.structured.generated import (
     decaying_moving_eddy_dataset,
     moving_eddy_dataset,
@@ -176,7 +179,9 @@ def SubmergeParticle(particles, fieldset):  # pragma: no cover
     kernels.append(DeleteParticle)
 
     pset = ParticleSet(fieldset=fieldset, lon=0.5, lat=0.5, z=0.9)
-    pset.execute(kernels, runtime=np.timedelta64(10, "s"), dt=np.timedelta64(1, "s"))
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", FieldEvalWarning)
+        pset.execute(kernels, runtime=np.timedelta64(10, "s"), dt=np.timedelta64(1, "s"))
 
     if direction == "up" and resubmerge_particle:
         np.testing.assert_allclose(pset.lon[0], 0.6, atol=1e-5)

tests/test_field.py

@@ -6,8 +6,11 @@
 from parcels import Field, VectorField
 from parcels._core.fieldset import FieldSet
 from parcels._core.model import StructuredModelData
+from parcels._core.warnings import FieldEvalWarning
+from parcels._datasets.structured.generated import simple_UV_dataset
 from parcels._datasets.structured.generic import T as T_structured
 from parcels._datasets.structured.generic import datasets as datasets_structured
+from parcels._datasets.unstructured.generic import _ux_constant_flow_face_centered_2D
 from parcels._datasets.unstructured.generic import datasets as datasets_unstructured
 from parcels.interpolators import (
     UxConstantFaceConstantZC,
@@ -180,6 +183,57 @@ def test_field_constant_in_time():
     assert np.isclose(P1, P2)
 
 
+@pytest.mark.parametrize(
+    "field_name, location, expected",
+    [
+        ("U", (0.0, 0.0, 0.0, 5e6), 0.0),
+        ("UV", (0.0, 0.0, 0.0, 5e6), [[0.0], [0.0]]),
+        ("U", (0.0, 0.0, 5e6, 0.0), 0.0),
+        ("UV", (0.0, 0.0, 5e6, 0.0), [[0.0], [0.0]]),
+        ("U", (0.0, 5e6, 0.0, 0.0), 0.0),
+        ("UV", (0.0, 5e6, 0.0, 0.0), [[0.0], [0.0]]),
+    ],
+)
+def test_field_eval_out_of_bounds_structured(field_name, location, expected):
+    """Test that Field.eval returns IndexError when queried outside the grid boundaries."""
+    # eval outside of bounds should return 0.0
+    ds = simple_UV_dataset(mesh="flat")
+    fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat")
+    fieldset.U.data[:] = 1.0
+    fieldset.V.data[:] = 2.0
+    field = getattr(fieldset, field_name)
+    with pytest.warns(
+        FieldEvalWarning,
+        match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.",
+    ):
+        np.testing.assert_allclose(field.eval(*location), expected)
+
+
+@pytest.mark.parametrize(
+    "field_name, location, expected",
+    [
+        ("U", (0.0, 0.0, 0.0, 5e6), 0.0),
+        ("UV", (0.0, 0.0, 0.0, 5e6), [[0.0], [0.0]]),
+        ("U", (0.0, 0.0, 5e6, 0.0), 0.0),
+        ("UV", (0.0, 0.0, 5e6, 0.0), [[0.0], [0.0]]),
+        ("U", (0.0, 5e6, 0.0, 0.0), 0.0),
+        ("UV", (0.0, 5e6, 0.0, 0.0), [[0.0], [0.0]]),
+    ],
+)
+def test_field_eval_out_of_bounds_unstructured(field_name, location, expected):
+    """Test that Field.eval returns IndexError when queried outside the grid boundaries."""
+    ds = _ux_constant_flow_face_centered_2D()
+    fieldset = FieldSet.from_ugrid_conventions(ds, mesh="flat")
+    fieldset.U.data[:] = 1.0
+    fieldset.V.data[:] = 2.0
+    field = getattr(fieldset, field_name)
+    with pytest.warns(
+        FieldEvalWarning,
+        match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.",
+    ):
+        np.testing.assert_allclose(field.eval(*location), expected)
+
+
 def test_field_unstructured_grid_creation(): ...