SciPy，非矩形网格中的2D插值问题 [英] Problem with 2D interpolation in SciPy, non-rectangular grid

问题描述

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

I've been trying to use scipy.interpolate.bisplrep() and scipy.interpolate.interp2d() to find interpolants for data on my (218x135) 2D spherical-polar grid. To these I pass 2D arrays, X and Y, of the Cartesian positions of my grid nodes. I keep getting errors like the following (for linear interp. with interp2d):

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

"Warning: No more knots can be added because the additional knot would coincide with an old one. Probably cause: s too small or too large a weight to an inaccurate data point. (fp>s) kx,ky=1,1 nx,ny=4,5 m=29430 fp=1390609718.902140 s=0.000000"

对于带有平滑参数s等默认值的双变量样条曲线，我得到了类似的结果.我的数据是平滑的.我在下面附加了我的代码，以防万一我做错了明显的事情.

I get a similar result for bivariate splines with the default value of the smoothing parameter s etc. My data are smooth. I've attached my code below in case I'm doing something obviously wrong.

有什么想法吗? 谢谢！ 凯尔

Any ideas? Thanks! Kyle

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

推荐答案

添加了27Aug:Kyle在 密码用户线程.

Added 27Aug: Kyle followed this up on a scipy-user thread.

8月30日:@Kyle，在Cartesion X，Y和Polar Xnew，Ynew之间似乎存在混淆. 请参见下面太长的注释中的极性".

30Aug: @Kyle, it looks as though there's a mixup between Cartesion X,Y and polar Xnew,Ynew. See "polar" in the too-long notes below.

# 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() 

关于2d插值，BivariateSpline与griddata的说明.

Notes on 2d interpolation, BivariateSpline vs. griddata.

scipy.interpolate.*BivariateSpline和matplotlib.mlab.griddata 两者都以1d数组作为参数:

scipy.interpolate.*BivariateSpline and matplotlib.mlab.griddata both take 1d arrays as arguments:

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表面可能是光滑的，也可能是光滑的+有点噪音， 或根本不平坦的火山山.

The inputs X,Y,Z describe a surface or cloud of points in 3-space: X,Y (or latitude,longitude or ...) points in a plane, and Z a surface or terrain above that. X,Y may fill most of the rectangle [Xmin .. Xmax] x [Ymin .. Ymax], or may be just a squiggly S or Y inside it. The Z surface may be smooth, or smooth + a bit of noise, or not smooth at all, rough volcanic mountains.

Xnew和Ynew通常也是1d，描述一个矩形网格 的| Xnew | x | Ynew |要插入或估算Z的点.
Znew = griddata(...)返回此网格上的二维数组np.meshgrid(Xnew，Ynew):

Xnew and Ynew are usually also 1d, describing a rectangular grid of |Xnew| x |Ynew| points where you want to interpolate or estimate Z.
Znew = griddata(...) returns a 2d array over this grid, np.meshgrid(Xnew,Ynew):

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

Xnew，Ynew指向远离任何输入X，Y的拼写问题. griddata对此进行检查:

Xnew,Ynew points far from any of the input X,Y s spell trouble. griddata checks this:

如果任何网格点在凸面之外，则返回蒙版数组 输入数据定义的船体(不进行推断).

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

(凸包"是虚构内的区域 橡皮筋围绕所有X，Y点伸展.)

("Convex hull" is the area inside an imaginary rubber band stretched around all the X,Y points.)

griddata通过首先构造Delaunay三角剖分来工作 X，Y的输入，然后做 自然邻居 插值.这是健壮且相当快的.

griddata works by first constructing a Delaunay triangulation of the input X,Y, then doing Natural neighbor interpolation. This is robust and quite fast.

BivariateSpline可以外推， 生成没有警告的疯狂秋千. 此外，Fitpack中的所有* Spline例程 对平滑参数S非常敏感 Dierckx的书(books.google isbn 019853440X第89页)说:
如果S太小，则样条近似值太摆动 并吸收了过多的噪音(过拟合)；
如果S太大，则样条曲线将太平滑 并且信号将丢失(欠拟合).

BivariateSpline, though, can extrapolate, generating wild swings without warning. Furthermore, all the *Spline routines in Fitpack are very sensitive to smoothing parameter S. Dierckx's book (books.google isbn 019853440X p. 89) says:
if S is too small, the spline approximation is too wiggly and picks up too much noise (overfit);
if S is too large the spline will be too smooth and signal will be lost (underfit).

对分散的数据进行插值很困难，平滑起来并不容易，两者加在一起确实很困难. 内插器在XY中有大孔或Z噪声很大时应如何处理? (如果要出售它，则必须对其进行描述.")

Interpolation of scattered data is hard, smoothing not easy, both together really hard. What should an interpolator do with big holes in XY, or with very noisy Z ? ("If you want to sell it, you're going to have to describe it.")

还有更多笔记，印刷精美:

Yet more notes, fine print:

1d与2d:某些内插器将X，Y，Z取为1d或2d. 其他人只取1d，因此在插值之前先展平:

1d vs 2d: Some interpolators take X,Y,Z either 1d or 2d. Others take 1d only, so flatten before interpolating:

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的点. 但是我不会打赌bozo numpy/scipy函数将完全起作用. 像这样检查X，Y的凸包之外的插值:

On masked arrays: matplotlib handles them just fine, plotting only unmasked / non-NaN points. But I wouldn't bet that that a bozo numpy/scipy functions would work at all. Check for interpolation outside the convex hull of X,Y like this:

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应该在同一空间中， 两种笛卡尔，或两者都以[rmin .. rmax] x [tmin .. tmax]表示.
在3d中绘制(r，theta，z)点:

On polar coordinates: X,Y and Xnew,Ynew should be in the same space, both Cartesion, or both in [rmin .. rmax] x [tmin .. tmax].
To plot (r, theta, z) points in 3d:

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 )

另请参阅(尚未尝试过):

See also (haven't tried this):

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

---

警惕的程序员的两个提示:

---

Two tips for the wary programmer:

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

check for outliers, or funny scaling:

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)

使用简单数据检查插值:

check interpolation with simple data:

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

这篇关于SciPy，非矩形网格中的2D插值问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！