python - SciPy 中的二维插值问题，非矩形网格

Question

我一直在尝试使用 scipy.interpolate.bisplrep() 和 scipy.interpolate.interp2d() 在我的 (218x135) 2D 球形极坐标网格上查找数据的插值。我将我的网格节点的笛卡尔位置的二维数组 X 和 Y 传递给它们。我不断收到如下错误(对于线性插值。使用 interp2d):

“警告:不能添加更多的结，因为额外的结会重合 和一个旧的。可能原因:重量太小或太大 到一个不准确的数据点。 (fp>s) kx,ky=1,1 nx,ny=4,5 m=29430 fp=1390609718.902140 s=0.000000"

使用默认的平滑参数 s 等，我得到了双变量样条的类似结果。我的数据是平滑的。我在下面附上了我的代码，以防我做错了什么。

有什么想法吗？ 谢谢! 凯尔

class Field(object):
  Nr = 0
  Ntheta = 0
  grid = np.array([])

  def __init__(self, Nr, Ntheta, f):
    self.Nr = Nr
    self.Ntheta = Ntheta
    self.grid = np.empty([Nr, Ntheta])
    for i in range(Nr):
      for j in range(Ntheta):
        self.grid[i,j] = f[i*Ntheta + j]


def calculate_lines(filename):
  ri,ti,r,t,Br,Bt,Bphi,Bmag = np.loadtxt(filename, skiprows=3,\
    usecols=(1,2,3,4,5,6,7,9), unpack=True)
  Nr = int(max(ri)) + 1
  Ntheta = int(max(ti)) + 1

  ### Initialise coordinate grids ###
  X = np.empty([Nr, Ntheta])
  Y = np.empty([Nr, Ntheta])
  for i in range(Nr):
    for j in range(Ntheta):
      indx = i*Ntheta + j
      X[i,j] = r[indx]*sin(t[indx])
      Y[i,j] = r[indx]*cos(t[indx])

  ### Initialise field objects ###
  Bradial = Field(Nr=Nr, Ntheta=Ntheta, f=Br)

  ### Interpolate the fields ###
  intp_Br = interpolate.interp2d(X, Y, Bradial.grid, kind='linear')

  #rbf_0 = interpolate.Rbf(X,Y, Bradial.grid, epsilon=2)

  return

Answer 1

8 月 27 日添加:Kyle 在 scipy-user thread .

8 月 30 日:@Kyle，Cartesion X,Y 和极地 Xnew,Ynew 之间似乎存在混淆。 请参阅下面太长的注释中的“极地”。

# griddata vs SmoothBivariateSpline
# http://stackoverflow.com/questions/3526514/
#   problem-with-2d-interpolation-in-scipy-non-rectangular-grid

# http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
# http://en.wikipedia.org/wiki/Natural_neighbor
# http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

from __future__ import division
import sys
import numpy as np
from scipy.interpolate import SmoothBivariateSpline  # $scipy/interpolate/fitpack2.py
from matplotlib.mlab import griddata

__date__ = "2010-10-08 Oct"  # plot diffs, ypow
    # "2010-09-13 Sep"  # smooth relative

def avminmax( X ):
    absx = np.abs( X[ - np.isnan(X) ])
    av = np.mean(absx)
    m, M = np.nanmin(X), np.nanmax(X)
    histo = np.histogram( X, bins=5, range=(m,M) ) [0]
    return "av %.2g  min %.2g  max %.2g  histo %s" % (av, m, M, histo)

def cosr( x, y ):
    return 10 * np.cos( np.hypot(x,y) / np.sqrt(2) * 2*np.pi * cycle )

def cosx( x, y ):
    return 10 * np.cos( x * 2*np.pi * cycle )

def dipole( x, y ):
    r = .1 + np.hypot( x, y )
    t = np.arctan2( y, x )
    return np.cos(t) / r**3

#...............................................................................
testfunc = cosx
Nx = Ny = 20  # interpolate random Nx x Ny points -> Newx x Newy grid
Newx = Newy = 100
cycle = 3
noise = 0
ypow = 2  # denser => smaller error
imclip = (-5., 5.)  # plot trierr, splineerr to same scale
kx = ky = 3
smooth = .01  # Spline s = smooth * z2sum, see note
    # s is a target for sum (Z() - spline())**2  ~ Ndata and Z**2;
    # smooth is relative, s absolute
    # s too small => interpolate/fitpack2.py:580: UserWarning: ier=988, junk out
    # grr error message once only per ipython session
seed = 1
plot = 0

