fix for issue #155

docs/api/changelog.rst

@@ -15,6 +15,8 @@ Fixed
   results in invalid MODFLOW6 input for 2D grid data, such as DIS top.
 - :meth:`imod.flow.ImodflowModel.write` no longer writes incorrect projectfiles
   for non-grid values with a time and layer dimension.
+- :func:`imod.evaluate.interpolate_value_boundaries`: Fix edge case when
+  successive values in z direction are exactly equal to the boundary value.
 
 Changed
 ~~~~~~~

imod/evaluate/boundaries.py

@@ -17,6 +17,9 @@ def _interpolate_value_boundaries_z1d(values, z, dz, threshold, out):
                 if (n % 2 == 0 and values[k, i, j] >= threshold) or (
                     n % 2 == 1 and values[k, i, j] <= threshold
                 ):  # exceeding (n even) or falling below threshold (n odd)
+                    if values[k, i, j] == values[k - 1, i, j]:
+                        # edge case: subsequent values equal threshold result in ZeroDivisionError
+                        continue
                     if not np.isnan(values[k - 1, i, j]):
                         # interpolate z coord of threshold
                         out[n, i, j] = z[k - 1] + (threshold - values[k - 1, i, j]) * (
@@ -42,6 +45,9 @@ def _interpolate_value_boundaries_z3d(values, z, dz, threshold, out):
                 if (n % 2 == 0 and values[k, i, j] >= threshold) or (
                     n % 2 == 1 and values[k, i, j] <= threshold
                 ):  # exceeding (n even) or falling below threshold (n odd)
+                    if values[k, i, j] == values[k - 1, i, j]:
+                        # edge case: subsequent values equal threshold result in ZeroDivisionError
+                        continue
                     if not np.isnan(values[k - 1, i, j]):
                         # interpolate z coord of threshold
                         out[n, i, j] = z[k - 1, i, j] + (
@@ -69,6 +75,10 @@ def _interpolate_value_boundaries_z3ddz1d(values, z, dz, threshold, out):
                 if (n % 2 == 0 and values[k, i, j] >= threshold) or (
                     n % 2 == 1 and values[k, i, j] <= threshold
                 ):  # exceeding (n even) or falling below threshold (n odd)
+                    # if (values[k, i, j] == threshold) and (values[k - 1, i, j] == threshold):
+                    if values[k, i, j] == values[k - 1, i, j]:
+                        # edge case: subsequent values equal threshold result in ZeroDivisionError
+                        continue
                     if not np.isnan(values[k - 1, i, j]):
                         # interpolate z coord of threshold
                         out[n, i, j] = z[k - 1, i, j] + (
@@ -111,18 +121,25 @@ def interpolate_value_boundaries(values, z, threshold):
 
     values = values.load()
     out = xr.full_like(values, np.nan)
+    # if z.dz is a scalar, conform to z
+    # TODO: broadcast dz to z, so _interpolate_value_boundaries_z3ddz1d is superfluous?
+    if len(z.dz.dims) == 0:
+        dz = np.ones_like(z) * z.dz.values
+    else:
+        dz = z.dz.values
+
     if len(z.dims) == 1:
         _interpolate_value_boundaries_z1d(
-            values.values, z.values, z.dz.values, threshold, out.values
+            values.values, z.values, dz, threshold, out.values
         )
     else:
-        if len(z.dz.dims) == 1:
+        if len(z.dz.dims) <= 1:
             _interpolate_value_boundaries_z3ddz1d(
-                values.values, z.values, z.dz.values, threshold, out.values
+                values.values, z.values, dz, threshold, out.values
             )
         else:
             _interpolate_value_boundaries_z3d(
-                values.values, z.values, z.dz.values, threshold, out.values
+                values.values, z.values, dz, threshold, out.values
             )
     out = out.rename({"layer": "boundary"})
     out = out.dropna(dim="boundary", how="all")

imod/tests/test_evaluate/test_boundaries.py

@@ -126,3 +126,61 @@ def test_interpolate_value_boundaries_voxels():
 
     assert exc_ref.round(5).equals(exc.round(5))
     assert fallb_ref.round(5).equals(fallb.round(5))
+
+
+def test_interpolate_value_boundaries_scalardz():
+    data = np.array([[[0.5], [1.2]], [[1.0], [0.5]], [[1.2], [3.0]]])
+    z = np.array([-1.0, -3.0, -5.0])
+    dz = np.array(2.0)
+    coords = {"layer": [1, 2, 3], "y": [0.5, 1.5], "x": [0.5]}
+    zcoords = {"layer": [1, 2, 3]}
+    dims = ("layer", "y", "x")
+    zdims = "layer"
+    da = xr.DataArray(data, coords, dims)
+    z = xr.DataArray(z, zcoords, zdims)
+    z = z.assign_coords({"dz": dz})
+
+    exc_ref = np.array([[[-3.0], [0.0]], [[np.nan], [-3.4]]])
+    fallb_ref = np.array([[[np.nan], [-1.57142857]]])
+    coords2 = {"boundary": [0, 1], "y": [0.5, 1.5], "x": [0.5]}
+    dims2 = ("boundary", "y", "x")
+    exc_ref = xr.DataArray(exc_ref, coords2, dims2)
+    coords2["boundary"] = [0]
+    fallb_ref = xr.DataArray(fallb_ref, coords2, dims2)
+
+    exc, fallb = imod.evaluate.interpolate_value_boundaries(da, z, 1.0)
+
+    assert exc_ref.round(5).equals(exc.round(5))
+    assert fallb_ref.round(5).equals(fallb.round(5))
+
+
+def test_interpolate_value_boundaries_samevalues():
+    # sequence of same data - not equal and equal to boundary
+    data = np.array([[[1.0], [1.0]], [[1.0], [1.0]], [[1.0], [3.0]]])
+    z = np.array([[[0.5], [0.5]], [[-1.0], [-1.0]], [[-5.0], [-5.0]]])
+    dz = np.array([[[1.0], [1.0]], [[2.0], [2.0]], [[6.0], [6.0]]])
+    coords = {"layer": [1, 2, 3], "y": [0.5, 1.5], "x": [0.5]}
+    dims = ("layer", "y", "x")
+    da = xr.DataArray(data, coords, dims)
+    z = xr.DataArray(z, coords, dims)
+    z = z.assign_coords({"dz": (("layer", "y", "x"), dz)})
+
+    exc_ref = np.array([[[1.0], [1.0]]])
+    fallb_ref = np.empty((0, 2, 1))
+    coords2 = {"boundary": [0], "y": [0.5, 1.5], "x": [0.5]}
+    dims2 = ("boundary", "y", "x")
+    exc_ref = xr.DataArray(exc_ref, coords2, dims2)
+    coords2["boundary"] = []
+    fallb_ref = xr.DataArray(fallb_ref, coords2, dims2)
+
+    # values not equal to boundary
+    exc, fallb = imod.evaluate.interpolate_value_boundaries(da, z, 0.9)
+
+    assert exc_ref.round(5).equals(exc.round(5))
+    assert fallb_ref.round(5).equals(fallb.round(5))
+
+    # values equal to boundary
+    exc, fallb = imod.evaluate.interpolate_value_boundaries(da, z, 1.0)
+
+    assert exc_ref.round(5).equals(exc.round(5))
+    assert fallb_ref.round(5).equals(fallb.round(5))

setup.cfg

@@ -81,10 +81,14 @@ pytest11 =
 
 [flake8]
 ignore =
-    E203 # whitespace before ':' - doesn't work well with black
-    E402 # module level import not at top of file
-    E501 # line too long - let black worry about that
-    W503 # line break before binary operator
+    # whitespace before ':' - doesn't work well with black
+    E203,
+    # module level import not at top of file
+    E402,
+    # line too long - let black worry about that
+    E501,
+    # line break before binary operator
+    W503, 
 per-file-ignores = 
     ./docs/conf.py:F401
     __init__.py:F401