用于 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[0], y = y0[0];当 x > x0[-1] 时，y = y0[-1]。这里，假设 x0 是升序排列。

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

作为引用，我当前的 extrapolated_spline_2D 版本非常慢。基本思路是:

(1)首先，给定一维数组x0、y0和二维数组z2d0作为输入，使nx0个extrapolated_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

我做过类似的工作，叫做 GlobalSpline2D here ，并且它在线性、三次或五次样条下都能完美地工作。

基本上是继承了interp2d ，并通过 InterpolatedUnivariateSpline 促进二维外推的使用.它们都是 scipy 内部函数。

其用法引用document以及 call method interp2d 的。