diff --git a/map2loop/interpolators.py b/map2loop/interpolators.py index 9648348e..6db5548f 100644 --- a/map2loop/interpolators.py +++ b/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 diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py index 483975d0..27549ace 100644 --- a/map2loop/thickness_calculator.py +++ b/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