DEP: replace deprecated scipy interp1d for cube attr

src/xtgeo/surface/_regsurf_cube_window_v3.py

@@ -6,7 +6,6 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from scipy.interpolate import interp1d
 
 from xtgeo.common.log import null_logger
 
@@ -141,8 +140,8 @@ def _create_depth_cube(cube: Cube) -> np.ndarray:
     return dcube
 
 
-def _refine_cubes_vertically(cvalues, dvalues, ndiv):
-    """Refine the cubes vertically for better resolution"""
+def _refine_cubes_vertically(cvalues, dvalues, ndiv=2):
+    """Refine the cubes vertically for better resolution."""
     if not ndiv:
         ndiv = 2  # default
     logger.info("Resampling vertically, according to ndiv = %s", ndiv)
@@ -151,28 +150,24 @@ def _refine_cubes_vertically(cvalues, dvalues, ndiv):
         return cvalues, dvalues
 
     logger.info("Original shape is %s", cvalues.shape)
-    cref = np.random.rand(cvalues.shape[0], cvalues.shape[1], ndiv * cvalues.shape[2])
-    dref = cref.copy()
-
-    num_points = cref.shape[-1]
-    # Create interpolation function for cube values
-    interp_func1 = interp1d(
-        np.arange(cvalues.shape[-1]),
-        cvalues,
-        axis=-1,
-        kind="linear",
-        fill_value="extrapolate",
-    )
-    interp_func2 = interp1d(
-        np.arange(dvalues.shape[-1]),
-        dvalues,
-        axis=-1,
-        kind="linear",
-        fill_value="extrapolate",
-    )
-    # Resample array2 to match the number of points in array1
-    cref = interp_func1(np.linspace(0, cvalues.shape[-1] - 1, num_points))
-    dref = interp_func2(np.linspace(0, dvalues.shape[-1] - 1, num_points))
+
+    # Determine the new shape
+    new_shape = (cvalues.shape[0], cvalues.shape[1], ndiv * cvalues.shape[2])
+    cref = np.zeros(new_shape)
+    dref = np.zeros(new_shape)
+
+    # Compute new indices for interpolation
+    new_x = np.linspace(0, cvalues.shape[2] - 1, new_shape[2])
+
+    # Interpolate for each layer and column
+    for i in range(cvalues.shape[0]):
+        for j in range(cvalues.shape[1]):
+            cref[i, j, :] = np.interp(
+                new_x, np.arange(cvalues.shape[2]), cvalues[i, j, :]
+            )
+            dref[i, j, :] = np.interp(
+                new_x, np.arange(dvalues.shape[2]), dvalues[i, j, :]
+            )
 
     logger.info("Resampling done, new shape is %s", cref.shape)
     return cref, dref

tests/test_surface/test_surface_cube_attrs.py

@@ -197,6 +197,64 @@ def test_various_attrs_algorithm3(loadsfile1):
     assert surfx.values.mean() == pytest.approx(529.328, abs=0.01)
 
 
+@pytest.mark.parametrize(
+    "ndiv, expected_mean",
+    (
+        [2, 177.35],
+        [4, 176.12],
+        [12, 176.43],
+    ),
+)
+def test_various_attrs_algorithm3_ndiv(loadsfile1, ndiv, expected_mean):
+    cube1 = loadsfile1
+    surf1 = xtgeo.surface_from_cube(cube1, 2540)
+    surf2 = xtgeo.surface_from_cube(cube1, 2548)
+
+    surfx = surf1.copy()
+    surfx.slice_cube_window(
+        cube1,
+        other=surf2,
+        other_position="below",
+        attribute="mean",
+        sampling="trilinear",
+        snapxy=True,
+        ndiv=ndiv,
+        algorithm=3,
+    )
+
+    assert surfx.values.mean() == pytest.approx(expected_mean, abs=0.1)
+
+
+@pytest.mark.bigtest
+@pytest.mark.parametrize(
+    "ndiv, expected_mean",
+    (
+        [2, 0.0013886],
+        [4, 0.0013843],
+        [12, 0.0013829],
+    ),
+)
+def test_various_attrs_algorithm3_ndiv_large(loadsfile2, ndiv, expected_mean):
+    cube1 = loadsfile2
+    surf1 = xtgeo.surface_from_cube(cube1, 1560)
+    surf2 = xtgeo.surface_from_cube(cube1, 1760)
+
+    surfx = surf1.copy()
+    surfx.slice_cube_window(
+        cube1,
+        other=surf2,
+        other_position="below",
+        attribute="mean",
+        sampling="trilinear",
+        snapxy=True,
+        ndiv=ndiv,
+        algorithm=3,
+    )
+
+    print(surfx.values.mean())
+    assert surfx.values.mean() == pytest.approx(expected_mean, abs=0.00001)
+
+
 def test_avg_surface(loadsfile1):
     cube1 = loadsfile1
     surf1 = xtgeo.surface_from_cube(cube1, 1100.0)