新增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 ?