feat: opt-in local-radius RBF for InterpolatedStructure dip interpolation

map2loop/interpolators.py

@@ -1,9 +1,9 @@
 from abc import ABC, abstractmethod
-from typing import Any, Union
+from typing import Any, Optional, Union
 import beartype
 import numpy
 from numpy import ndarray
-from scipy.interpolate import Rbf, LinearNDInterpolator
+from scipy.interpolate import Rbf, LinearNDInterpolator, RBFInterpolator
 from sklearn.cluster import DBSCAN
 import pandas
 
@@ -303,9 +303,17 @@ class DipDipDirectionInterpolator(Interpolator):
         Interpolator(ABC): Derived from Abstract Base Class
     """
 
-    def __init__(self, data_type=None):
+    def __init__(self, data_type=None, neighbors: Optional[int] = None):
         """
         Initialiser of for IDWInterpolator
+
+        Args:
+            data_type: which of "dip"/"dipdir" to interpolate, defaults to both
+            neighbors: if set, and RBFInterpolator is requested via interpolate(),
+                restrict each grid point's fit to its `neighbors` nearest data
+                points instead of blending every measurement across the whole
+                map. This keeps a fold hinge or fault-bounded change in dip from
+                being smoothed into neighbouring, structurally unrelated areas.
         """
         if data_type is None:
             self.data_type = ["dip", "dipdir"]
@@ -318,6 +326,7 @@ def __init__(self, data_type=None):
         self.dip = None
         self.dipdir = None
         self.cell_size = None
+        self.neighbors = neighbors
         self.interpolator_label = "DipDipDirectionInterpolator"
 
     def type(self):
@@ -392,7 +401,10 @@ def interpolate(self, ni: Union[ndarray, list], interpolator: Any = Rbf):
 
         Args:
             ni (int): value to interpolate
-            interpolator: type of interpolator to use by default SciPy Rbf interpolator
+            interpolator: type of interpolator to use by default SciPy Rbf interpolator.
+                Pass scipy.interpolate.RBFInterpolator to fit each grid point using only
+                its `self.neighbors` nearest data points (a local radial basis function)
+                instead of a single surface fit to every point across the whole map.
 
         Returns:
             Rbf: radial basis function object
@@ -401,6 +413,16 @@ def interpolate(self, ni: Union[ndarray, list], interpolator: Any = Rbf):
             rbf = Rbf(self.x, self.y, ni, function="linear")
             return rbf(self.xi, self.yi)
 
+        if interpolator is RBFInterpolator:
+            points = numpy.column_stack([self.x, self.y])
+            query_points = numpy.column_stack([self.xi, self.yi])
+            # neighbors=None (the RBFInterpolator default) fits one global surface,
+            # equivalent in spirit to Rbf; a finite value restricts each query point's
+            # fit to its nearest neighbours, keeping the interpolation local.
+            neighbors = self.neighbors if self.neighbors is None else min(self.neighbors, len(self.x))
+            rbf = RBFInterpolator(points, numpy.asarray(ni), neighbors=neighbors, kernel="linear")
+            return rbf(query_points)
+
         if interpolator is LinearNDInterpolator:
             lnd_interpolator = LinearNDInterpolator(list(zip(self.x, self.y)), ni)
             return lnd_interpolator(self.xi, self.yi)
@@ -415,7 +437,8 @@ def __call__(
         Args:
             bounding_box (dict): a dictionary containing the bounding box of the map data
             structure_data (pandas.DataFrame): sampled structural data
-            interpolator (Union[Rbf, LinearNDInterpolator]): type of interpolator to use by default SciPy Rbf interpolator
+            interpolator (Union[Rbf, RBFInterpolator, LinearNDInterpolator]): type of interpolator to
+                use, by default SciPy Rbf interpolator
 
         Returns:
             numpy.ndarray: interpolated dip and dip direction values

map2loop/thickness_calculator.py

@@ -229,18 +229,29 @@ class InterpolatedStructure(ThicknessCalculator):
     """
 
     def __init__(
-        self, 
-        dtm_data: Optional[gdal.Dataset] = None, 
-        bounding_box: Optional[dict] = None, 
+        self,
+        dtm_data: Optional[gdal.Dataset] = None,
+        bounding_box: Optional[dict] = None,
         max_line_length: Optional[float] = None,
-        is_strike: Optional[bool] = False
+        is_strike: Optional[bool] = False,
+        local_interpolation_neighbors: Optional[int] = None,
         ):
         """
         Initialiser for interpolated structure version of the thickness calculator
+
+        Args:
+            local_interpolation_neighbors: if set, interpolate dip using a local
+                radial basis function fit from only this many nearest structure
+                points per grid location, instead of one surface fit to every
+                structure point across the whole map. This avoids smoothing dip
+                across fold hinges, faults or unrelated structural domains, at
+                the cost of a less smooth interpolated surface. Defaults to None
+                (whole-map interpolation, matching previous behaviour).
         """
         super().__init__(dtm_data, bounding_box, max_line_length, is_strike)
         self.thickness_calculator_label = "InterpolatedStructure"
         self.lines = None
+        self.local_interpolation_neighbors = local_interpolation_neighbors
 
     @beartype.beartype
     def compute(
@@ -317,9 +328,15 @@ def compute(
         if 'Z' in contacts.columns:
             contacts = contacts[["X", "Y", "Z", "geometry", "basal_unit"]].copy()
         # Interpolate the dip of the contacts
-        interpolator = DipDipDirectionInterpolator(data_type="dip")
-        # Interpolate the dip of the contacts
-        dip = interpolator(self.bounding_box, structure_data, interpolator=scipy.interpolate.Rbf)
+        interpolator = DipDipDirectionInterpolator(
+            data_type="dip", neighbors=self.local_interpolation_neighbors
+        )
+        if self.local_interpolation_neighbors is not None:
+            dip = interpolator(
+                self.bounding_box, structure_data, interpolator=scipy.interpolate.RBFInterpolator
+            )
+        else:
+            dip = interpolator(self.bounding_box, structure_data, interpolator=scipy.interpolate.Rbf)
         # create a GeoDataFrame of the interpolated orientations
         interpolated_orientations = geopandas.GeoDataFrame()
         # add the dip and dip direction to the GeoDataFrame