使用LinearNDInterpolator进行外推

Question

我有一个3D数据集，希望能够进行线性插值和外推。可以使用scipy.interpolate.LinearNDInterpolator轻松进行插值。该模块只能在参数范围之外的值上填充常数/NaN，但我不知道为什么它不提供打开外推选项的选项。

看代码时，我发现该模块是用cython编写的。由于没有cython经验，很难通过代码来实现外推。我可以用纯python代码编写它，但也许这里还有其他更好的想法？我的特殊情况涉及一个恒定的xy网格，但z值变化很大（-100,000），因此插值必须快速进行，因为每次z值改变时都会运行插值。

为了举例说明，假设我有一个像这样的网格：

xyPairs = [[-1.0, 0.0], [-1.0, 4.0],
           [-0.5, 0.0], [-0.5, 4.0],
           [-0.3, 0.0], [-0.3, 4.0],
           [+0.0, 0.0], [+0.0, 4.0],
           [+0.2, 0.0], [+0.2, 4.0]]

假设我想在x = -1.5, -0.8, +0.5和y = -0.2, +0.2, +0.5处计算数值。目前，我正在对每个y值沿着x轴执行1d插值/外推，然后对每个x值沿着y轴执行插值/外推。外推由ryggyr's answer中的第二个函数完成。

- oschoudhury

---

哦，一个月前问过了，现在编辑一下 - 你能发一些代码吗？比如你的数据集长什么样子，你现在使用什么，它给出什么结果，以及你想让它看起来像什么？ - usethedeathstar

Answer 1

使用最近邻和线性插值的组合。 如果LinearNDInterpolator无法插值，则返回np.nan， 否则它将返回大小为1的数组。 NearestNDInterpolator返回一个浮点数。

import scipy.interpolate
import numpy
class LinearNDInterpolatorExt(object):
  def __init__(self, points,values):
    self.funcinterp=scipy.interpolate.LinearNDInterpolator(points,values)
    self.funcnearest=scipy.interpolate.NearestNDInterpolator(points,values)
  def __call__(self,*args):
    t=self.funcinterp(*args)
    if not numpy.isnan(t):
      return t.item(0)
    else:
      return self.funcnearest(*args)

- Keith Williams

---

我喜欢这个想法，但是你的__call__方法有点像全或无解决方案。我的猜测是OP想要进行插值，尽可能填充结果边缘NaN值与最近的值。不过，这仍然是一个很好的起点。 - medley56

Answer 2

我提出一种方法，代码很糟糕，但我希望它能对你有所帮助。这个想法是，如果你预先知道需要外推的范围，就可以在数组的边缘添加线性外推值的额外列/行，然后在新数组上进行插值。以下是一个示例，其中包含一些数据，将进行外推直到x = ±50和y = ±40：

import numpy as np
x,y=np.meshgrid(np.linspace(0,6,7),np.linspace(0,8,9)) # create x,y grid
z=x**2*y # and z values
# create larger versions with two more columns/rows
xlarge=np.zeros((x.shape[0]+2,x.shape[1]+2))
ylarge=np.zeros((x.shape[0]+2,x.shape[1]+2))
zlarge=np.zeros((x.shape[0]+2,x.shape[1]+2))
xlarge[1:-1,1:-1]=x # copy data on centre
ylarge[1:-1,1:-1]=y
zlarge[1:-1,1:-1]=z
# fill extra columns/rows
xmin,xmax=-50,50
ymin,ymax=-40,40
xlarge[:,0]=xmin;xlarge[:,-1]=xmax # fill first/last column
xlarge[0,:]=xlarge[1,:];xlarge[-1,:]=xlarge[-2,:] # copy first/last row
ylarge[0,:]=ymin;ylarge[-1,:]=ymax
ylarge[:,0]=ylarge[:,1];ylarge[:,-1]=ylarge[:,-2]
# for speed gain: store factor of first/last column/row
first_column_factor=(xlarge[:,0]-xlarge[:,1])/(xlarge[:,1]-xlarge[:,2]) 
last_column_factor=(xlarge[:,-1]-xlarge[:,-2])/(xlarge[:,-2]-xlarge[:,-3])
first_row_factor=(ylarge[0,:]-ylarge[1,:])/(ylarge[1,:]-ylarge[2,:])
last_row_factor=(ylarge[-1,:]-ylarge[-2,:])/(ylarge[-2,:]-ylarge[-3,:])
# extrapolate z; this operation only needs to be repeated when zlarge[1:-1,1:-1] is updated
zlarge[:,0]=zlarge[:,1]+first_column_factor*(zlarge[:,1]-zlarge[:,2]) # extrapolate first column
zlarge[:,-1]=zlarge[:,-2]+last_column_factor*(zlarge[:,-2]-zlarge[:,-3]) # extrapolate last column
zlarge[0,:]=zlarge[1,:]+first_row_factor*(zlarge[1,:]-zlarge[2,:]) # extrapolate first row
zlarge[-1,:]=zlarge[-2,:]+last_row_factor*(zlarge[-2,:]-zlarge[-3,:]) #extrapolate last row

那么你可以在（xlarge，ylarge，zlarge）上进行插值。由于所有操作都是numpy切片操作，所以我希望它对您来说足够快。当z数据更新时，请将它们复制到zlarge [1：-1，1：-1]中并重新执行最后4行。

- JPG

---

哎呀，它看起来很丑，但确实很快。 - oschoudhury

不适用于N维插值和外推。 - denfromufa

Answer 3

谢谢@Keith William!

进行一些小修改，也可以使其与np.ndarray一起使用:

from typing import Union

import numpy as np
import scipy.interpolate


class LinearNDInterpolatorExtrapolator:

    def __init__(self, points: np.ndarray, values: np.ndarray):
        """ Use ND-linear interpolation over the convex hull of points, and nearest neighbor outside (for
            extrapolation)

            Idea taken from https://dev59.com/FHrZa4cB1Zd3GeqPyiBd
        """
        self.linear_interpolator = scipy.interpolate.LinearNDInterpolator(points, values)
        self.nearest_neighbor_interpolator = scipy.interpolate.NearestNDInterpolator(points, values)

    def __call__(self, *args) -> Union[float, np.ndarray]:
        t = self.linear_interpolator(*args)
        t[np.isnan(t)] = self.nearest_neighbor_interpolator(*args)[np.isnan(t)]
        if t.size == 1:
            return t.item(0)
        return t

- GabCaz

Answer 4

我稍微修改了@Keith Williams的答案，它对我很有效（请注意它不是线性外推 - 它只使用最近邻）：

import numpy as np
from scipy.interpolate import LinearNDInterpolator as linterp
from scipy.interpolate import NearestNDInterpolator as nearest

class LinearNDInterpolatorExt(object):
    def __init__(self, points, values):
        self.funcinterp = linterp(points, values)
        self.funcnearest = nearest(points, values)
    
    def __call__(self, *args):
        z = self.funcinterp(*args)
        chk = np.isnan(z)
        if chk.any():
            return np.where(chk, self.funcnearest(*args), z)
        else:
            return z

- Christo