ダサダサの備忘録です。
scipy.interpolate.RegularGridInterpolatorの引数やグリッドの指定の順番でいつも混乱するので、テストコードを書いてみた。
なんで混乱するかというと、scipy.interpolate.interp2dとかと、ぱっと見引数の順番が違うように見えてしまうからである。
いや、interp2dのグリッド用の引数は１次元配列二つで、RegularGridInterpolatorの第１引数は複数次元配列でしょ？ モノがちがうじゃん、と言われればそのとおりで、かつ、複数次元配列の場合、x(l個）、y(m個）、z(n個）ならshapeは(n, m, l)じゃん、と言われたら返す言葉もないのだが、、つい、逆に書いてしまうのである。＜昔から行列弱い
なので、こうやって恥をさらしておけば、次は間違わないかなー、と祈りつつ備忘録。
import numpy as np
from scipy.interpolate import RegularGridInterpolator
x_iter = np.array([0,1,2,3,4])
y_iter = np.array([5,6,7])
z_iter = np.array([8,9])
print("number of x_iter l = %d" % (len(x_iter)))
print("number of y_iter m = %d" % (len(y_iter)))
print("number of z_iter n = %d" % (len(z_iter)))
def data_func(x, y, z) :
return x + y*10 + z*100
data = []
for z in z_iter :
buf2 = []
for y in y_iter :
buf = []
for x in x_iter :
buf.append(data_func(x, y, z))
buf2.append(buf)
data.append(buf2)
data = np.array(data)
print("data", data)
print("shape of data, n x m x l array", data.shape)
data_grid = tuple([z_iter, y_iter, x_iter])
print("data grid with a format tuple(z_iter, y_iter, x_iter) is", data_grid)
f = RegularGridInterpolator(data_grid, data)
point_of_interest = tuple([8.3, 5.5, 1.1])
print("point_of_interest with a format tuple(z, y, x) is", point_of_interest)
interp_value = f(point_of_interest)
print("interpolated value is", interp_value)
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import colors
figid = 1
plt.figure(figid)
ax1 = plt.subplot(1, 1, 1)
xx, yy = np.meshgrid(x_iter, y_iter)
data_2D = data[0, :, :]
print("shape of data_2D for z_bin = 0", data_2D.shape)
phist = ax1.pcolormesh(xx, yy, data_2D, norm=colors.Normalize())
plt.colorbar(phist)
ax1.set_ylabel("y")
ax1.set_xlabel("x")
plt.savefig("interp_test.png")
forでまわしながらデータを作成する際に、回す順番に注意。
Defaultでは、x,y,zの値に、範囲外の値を指定するとエラーする。この振る舞いはbounds_error引数とfill_value変数でコントロール可能（デフォルトではbounds_error=True, fill_value=nanになっている）
ちなみに、scipy.interpolation.interp2Dを使う場合はこんな感じ。
Defaultでは、xとyの値に、範囲外の値を指定しても死なない（data table内の一番端の値を返す）。この振る舞いはbounds_error引数とfill_value変数でコントロール可能（デフォルトではbounds_error=False, fill_value=Noneになっている）
import numpy as np
from scipy.interpolate import interp2d
x_iter = np.array([0,1,2,3,4])
y_iter = np.array([5,6,7])
print("number of x_iter l = %d" % (len(x_iter)))
print("number of y_iter m = %d" % (len(y_iter)))
def data_func(x, y) :
return x + y*10
data = []
for y in y_iter :
buf = []
for x in x_iter :
buf.append(data_func(x, y))
data.append(buf)
data = np.array(data)
print("data", data)
print("shape of data, m x l array", data.shape)
# order of arguments: x-axis, y-axis, data
f = interp2d(x_iter, y_iter, data)
x = 1.5
y = 5.5
print("datapoint is (%f, %f), data_func(%f, %f) is %f" % (x, y, x, y, data_func(x, y)))
interp_value = f(x, y)
print("interpolated value is", interp_value)
# test2 : what happens if x1 is out of range?
x1 = 7 # x-range is (0, 4)
y1 = 10 # y-range is (5, 7)
print("datapoint is (%f, %f), data_func(%f, %f) is %f" % (x1, y1, x1, y1, data_func(x1, y1)))
interp_value = f(x1, y1)
print("interpolated value is", interp_value)
# this value is
print("interpolated value for (%f, %f) is %f" % (x_iter[-1], y_iter[-1], f(x_iter[-1], y_iter[-1])))
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import colors
figid = 1
plt.figure(figid)
ax1 = plt.subplot(1, 1, 1)
xx, yy = np.meshgrid(x_iter, y_iter)
phist = ax1.pcolormesh(xx, yy, data, norm=colors.Normalize())
plt.colorbar(phist)
ax1.set_ylabel("y")
ax1.set_xlabel("x")
plt.savefig("interp_test2D.png")