```
# 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 interpolation2d.py 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
```