[BUG/ENH] Making cross sections compatible with new scipy

Author: Miguel de la Varga

gempy/core/data/core_utils.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.interpolate import RegularGridInterpolator, RectBivariateSpline
 
 from gempy.optional_dependencies import require_scipy
 
@@ -20,39 +21,22 @@ def distance_2_points(p1, p2):
     return np.sqrt(np.diff((p1[0], p2[0])) ** 2 + np.diff((p1[1], p2[1])) ** 2)
 
 
-def interpolate_zvals_at_xy(xy, topography, method='interp2d'):
+def interpolate_zvals_at_xy(xy, topography, method='DEP'):
     """
-    Interpolates DEM values on a defined section
+    Interpolates DEM values on a defined section.
 
     Args:
-        xy: x (EW) and y (NS) coordinates of the profile
-        topography (:class:`gempy.core.grid_modules.topography.Topography`)
-        method: interpolation method, 'interp2d' for cubic scipy.interpolate.interp2d
-                                         'spline' for scipy.interpolate.RectBivariateSpline
+        xy (np.ndarray): Array of shape (n, 2) containing x (EW) and y (NS) coordinates of the profile.
+        topography (Topography): An instance of Topography containing the DEM data.
 
     Returns:
-        numpy.ndarray: z values, i.e. topography along the profile
-
+        np.ndarray: z values, i.e., topography along the profile.
     """
-
     xj = topography.values_2d[:, 0, 0]
     yj = topography.values_2d[0, :, 1]
     zj = topography.values_2d[:, :, 2]
-    scipy = require_scipy()
-    if method == 'interp2d':
-        f = scipy.interpolate.interp2d(xj, yj, zj.T, kind='cubic')
-        zi = f(xy[:, 0], xy[:, 1])
-        if xy[:, 0][0] <= xy[:, 0][-1] and xy[:, 1][0] <= xy[:, 1][-1]:
-            return np.diag(zi)
-        else:
-            return np.flipud(zi).diagonal()
-    else:
-        assert xy[:, 0][0] <= xy[:, 0][
-            -1], 'The xy values of the first point must be smaller than second.' \
-                 'Please use interp2d as method argument. Will be fixed.'
-        assert xy[:, 1][0] <= xy[:, 1][
-            -1], 'The xy values of the first point must be smaller than second.' \
-                 'Please use interp2d as method argument. Will be fixed.'
-        f = scipy.interpolate.RectBivariateSpline(xj, yj, zj)
-        zi = f(xy[:, 0], xy[:, 1])
-        return np.flipud(zi).diagonal()
+
+    spline = RectBivariateSpline(xj, yj, zj)
+    zi = spline.ev(xy[:, 0], xy[:, 1])
+    
+    return zi