新增27Aug：凯尔跟着这件事上 scipy-user thread。
30Aug：@凯尔，看起来好像Cartesion X，Y和极地Xnew，Ynew之间有混淆。 请参阅下面的过长说明中的“polar”。
# 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对比的GridData。
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描述在三维空间的点的表面或云： X,Y（或纬度，经度或...）分在飞机上， 和Z上面的表面或地形。 X,Y可以填充矩形的大部分[Xmin .. Xmax] x [Ymin .. Ymax]， 或者可能只是它内部的一个波浪状的S或Y. 该Z表面可能是光滑的，或光滑+有点噪音， 或不光滑，粗糙的火山山脉。
Xnew和Ynew通常也是1d，描述了一个矩形网格 | Xnew | x | Ynew |
Xnew，Ynew点远：要插值点或估计Z.
Znew =的GridData（...）在这个格，np.meshgrid（Xnew，Ynew）返回一个二维数组任何输入X，Y的拼写都有问题。 griddata检查此：
如果任何网格点以外凸 船体由输入数据所定义（没有外插完成）甲屏蔽数组被返回。
（“凸包”是区域的假想 橡皮筋拉伸周围所有的X，Y点。内）
griddata作品通过首先构建Delaunay三角网的输入X，Y的 ，然后做 Natural neighbor 插值。这是强大和快速的。
但是，BivariateSpline可以推断出 产生大幅波动而不会发出警告。 此外，所有在* Fitpack 样条程序是平滑参数S. Dierckx的书很敏感的说（books.google ISBN 019853440X第89页。）：
如果S是太小了，样条近似太蠕动 和拾起太多噪音（过度配合）;
如果S太大，则样条曲线太平滑 并且信号将丢失（欠调）。
分散数据的插值很难，平滑不容易，两者在一起真的很难。 插值器应该在XY中用大孔还是用非常嘈杂的Z做什么？ （“如果你想卖掉它，你将不得不来形容它”）
然而，更多的音符，小字：
1D VS 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的点。 但我不敢肯定，一个bozo 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应该是在同一个空间， 两个笛卡尔，或者两者[RMIN .. rmax] x [tmin .. tmax]。
要在3D绘图（R，θ，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 ?