python中griddata的外插值_griddata二维插值

"""Simple N-D interpolation

.. versionadded:: 0.9"""

#

#Copyright (C) Pauli Virtanen, 2010.#

#Distributed under the same BSD license as Scipy.#

#

#Note: this file should be run through the Mako template engine before#feeding it to Cython.#

#Run ``generate_qhull.py`` to regenerate the ``qhull.c`` file#cimport cythonfromlibc.float cimport DBL_EPSILONfromlibc.math cimport fabs, sqrtimportnumpy as npimportscipy.spatial.qhull as qhull

cimport scipy.spatial.qhull as qhullimportwarnings#------------------------------------------------------------------------------#Numpy etc.#------------------------------------------------------------------------------

cdef externfrom "numpy/ndarrayobject.h":

cdef enum:

NPY_MAXDIMS

ctypedef fused double_or_complex:

double

double complex#------------------------------------------------------------------------------#Interpolator base class#------------------------------------------------------------------------------

classNDInterpolatorBase(object):"""Common routines for interpolators.

.. versionadded:: 0.9"""

def __init__(self, points, values, fill_value=np.nan, ndim=None,

rescale=False, need_contiguous=True, need_values=True):"""Check shape of points and values arrays, and reshape values to

(npoints, nvalues). Ensure the `points` and values arrays are

C-contiguous, and of correct type."""

ifisinstance(points, qhull.Delaunay):#Precomputed triangulation was passed in

ifrescale:raise ValueError("Rescaling is not supported when passing"

"a Delaunay triangulation as ``points``.")

self.tri=points

points=points.pointselse:

self.tri=None

points=_ndim_coords_from_arrays(points)

values=np.asarray(values)

_check_init_shape(points, values, ndim=ndim)ifneed_contiguous:

points= np.ascontiguousarray(points, dtype=np.double)ifneed_values:

self.values_shape= values.shape[1:]if values.ndim == 1:

self.values=values[:,None]elif values.ndim == 2:

self.values=valueselse:

self.values=values.reshape(values.shape[0],

np.prod(values.shape[1:]))#Complex or real?

self.is_complex =np.issubdtype(self.values.dtype, np.complexfloating)ifself.is_complex:ifneed_contiguous:

self.values=np.ascontiguousarray(self.values,

dtype=np.complex128)

self.fill_value=complex(fill_value)else:ifneed_contiguous:

self.values= np.ascontiguousarray(self.values, dtype=np.double)

self.fill_value=float(fill_value)if notrescale:

self.scale=None

self.points=pointselse:#scale to unit cube centered at 0

self.offset = np.mean(points, axis=0)

self.points= points -self.offset

self.scale= self.points.ptp(axis=0)

self.scale[~(self.scale > 0)] = 1.0 #avoid division by 0

self.points /=self.scaledef_check_call_shape(self, xi):

xi=np.asanyarray(xi)if xi.shape[-1] != self.points.shape[1]:raise ValueError("number of dimensions in xi does not match x")returnxidef_scale_x(self, xi):if self.scale isNone:returnxielse:return (xi - self.offset) /self.scaledef __call__(self, *args):"""interpolator(xi)

Evaluate interpolator at given points.

Parameters

----------

x1, x2, ... xn: array-like of float

Points where to interpolate data at.

x1, x2, ... xn can be array-like of float with broadcastable shape.

or x1 can be array-like of float with shape ``(..., ndim)``"""xi= _ndim_coords_from_arrays(args, ndim=self.points.shape[1])

xi=self._check_call_shape(xi)

shape=xi.shape

xi= xi.reshape(-1, shape[-1])

xi= np.ascontiguousarray(xi, dtype=np.double)

xi=self._scale_x(xi)ifself.is_complex:

r=self._evaluate_complex(xi)else:

r=self._evaluate_double(xi)return np.asarray(r).reshape(shape[:-1] +self.values_shape)

cpdef _ndim_coords_from_arrays(points, ndim=None):"""Convert a tuple of coordinate arrays to a (..., ndim)-shaped array."""cdef ssize_t j, nif isinstance(points, tuple) and len(points) == 1:#handle argument tuple

points =points[0]ifisinstance(points, tuple):

p= np.broadcast_arrays(*points)

n=len(p)for j in range(1, n):if p[j].shape !=p[0].shape:raise ValueError("coordinate arrays do not have the same shape")

points= np.empty(p[0].shape + (len(points),), dtype=float)for j, item inenumerate(p):

points[...,j]=itemelse:

points=np.asanyarray(points)if points.ndim == 1:if ndim isNone:

points= points.reshape(-1, 1)else:

points= points.reshape(-1, ndim)returnpoints

cdef _check_init_shape(points, values, ndim=None):"""Check shape of points and values arrays"""

if values.shape[0] !=points.shape[0]:raise ValueError("different number of values and points")if points.ndim != 2:raise ValueError("invalid shape for input data points")if points.shape[1] < 2:raise ValueError("input data must be at least 2-D")if ndim is not None and points.shape[1] !=ndim:raise ValueError("this mode of interpolation available only for"

"%d-D data" %ndim)#------------------------------------------------------------------------------#Linear interpolation in N-D#------------------------------------------------------------------------------

classLinearNDInterpolator(NDInterpolatorBase):"""LinearNDInterpolator(points, values, fill_value=np.nan, rescale=False)

Piecewise linear interpolant in N dimensions.

.. versionadded:: 0.9

Methods

-------

__call__

Parameters

----------

points : ndarray of floats, shape (npoints, ndims); or Delaunay

Data point