exec "\n".join( sys.argv[1:] )  # run this.py N= ...
np.random.seed(seed)
np.set_printoptions( 1, threshold=100, suppress=True )  # .1f

print 80 * "-"
print "%s  Nx %d Ny %d -> Newx %d Newy %d  cycle %.2g noise %.2g  kx %d ky %d smooth %s" % (
    testfunc.__name__, Nx, Ny, Newx, Newy, cycle, noise, kx, ky, smooth)

#...............................................................................

    # interpolate X Y Z to xnew x ynew --
X, Y = np.random.uniform( size=(Nx*Ny, 2) ) .T
Y **= ypow
    # 1d xlin ylin -> 2d X Y Z, Ny x Nx --
    # xlin = np.linspace( 0, 1, Nx )
    # ylin = np.linspace( 0, 1, Ny )
    # X, Y = np.meshgrid( xlin, ylin )
Z = testfunc( X, Y )  # Ny x Nx
if noise:
    Z += np.random.normal( 0, noise, Z.shape )
# print "Z:\n", Z
z2sum = np.sum( Z**2 )

xnew = np.linspace( 0, 1, Newx )
ynew = np.linspace( 0, 1, Newy )
Zexact = testfunc( *np.meshgrid( xnew, ynew ))
if imclip is None:
    imclip = np.min(Zexact), np.max(Zexact)
xflat, yflat, zflat = X.flatten(), Y.flatten(), Z.flatten()

#...............................................................................
print "SmoothBivariateSpline:"
fit = SmoothBivariateSpline( xflat, yflat, zflat, kx=kx, ky=ky, s = smooth * z2sum )
Zspline = fit( xnew, ynew ) .T  # .T ??

splineerr = Zspline - Zexact
print "Zspline - Z:", avminmax(splineerr)
print "Zspline:    ", avminmax(Zspline)
print "Z:          ", avminmax(Zexact)
res = fit.get_residual()
print "residual %.0f  res/z2sum %.2g" % (res, res / z2sum)
# print "knots:", fit.get_knots()
# print "Zspline:", Zspline.shape, "\n", Zspline
print ""

#...............................................................................
print "griddata:"
Ztri = griddata( xflat, yflat, zflat, xnew, ynew )
        # 1d x y z -> 2d Ztri on meshgrid(xnew,ynew)

nmask = np.ma.count_masked(Ztri)
if nmask > 0:
    print "info: griddata: %d of %d points are masked, not interpolated" % (
        nmask, Ztri.size)
    Ztri = Ztri.data  # Nans outside convex hull
trierr = Ztri - Zexact
print "Ztri - Z:", avminmax(trierr)
print "Ztri:    ", avminmax(Ztri)
print "Z:       ", avminmax(Zexact)
print ""

#...............................................................................
if plot:
    import pylab as pl
    nplot = 2
    fig = pl.figure( figsize=(10, 10/nplot + .5) )
    pl.suptitle( "Interpolation error: griddata - %s, BivariateSpline - %s" % (
        testfunc.__name__, testfunc.__name__ ), fontsize=11 )

    def subplot( z, jplot, label ):
        ax = pl.subplot( 1, nplot, jplot )
        im = pl.imshow(
            np.clip( z, *imclip ),  # plot to same scale
            cmap=pl.cm.RdYlBu,
            interpolation="nearest" )
                # nearest: squares, else imshow interpolates too
                # todo: centre the pixels
        ny, nx = z.shape
        pl.scatter( X*nx, Y*ny, edgecolor="y", s=1 )  # for random XY
        pl.xlabel(label)
        return [ax, im]

    subplot( trierr, 1,
        "griddata, Delaunay triangulation + Natural neighbor: max %.2g" %
        np.nanmax(np.abs(trierr)) )

    ax, im = subplot( splineerr, 2,
        "SmoothBivariateSpline kx %d ky %d smooth %.3g: max %.2g" % (
        kx, ky, smooth, np.nanmax(np.abs(splineerr)) ))

    pl.subplots_adjust( .02, .01, .92, .98, .05, .05 )  # l b r t
    cax = pl.axes([.95, .05, .02, .9])  # l b w h
    pl.colorbar( im, cax=cax )  # -1.5 .. 9 ??
    if plot >= 2:
        pl.savefig( "tmp.png" )
    pl.show() 

关于二维插值、BivariateSpline 与网格数据的说明。

scipy.interpolate.*BivariateSpline 和 matplotlib.mlab.griddata 都将一维数组作为参数:

Znew = griddata( X,Y,Z, Xnew,Ynew )
    # 1d X Y Z Xnew Ynew -> interpolated 2d Znew on meshgrid(Xnew,Ynew)
