用于 2D 外推样条函数的 Python Scipy？

Question

我想为二维矩阵编写一个外推样条函数。我现在拥有的是一维数组的外推样条函数，如下所示。使用scipy.interpolate.InterpolatedUnivariateSpline()。

import numpy as np 
import scipy as sp 

def extrapolated_spline_1D(x0,y0):
    x0 = np.array(x0)
    y0 = np.array(y0)
    assert x0.shape == y.shape 

    spline = sp.interpolate.InterpolatedUnivariateSpline(x0,y0)
    def f(x, spline=spline):
        return np.select(
            [(x<x0[0]),              (x>x0[-1]),              np.ones_like(x,dtype='bool')], 
            [np.zeros_like(x)+y0[0], np.zeros_like(x)+y0[-1], spline(x)])

    return f

它需要 x0，这是定义函数的位置，以及 y0，这是相应的值。当 x x0[-1] 时，y = y0[-1]。这里，假设 x0 是按升序排列的。

我想要一个类似的外插样条函数来处理使用np.select() 的二维矩阵，就像在 extrapolated_spline_1D 中一样。我认为scipy.interpolate.RectBivariateSpline() 可能会有所帮助，但我不知道该怎么做。

作为参考，我当前的 extrapolated_spline_2D 版本非常慢。 基本思路是：

(1)首先，给定一维数组x0，y0和二维数组z2d0作为输入，制作nx0extrapolated_spline_1D函数，y0_spls，每一个代表定义在y0上的一层z2d0；

(2) 第二，对于不在网格上的点(x,y)，计算nx0个值，每个等于y0_spls[i](y)；

(3) 第三个，用extrapolated_spline_1D 拟合(x0, y0_spls[i](y)) 到x_spl 并返回x_spl(x) 作为最终结果。

def extrapolated_spline_2D(x0,y0,z2d0): 
    '''    
    x0,y0 : array_like, 1-D arrays of coordinates in strictly monotonic order. 
    z2d0  : array_like, 2-D array of data with shape (x.size,y.size).
    '''    
    nx0 = x0.shape[0]
    ny0 = y0.shape[0]
    assert z2d0.shape == (nx0,ny0)

    # make nx0 splines, each of which stands for a layer of z2d0 on y0 
    y0_spls = [extrapolated_spline_1D(y0,z2d0[i,:]) for i in range(nx0)]

    def f(x, y):     
        '''
        f takes 2 arguments at the same time --> x, y have the same dimention
        Return: a numpy ndarray object with the same shape of x and y
        '''
        x = np.array(x,dtype='f4')
        y = np.array(y,dtype='f4') 
        assert x.shape == y.shape        
        ndim = x.ndim 

        if ndim == 0:    
            '''
            Given a point on the xy-plane. 
            Make ny = 1 splines, each of which stands for a layer of new_xs on x0
            ''' 
            new_xs = np.array([y0_spls[i](y) for i in range(nx0)]) 
            x_spl  = extrapolated_spline_1D(x0,new_xs)
            result = x_spl(x)

        elif ndim == 1:
            '''
            Given a 1-D array of points on the xy-plane. 
            '''
            ny     = len(y)            
            new_xs = np.array([y0_spls[i](y)                 for i in range(nx0)]) # new_xs.shape = (nx0,ny)       
            x_spls = [extrapolated_spline_1D(x0,new_xs[:,i]) for i in range(ny)]
            result = np.array([x_spls[i](x[i])               for i in range(ny)])

        else:
            '''
            Given a multiple dimensional array of points on the xy-plane.  
            '''
            x_flatten = x.flatten()
            y_flatten = y.flatten()     
            ny = len(y_flatten)       
            new_xs = np.array([y0_spls[i](y_flatten)         for i in range(nx0)])         
            x_spls = [extrapolated_spline_1D(x0,new_xs[:,i]) for i in range(ny)]
            result = np.array([x_spls[i](x_flatten[i])       for i in range(ny)]).reshape(y.shape)
        return result      
    return f

Answer 1

我想我自己想出了一个答案，它使用scipy.interpolate.RectBivariateSpline()，比我的旧答案快 10 倍以上。

这是函数extrapolated_spline_2D_new。

def extrapolated_spline_2D_new(x0,y0,z2d0):
    '''    
    x0,y0 : array_like，1-D arrays of coordinates in strictly ascending order. 
    z2d0  : array_like，2-D array of data with shape (x.size,y.size).
    ''' 
    assert z2d0.shape == (x0.shape[0],y0.shape[0])

    spline = scipy.interpolate.RectBivariateSpline(x0,y0,z2d0,kx=3,ky=3)
    '''
    scipy.interpolate.RectBivariateSpline
    x,y : array_like, 1-D arrays of coordinates in strictly ascending order.
    z   : array_like, 2-D array of data with shape (x.size,y.size).
    '''  
    def f(x,y,spline=spline):
        '''
        x and y have the same shape with the output. 
        '''
        x = np.array(x,dtype='f4')
        y = np.array(y,dtype='f4') 
        assert x.shape == y.shape 
        ndim = x.ndim   
        # We want the output to have the same dimension as the input, 
        # and when ndim == 0 or 1, spline(x,y) is always 2D. 
        if   ndim == 0: result = spline(x,y)[0][0]
        elif ndim == 1: 
            result = np.array([spline(x[i],y[i])[0][0] for i in range(len(x))])
        else:           
            result = np.array([spline(x.flatten()[i],y.flatten()[i])[0][0] for i in range(len(x.flatten()))]).reshape(x.shape)         
        return result
    return f

注意： 在上面的版本中，我是一一计算值，而不是使用下面的代码。

def f(x,y,spline=spline):
    '''
    x and y have the same shape with the output. 
    '''
    x = np.array(x,dtype='f4')
    y = np.array(y,dtype='f4') 
    assert x.shape == y.shape 
    ndim = x.ndim
    if   ndim == 0: result = spline(x,y)[0][0]
    elif ndim == 1: 
         result = spline(x,y).diagonal()
    else:           
         result = spline(x.flatten(),y.flatten()).diagonal().reshape(x.shape)       
    return result

因为当我尝试使用下面的代码进行计算时，它有时会给出错误消息：

<ipython-input-65-33285fd2319d> in f(x, y, spline)
 29         if   ndim == 0: result = spline(x,y)[0][0]
 30         elif ndim == 1:
---> 31             result = spline(x,y).diagonal()
 32         else:
 33             result = spline(x.flatten(),y.flatten()).diagonal().reshape(x.shape)

/usr/local/lib/python2.7/site-packages/scipy/interpolate/fitpack2.pyc in __call__(self, x, y, mth, dx, dy, grid)
826                 z,ier = dfitpack.bispev(tx,ty,c,kx,ky,x,y)
827                 if not ier == 0:
--> 828                     raise ValueError("Error code returned by bispev: %s" % ier)
829         else:
830             # standard Numpy broadcasting

ValueError: Error code returned by bispev: 10

我不知道这是什么意思。

Answer 2

我已经完成了一项名为 GlobalSpline2D here 的类似工作，它在线性样条、三次样条或五次样条下都能完美运行。

基本上它继承了interp2d，并通过InterpolatedUnivariateSpline 将其推广到2D 外推。它们都是 scipy 内部函数。

其用法参考interp2d的document和call method。