geoutils.operators.Interpolator

class geoutils.operators.Interpolator(neighborhood=None)

Interpolator to continuously predict values at new coordinates from nearby observations.

An Interpolator assumes a spatial field that can be predicted continuously at a point. It is tied to a neighborhood, a nodata propagation scheme, and has optional support for sample weights, when the method accepts them. Use a Reducer when the prediction should reduce observations across a neighborhood instead (e.g., average over an area). Fractional area support exists only for Reducers, not Interpolators.

To define your own interpolator, implement predict(). It will receive a LocalData object with nearby values, their X/Y coordinates and the target X/Y coordinate. Implement coefficients() instead when the prediction is a weighted sum that can also be used for uncertainty propagation.

For example, the custom interpolator below fits a plane to nearby raster cell values:

import numpy as np
import geoutils as gu

class LocalPlane(gu.operators.Interpolator):
    minimum_inputs = 3

    def predict(self, data):
        # Calculate the X and Y differences between each source cell and the requested point
        x_from_target = data.coordinates[:, 0] - data.target[0]
        y_from_target = data.coordinates[:, 1] - data.target[1]

        # Fit value = value_at_target + x_slope * x_from_target + y_slope * y_from_target
        plane_terms = np.column_stack((np.ones(len(data.values)), x_from_target, y_from_target))
        value_at_target, x_slope, y_slope = np.linalg.lstsq(plane_terms, data.values, rcond=None)[0]
        return float(value_at_target)

predicted = raster.interp_at_points((x, y), method=LocalPlane(), as_array=True)

A linear interpolator should implement coefficients() to support weights and uncertainty propagation. As an example below, a simple inverse-distance interpolation can be naturally written as a weighted sum:

class SimpleInverseDistance(gu.operators.Interpolator):
    '''Predict from nearby cells, giving closer cells more weight.'''

    def coefficients(self, data):
        '''Return normalized inverse-distance weights for each cell.'''

        # A cell at the target takes precedence and avoids division by zero
        exact = data.distances == 0
        if np.any(exact):
            weights = exact.astype(float)
        else:
            weights = 1 / data.distances
        return gu.operators.LinearCoefficients(weights / weights.sum())

errors = gu.ErrorStructure([gu.ErrorComponent("measurement", 2.0)])
summary = gu.uncertainty.propagate(
    raster.interp_at_points, error_structure=errors,
    operation_kwargs=dict(points=(x, y), method=SimpleInverseDistance(), as_array=True),
)
predicted = summary.estimate

An Interpolator can also limit which nearby observations it uses. The following example uses at most eight points within 50 coordinate units of each output cell:

nearby = gu.operators.PointNeighbours(max_points=8, radius=50)
gridded = point_cloud.grid(
    ref=reference_raster, resampling=SimpleInverseDistance(neighborhood=nearby)
)

__init__(neighborhood=None)

Choose how nearby points or raster cells are selected.

Parameters:

neighborhood (PointNeighbours | GridNeighbours | None) – PointNeighbours for a point cloud or GridNeighbours for a raster. By default, the input type determines which to use (eight nearest points or a 3 x 3 window for a custom interpolator).

Methods

__init__([neighborhood])	Choose how nearby points or raster cells are selected.
coefficients(data)	Return the weights and optional constant offset for a prediction.
evaluate(data, *[, nodata_propagation])	Check nodata and weights before predicting a value.
predict(data)	Predict a value from the nearby source observations.
predict_batch(data, *[, nodata_propagation])	Predict several values by calling evaluate() for each group of source observations.

Attributes

accepts_sample_weights
default_neighborhood
default_nodata_propagation
error_structure
interpolation_order
minimum_inputs