assert X.ndim == Y.ndim == Z.ndim == 1  and  len(X) == len(Y) == len(Z)

输入 X,Y,Z 描述 3 空间中的表面或点云: X,Y(或纬度、经度或...)平面中的点， 和 Z 上面的表面或地形。 X,Y 可能会填满矩形 [Xmin .. Xmax] x [Ymin .. Ymax] 的大部分， 或者可能只是其中的波浪形 S 或 Y。 Z 表面可能是光滑的，或光滑 + 一点噪音， 或者根本不平坦，粗糙的火山山脉。

Xnew和Ynew通常也是1d，描述一个矩形网格 的 |Xnew| x |Y新|您想要插值或估计 Z 的点。
Znew = griddata(...) 在这个网格上返回一个二维数组，np.meshgrid(Xnew,Ynew):

Znew[Xnew0,Ynew0], Znew[Xnew1,Ynew0], Znew[Xnew2,Ynew0] ...
Znew[Xnew0,Ynew1] ...
Znew[Xnew0,Ynew2] ...
...

Xnew,Ynew 点远离任何输入X,Y 的拼写麻烦。 griddata 检查这个:

A masked array is returned if any grid points are outside convex hull defined by input data (no extrapolation is done).

("Convex hull"是一个假想的里面的区域 橡皮筋围绕所有 X、Y 点拉伸(stretch)。)

griddata 首先构建一个 Delaunay 三角剖分 输入 X,Y，然后做 Natural neighbor 插值。这是稳健且非常快速的。

不过，BivariateSpline 可以外推， 在没有警告的情况下产生剧烈波动。 此外，Fitpack 中的所有 *Spline 例程 对平滑参数 S 非常敏感。 Dierckx 的书(books.google isbn 019853440X p. 89)说:
如果 S 太小，则样条逼近太摇摆不定 并拾取太多噪音(过拟合)；
如果 S 太大，样条会太平滑 并且信号将丢失(欠拟合)。

散乱数据的插值很难，平滑不容易，两者结合起来真的很难。 插值器应该如何处理 XY 中的大孔或非常嘈杂的 Z ？ (“如果你想出售它，你将不得不描述它。”)

还有更多注意事项，细则:

1d 与 2d:一些插值器将 X、Y、Z 取为 1d 或 2d。 其他人只取 1d，所以在插值之前展平:

Xmesh, Ymesh = np.meshgrid( np.linspace(0,1,Nx), np.linspace(0,1,Ny) )
Z = f( Xmesh, Ymesh )  # Nx x Ny
Znew = griddata( Xmesh.flatten(), Ymesh.flatten(), Z.flatten(), Xnew, Ynew )

在屏蔽数组上:matplotlib 处理得很好， 仅绘制未屏蔽/非 NaN 点。 但我不敢打赌笨蛋 numpy/scipy 函数会起作用。 检查 X、Y 的凸包外的插值，如下所示:

Znew = griddata(...)
nmask = np.ma.count_masked(Znew)
if nmask > 0:
    print "info: griddata: %d of %d points are masked, not interpolated" % (
        nmask, Znew.size)
    # Znew = Znew.data  # array with NaNs

关于极坐标: X,Y和Xnew,Ynew应该在同一个空间， 两个 Cartesion，或两者都在 [rmin .. rmax] x [tmin .. tmax] 中。
在 3d 中绘制 (r, theta, z) 点:

from mpl_toolkits.mplot3d import Axes3D
Znew = griddata( R,T,Z, Rnew,Tnew )
ax = Axes3D(fig)
ax.plot_surface( Rnew * np.cos(Tnew), Rnew * np.sin(Tnew), Znew )

另请参阅(还没有尝试过):

ax = subplot(1,1,1, projection="polar", aspect=1.)
ax.pcolormesh(theta, r, Z)

---

给谨慎的程序员的两个提示:

检查异常值或有趣的缩放比例:

def minavmax( X ):
    m = np.nanmin(X)
    M = np.nanmax(X)
    av = np.mean( X[ - np.isnan(X) ])  # masked ?
    histo = np.histogram( X, bins=5, range=(m,M) ) [0]
    return "min %.2g  av %.2g  max %.2g  histo %s" % (m, av, M, histo)

for nm, x in zip( "X Y Z  Xnew Ynew Znew".split(),
                (X,Y,Z, Xnew,Ynew,Znew) ):
    print nm, minavmax(x)

检查简单数据的插值:

interpolate( X,Y,Z, X,Y )  -- interpolate at the same points
interpolate( X,Y, np.ones(len(X)), Xnew,Ynew )  -- constant 1 ?