2次元データの解析（20211221更新, 書きかけ）

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')

と変更すれば、違う色でプロットができます。

のようにすれば変更可能です。使えるキーワードは https://matplotlib.org/stable/gallery/color/colormap_reference.html を参照のこと。 また、同じカラーマップで、数値と色の関係を変えたい場合はvmin, vmaxオプションを使います。たとえば

myimg = ax.imshow(my2Ddata, cmap = 'plasma', vmin=0, vmax=5000)
fig.colorbar(myimg)

のようにすれば、

plasmaカラーマップによるプロットの例

のような画像が得られるはずです。 ここで明るい色が指定した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')

すると

1次元化した実験データ

のようなグラフが得られます。もとのmy2Ddataとの対応関係を図示すると、

ravel()実行前後の対応関係

となります。つまり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()関数

• ravel()関数

• np.vstack()関数

を使います。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を作るためのものですが、イマイチわかりにくいので、何をやっているのか図示すると

meshgridの機能

となります。

次に、フィッティングで使う二次元ガウス関数

$$z(x,y)=A\exp(-(x-x_0)^2/\sigma_1^2)\exp(-(y-y_0)^2/\sigma_2^2)+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次元化された実験データとフィッティングの比較

のようになり、たしかにフィッティングは合っていそうですが、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')

とすれば、

実験データとフィッティングの比較

のように、かなりよく、フィッティングできていることがわかります。