Source code for safe.gis.interpolation1d

# coding=utf-8
"""**Module for 1D interpolation**

This module:

* provides piecewise constant (nearest neighbour) and bilinear interpolation
* is fast (based on numpy vector operations)
* depends only on numpy
* guarantees that interpolated values never exceed the two nearest neighbours
* handles missing values in domain sensibly using NaN
* is unit tested with a range of common and corner cases

See for documentation of the mathematical derivation used.

import numpy

from safe.gis.interpolation import validate_inputs, validate_mode
# pylint: disable=W0105

# noinspection PyArgumentEqualDefault,PyTypeChecker
[docs]def interpolate1d(x, z, points, mode='linear', bounds_error=False): """Fundamental 1D interpolation routine. :param x: 1D array of x-coordinates on which to interpolate :type x: numpy.ndarray :param z: 1D array of values for each x :type z: numpy.ndarray :param points: 1D array of coordinates where interpolated values are sought :type points: numpy.ndarray :param mode: Determines the interpolation order. Options are: * 'constant' - piecewise constant nearest neighbour interpolation * 'linear' - bilinear interpolation using the two nearest \ neighbours (default) :type mode: str :param bounds_error: Flag to indicate whether an exception will be raised when interpolated values are requested outside the domain of the input data. If False, nan is returned for those values. :type bounds_error: bool :returns: 1D array with same length as points with interpolated values :rtype: numpy.ndarry :raises: RuntimeError ..note:: Input coordinates x are assumed to be monotonically increasing, but need not be equidistantly spaced. z is assumed to have dimension M where M = len(x). """ # Input checks validate_mode(mode) # pylint: disable=unbalanced-tuple-unpacking x, z, xi = validate_inputs( x=x, z=z, points=points, bounds_error=bounds_error) # Identify elements that are outside interpolation domain or NaN outside = (xi < x[0]) + (xi > x[-1]) outside += numpy.isnan(xi) inside = -outside xi = xi[inside] # Find upper neighbours for each interpolation point idx = numpy.searchsorted(x, xi, side='left') # Internal check (index == 0 is OK) msg = ('Interpolation point outside domain. This should never happen. ' 'Please email [email protected]') if len(idx) > 0: if not max(idx) < len(x): raise RuntimeError(msg) # Get the two neighbours for each interpolation point x0 = x[idx - 1] x1 = x[idx] z0 = z[idx - 1] z1 = z[idx] # Coefficient for weighting between lower and upper bounds alpha = (xi - x0) / (x1 - x0) if mode == 'linear': # Bilinear interpolation formula dx = z1 - z0 zeta = z0 + alpha * dx else: # Piecewise constant (as verified in input_check) # Set up masks for the quadrants left = alpha < 0.5 # Initialise result array with all elements set to right neighbour zeta = z1 # Then set the left neigbours zeta[left] = z0[left] # Self test if len(zeta) > 0: mzeta = numpy.nanmax(zeta) mz = numpy.nanmax(z) # noinspection PyStringFormat msg = ('Internal check failed. Max interpolated value %.15f ' 'exceeds max grid value %.15f ' % (mzeta, mz)) if not(numpy.isnan(mzeta) or numpy.isnan(mz)): if not mzeta <= mz: raise RuntimeError(msg) # Populate result with interpolated values for points inside domain # and NaN for values outside r = numpy.zeros(len(points)) r[inside] = zeta r[outside] = numpy.nan return r # Mathematical derivation of the interpolation formula used # noinspection PyStatementEffect
""" Bilinear interpolation is based on the standard 1D linear interpolation formula: Given points (x0, y0) and (x1, x0) and a value of x where x0 <= x <= x1, the linearly interpolated value y at x is given as alpha*(y1-y0) + y0 or alpha*y1 + (1-alpha)*y0 (1) where alpha = (x-x0)/(x1-x0) (1a) Piecewise constant interpolation can be implemented using the same coefficients (1a) that are used for bilinear interpolation as they are a measure of the relative distance to the left and lower neigbours. A value of 0 will pick the lower bound whereas a value of 1 will pick the higher bound. Hence z can be assigned to its nearest neigbour as follows | z0 alpha < 0.5 # lower value z = | | z1 alpha >= 0.5 # higher value """ # pylint: enable=W0105