2次元データのグラフ
、2次元のasciiデータや画像データをnumpy arrayとして読み込むことができていれば、この2次元データを解析することができます。まずはじめに、二次元データを可視化するための手法を一つ試してみましょう。
my2Ddataという2次元のnumpy arrayがある場合、
fig, ax = plt.subplots()
myimg = ax.imshow(my2Ddata)
fig.colorbar(myimg) #カラーバーを表示
と実行すると下図のようなイメージプロットを作ることができます。
違うカラーマップを使いたい場合はimshowのcmapオプションを使います。つまり、上のコードのax.imshowのところを
myimg = ax.imshow(my2Ddata, cmap = 'plasma')
と変更すれば、違う色でプロットができます。
myimg = ax.imshow(my2Ddata, cmap = 'plasma', vmin=0, vmax=5000)
fig.colorbar(myimg)
のようにすれば、
のような画像が得られるはずです。 ここで明るい色が指定したvmaxの値、5000に対応していることに注意してください。
2次元のグラフで横軸, 縦軸を変更したい
上の画像で、ピクセルではなく、違う軸で表示したい場合があるかもしれません。例えばピクセルが光の波長に対応したり、実空間に対応したりする場合には、ピクセルの値を変換して、変換後の値で表示した方がわかりやすくなります。
ここでは、例えば1 pixel = 0.0025 µmだとします。その場合、上のコードの該当部分を
x_calibration=0.0025
y_calibration=0.0025
myimg = ax.imshow(my2Ddata, cmap = 'plasma',extent=(-0.5*y_calibration,y_calibration*my2Ddata.shape[1],-0.5*x_calibration,x_calibration*my2Ddata.shape[0]))
のように変更してextentオプションを使うことで
のように、自分の指定した軸を使ってプロットが可能です。
2次元データのフィッティング
次に、2次元のarrayを関数でフィッティングしたい場合を考えます。例えば、2次元arrayのmy2Ddataと、その2つの軸を表す1次元arrayのxdataとydataを、2次元ガウス関数でフィッティングしたいとします。 先程と同様に、まずフィッティングしたい2次元の実験データをプロットしてみましょう。
fig, ax = plt.subplots(figsize=(4.8, 3.4),dpi=200)
myimg = ax.imshow(my2Ddata)
fig.colorbar(myimg)
ここでは説明をわかりやすくするために、縦軸も横軸もデフォルトのまま、つまりデータ1つごとに目盛り1が対応しています。実際には、必要な縦軸横軸を使ってください。
次に、1変数のフィッティングで用いたcurve_fit関数は、1次元のデータしか扱えないので、この2次元データをravel()関数で無理やり1次元化します。ravel()関数の役割を説明するために、実際にmy2Ddata.ravel()をプロットしてみましょう。
fig, ax = plt.subplots(figsize=(4.8, 3.4),dpi=200)
xdata2=np.arange(0,np.size(my2Ddata.ravel()))
curve1 = ax.plot(xdata2,my2Ddata.ravel(),marker='o',label='Inline label',c='b')
すると
のようなグラフが得られます。もとのmy2Ddataとの対応関係を図示すると、
となります。つまりmy2Ddata[0,0:20]とmy2Ddata[1,0:20]と......とmy2Ddata[9,0:20]をつなげたような1次元のarrayになっています。フィッティングするためには、この1次元化したデータの横軸に対応する長さ（=全データ数=今回の場合20×10=200）のarrayを作っておく必要があります。が、もともとmy2Ddataは2次元なので、一つのデータ値につき、行番号をあらわす座標と列番号を表す座標の2つの情報が必要なので、実際に用意するのは、2×200の2次元配列です。この配列を作るために、
を使います。np.meshgrid()関数は、1次元配列2つから、格子点の座標を作り出す関数で、np.vstack()関数は2つの1次元配列をもとに2次元配列をつくる関数です。コードは以下のようになります。
xdata=np.arange(0,my2Ddata.shape[0]) #すでにある場合は作る必要なし。
ydata=np.arange(0,my2Ddata.shape[1]) #すでにある場合は作る必要なし。
meshx, meshy = np.meshgrid(xdata, ydata,indexing='ij')
xydata = np.vstack((meshx.ravel(),meshy.ravel()))
このコード、データがどの行、どの列のデータ7日をあらわすarrayを作るためのものですが、イマイチわかりにくいので、何をやっているのか図示すると
となります。
次に、フィッティングで使う二次元ガウス関数
z(x,y)=Aexp(−(x−x0)2/σ12)exp(−(y−y0)2/σ22)+O 自体を、上記の1次元化に対応できるようにしておく必要があります。つまり
def gaussian2d(xydata,x0,y0,sigma1,sigma2,amp,offset):
return amp*np.exp(-((xydata[0]-x0)**2/(2*sigma1**2)+(xydata[1]-y0)**2/(2*sigma2**2)))+ offset
となります。
ここでxydata[0]は（1次元化された）データ点の行番号を表す1次元numpy arrayが、xydata[1]は（1次元化された）データ点の列番号を表す1次元numpy arrayが、入力されることが想定されています。
この関数を使って、初期値を設定しフィッティングを行います。
guess = [5,11,3,3,2500,19]
popt, pcov = curve_fit(gaussian2d, xydata, my2Ddata.ravel(), guess)
これまでのコードをすべてまとめると
from scipy.optimize import curve_fit
def gaussian2d(xydata,x0,y0,sigma1,sigma2,amp,offset):
return amp*np.exp(-((xydata[0]-x0)**2/(2*sigma1**2)+(xydata[1]-y0)**2/(2*sigma2**2)))+ offset
xdata=np.arange(0,my2Ddata.shape[0]) #すでにある場合は作る必要なし。
ydata=np.arange(0,my2Ddata.shape[1]) #すでにある場合は作る必要なし。
meshx, meshy = np.meshgrid(xdata, ydata,indexing='ij')
xydata = np.vstack((meshx.ravel(),meshy.ravel()))
guess = [5,11,3,3,2500,19]
popt, pcov = curve_fit(gaussian2d, xydata, my2Ddata.ravel(), guess)
です。
最後に、フィッティングが適切かどうか、実際に比べてみましょう。
fig, ax = plt.subplots(figsize=(4.8, 3.4),dpi=200)
curve1 = ax.plot(my2Ddata.ravel(),marker='None',label='Inline label',c='r')
curve2 = ax.plot(gaussian2d(xydata,*popt),marker='o',c='b',ls='None')
結果は、
のようになり、たしかにフィッティングは合っていそうですが、1次元化した後での比較なので、直感的によくわかりません。 そこで、フィッティング結果をもとに2次元のイメージプロットを書いてみましょう。 そのためには、reshape関数を使って、フィッティング結果を2次元にreshapeしたうえでプロットをします。つまり
fig, ax = plt.subplots(figsize=(4.8, 3.4),dpi=200)
myimg = ax.imshow(gaussian2d(xydata,*popt).reshape(my2Ddata.shape))
fig.colorbar(myimg)
とすると、結果として
が得られて、フィッティングが良さそうであることがわかります。より定量的にフィッティングの正当性をチェックしたい場合には、適当な場所で断面をとってみます。例えば
fig, ax = plt.subplots(figsize=(4.8, 3.4),dpi=200)
ax.tick_params(direction="in")
ax.tick_params(bottom=True, top=True, left=True, right=True)
curve1 = ax.plot(my2Ddata[5],marker='None',label='Inline label',c='r')
curve2 = ax.plot(gaussian2d(xydata,*popt).reshape(my2Ddata.shape)[5],marker='o',c='b',ls='None')
とすれば、
のように、かなりよく、フィッティングできていることがわかります。
