Python を用いて実験データを読み込み、可視化する

# まえがき

この記事は、実験等で得られる離散的な数値データをndarray として読み込み、matplotlib でグラフや強度マップを描く流れをざっくりまとめたものです。実験データはテキストファイルとします。（バイナリファイルでもフォーマットさえ分かっていれば読み込めますが、ちょっと面倒です。）

この記事で示すコードでは、事前に以下のようにライブラリがimport してあるものとします。

```Python
import numpy as np
import matplotlib.pyplot as plt
```

また、Jupyter Notebook を使うことをおすすめします。その場合は、最初のセルで以下のようなマジックコマンドを実行すると良いかと思います。

```
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
```

# データの読み込み

何らかの装置でデータを取得した場合、出力されるテキストファイルは以下のような形式になっていることが多いでしょう。

```
# comment
# axis1 axis2 value1 value2
 0.1   10.0   401   102
 0.2   11.0   392   232
 0.3   12.0   397   574
 0.4   13.0   404   896
...
```

すなわち、**何行かのコメント（ヘッダー）の後に、数値データが何列かで記載されている**というものです。このような形式のデータは、np.loadtxt もしくはnp.genfromtxt でndarray として読み込むことができます。

## np.loadtxt

まず、np.loadtxt によるデータの読み込み方を確認します。以下のような内容のdata01.txt というファイルがPython のスクリプトファイルもしくはJupyter Notebook ファイルと同じディレクトリにあるとします。

```:data01.txt
! comment line1
! comment line2
0  0
1  1
2  4
3  9
4  16
5  25
```

このファイルをnp.loadtxt で読み込むときは、以下のようにします。comments にコメント行の先頭文字（ここでは!）を指定しています。（デフォルトでは'#' に設定されているようです。）

```Python
data01 = np.loadtxt("./data01.txt", comments='!')
```

読み込んだデータは以下のようなndarray になっています。dtype は特に指定しなければ、np.float64 になります。

```Python
array([[  0.,   0.],
       [  1.,   1.],
       [  2.,   4.],
       [  3.,   9.],
       [  4.,  16.],
       [  5.,  25.]])
```

元のテキストデータの1列目は```data01[:, 0]```、2列目は```data01[:, 1]``` でアクセスできますが、ちょっと気持ち悪いですね。データ列ごとに一つの1 次元ndarray になってくれると直観的に扱えそうです。これは、unpack というパラメータをTrue にすることで実現できます。

```Python
data01_axis1, data01_value1 = np.loadtxt("./data01.txt", comments='!', unpack=True)
```

このようにすると、data01_axis1 に1列目のデータが、data01_value1 に2列目のデータが、それぞれ1次元ndarray として入ります。

実際の実験データのファイルは、上記のdata01.txt のように整ったものとは限りません。以下では、形式がいまいち整っていないファイルを読み込む方法を確認します。

たとえば、以下ののようにコメント行に特定の先頭文字がない場合があります。これをdata02.txt としましょう。

```:data02.txt
comment line1
comment line2
0  0
1  1
2  4
3  9
4  16
5  25
```

このようなファイルは、comments では対応できません。その代わり、skiprows を用いることで、指定行数だけ先頭をスキップして読み込むことが可能です。

```Python
data02_axis1, data02_value1 = np.loadtxt("./data02.txt", skiprows=2, unpack=True)
```

また、区切り文字がある場合にも対応できます。たとえば、以下のようなカンマ区切りのdata03.txt があるとしましょう。

```:data03.txt
0,  0
1,  1
2,  4
3,  9
4,  16
5,  25
```

このようなファイルを読み込む際には、区切り文字をdelimiter に指定する必要があります。

```Python
data03_axis1, data03_value1 = np.loadtxt("./data03.txt", delimiter=',', unpack=True)
```

delimiter は文字列で指定するようになっており、デフォルトは空白文字です。また、data03.txt のように、区切り文字（この場合はカンマ）と空白文字が混ざっている場合は、文字列間の空白文字まで含めて区切り文字と認識するようです。

より複雑な場合として、欠損値があったり文字列を含む列がある場合なども考えられますが、そのような場合にはnp.loadtxt では対応できません。代わりに、np.genfromtxt を利用しましょう。

## np.genfromtxt

np.genfromtxt は、np.loadtxt よりも複雑な条件付けができるデータ読み込み用のメソッドです。主に欠損値に対応することを目的としているようです。パラメータはnp.loadtxt と共通のものが多いので、新しく覚えなければいけないことは、そう多くはありません。

まず、欠損値の対応の仕方を見てみましょう。以下のdata04.txt を用意します。

```:data04.txt
0,  0
1,  1
2,  4
3,  
4,  16
5,  25
```

4行目にデータの欠損があります。これをgenfromtxt で読み込みます。

```Python
data04 = np.genfromtxt("./data04.txt", delimiter=',')
```

data04 は以下のようになっており、欠損値はnan(not a number) で埋まっていることがわかります。

```
array([[  0.,   0.],
       [  1.,   1.],
       [  2.,   4.],
       [  3.,  nan],
       [  4.,  16.],
       [  5.,  25.]])
```

もちろん、unpack した場合も、対応する位置のみnan で埋まります。ちなみに、文字列が混ざっている場合にも欠損値として扱われるので、その部分はnan で埋まります。

読み込んだ後の処理によっては、欠損値部分がnan だと困ることがありそうです。その場合は、filling_values を指定することで、欠損値をnan 以外の値で埋めることができます。たとえば、

```Python
data04 = np.genfromtxt("./data04.txt", delimiter=',', filling_values=0)
```

とすると、data04 は

```
array([[  0.,   0.],
       [  1.,   1.],
       [  2.,   4.],
       [  3.,   0.],
       [  4.,  16.],
       [  5.,  25.]])
```

と欠損値が0 で埋まって読み込まれます。

delimiter に整数値が指定できるようになっていることも、np.loadtxt からの重要な変更点です。これによって、列の区切り方を文字数で指定することができます。以下のような各列が3文字のフィールドになっているdata05.txt があるとしましょう。6行目では文字がつながっているので、空白文字区切りでは分割できず、エラーとなります。

```:data05.txt
  0  0
  1  1
  2  8
  3 27
  4 64
  5125
```

この場合は、

```Python
data05 = np.genfromtxt("./data05.txt", delimiter=3)
```

とすれば読み込めます。列ごとにフィールドの長さ（文字数）が違うときは、delimiter にint の配列を指定すれば対応可能です。

## より複雑なフォーマットのデータの読み込み

np.genfromtxt でも対応できないようなフォーマットのデータは、既存のメソッドでndarray として読み込むことは困難です。

そのような場合には、ファイルを読み込んで特定のルールに従って数値に変換する処理を自分で作る必要があります。基本的な流れとしては、ルールに従って文字列から数値配列を作り、それをnp.array でndarray へ変換という感じになります。

# データの可視化

ここからはmatplotlib によるデータの可視化の方法を確認します。

matplotlib には、pylab スタイルとオブジェクト指向スタイルという二種類のコーディング形式があります。このあたりの話は色々わかりやすい記事[^1]があるので省略しますが、ここではオブジェクト指向スタイルに統一します。

[^1]: たとえば、「matplotlibをオブジェクト指向スタイルで使う その1」 (http://minus9d.hatenablog.com/entry/2015/01/26/225731) など。

matplotlib ではさまざまな方法・手順でグラフを描画できますが、この記事では、

1. plt.figure() でfigure を作る。
1. 作ったfigure に対して、一つもしくは複数のグラフ描画領域(axes) をplt.add_subplot() で追加する。
1. 描画領域にデータをプロットする。

という流れを採用します。

## 1次元データの可視化（ラインプロファイル）

1次元データのグラフ化には、plot() メソッドを用います。data01 のvalue1 をaxis1 に対してプロットしてみましょう。

```Python
data01_axis1, data01_value1 = np.loadtxt("./data01.txt", comments='!', unpack=True)

fig = plt.figure(figsize=(4, 6))
ax = fig.add_subplot(111)
ax.plot(data01_axis1, data01_value1, "o-", color="k", label="value1 of data01")
ax.set_xlim(-0.5, 5.5)
ax.set_ylim(-0.5, 25.5)
ax.set_xlabel("axis1")
ax.set_ylabel("value1")
ax.legend(loc="upper left")
plt.show()
```

<img src="https://qiita-image-store.s3.amazonaws.com/0/169980/0453695a-88f6-fb8b-eebc-a2fcd685ff99.png" width=400>

グラフが出ました。詳しく流れを追ってみましょう。まず、

```Python
fig = plt.figure(figsize=(4, 6))
```

でfigure クラスのインスタンスであるfig を作っています。figsize に指定しているタプルの中身は、グラフのサイズの縦と横で、単位はインチです。

次に、

```Python
ax = fig.add_subplot(111)
```

でfig にグラフ描画領域であるaxes を追加しています。このとき、引数に111 を入れていますが、ここには任意の三桁の数字を指定することができます。この数字をmni とすると、図の領域を縦m × 横n の区間に分けて、左上から数えてi 番目にグラフ描画領域を追加する、という意味になります。今回は111 なので、図は分割せず、一つのグラフ描画領域だけを設定することになります。

そして、

```Python
ax.plot(data01_axis1, data01_value1, "o-", color="k", label="value1 of data01")
```

でデータを描画領域であるax にプロットします。上記の場合、data01_axis1 がx 軸に、data01_value1 がy 軸に対応します。"o-" はマーカーと線の種類(丸のマーカーと実線)を表しており、color は色、label は凡例に表示される文字列です。

その後、

```Python
ax.set_xlim(-0.5, 5.5)
ax.set_ylim(-0.5, 25.5)
ax.set_xlabel("axis1")
ax.set_ylabel("value1")
ax.legend(loc="upper left")
```

でx,y の表示範囲およびラベル、凡例の設定を行います。

最後に

```Python
plt.show()
```

を実行するとグラフが表示されます。(Jupyter Notebook を利用している場合はこの一行は必須ありませんが、もちろん実行しても問題ありません。)

線の色や種類の変更、凡例の位置やフォントの指定、任意の位置に任意の文字列を追加等々、だいたい何でも出来ますので必要に応じてドキュメントを参照してください。

### X 線回折プロファイルの例

例として、多結晶Si 粉末のX線回折(X-Ray Diffraction; XRD) の測定結果[^2]をグラフにしたものを以下に示します。
データは"Poly_Si_XRD.dat" というファイル名で、以下のようなカンマ区切りの2列データのテキストファイルになっています。
一列目が横軸$2\theta (^\circ{})$ で、二列目が縦軸の強度(counts) です。

[^2]: これ以外ちょうど良いデータが手持ちになかったのですが、そもそもXRDってなんやねんという方はWikipedia (http://ja.wikipedia.org/wiki/X線回折) をどうぞ。

```:Poly_Si_XRD.dat
 10.00 ,     11
 10.02 ,     19
 10.04 ,     22
 10.06 ,     25
 10.08 ,     18
...
```

numpy を利用する利点の一つは、豊富な数値計算メソッドを使って読み込んだデータの変換が容易にできることです。ここでは、以下の式で横軸を変換した結果を3 つ並べたグラフを作成してみました。[^3]

[^3]: ここでは、$\lambda = 1.5406$ Å としました。

$$
q = \frac{4\pi \sin\theta}{\lambda}
$$

$$
d = \frac{2\pi}{q}
$$

手順は簡単で、読み込んだ横軸のtwotheta から、変換式に従って新たなndarray であるq, d を計算するだけです。

また、matplotlib では、**軸ラベル等のテキストを記述する際にLaTeX 表記が使えます** 。やり方は、文字列の中のLaTeX 記法で入力したいところを$ で囲むだけです。これによって、ギリシャ文字や特殊記号を簡単に表示させることが出来ます。

一つ注意が必要なのが、今回のようにfig.add_subplot() によって一つの図に複数のグラフを表示させるときは、グラフ同士の一部が重なることを避けるためにfig.tight_layout() を使う必要があるということです。これを怠ると、軸のラベルの一部が重なったりします。

```Python
twotheta, intensity = np.genfromtxt("./data/Poly_Si_XRD.dat", delimiter=',', unpack=True)
lmd = 1.5406
q = 4 * np.pi * np.sin(np.radians(twotheta/2)) / lmd
d = 2 * np.pi / q

fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(311)
ax1.plot(twotheta, intensity, '-', color='r')
ax1.set_xlabel("$2 \\theta$ / $^{\circ}$"); ax1.set_ylabel("Intensity / counts")
ax1.set_ylim(-500, 7000); ax1.set_xlim(10, 90)
ax2 = fig.add_subplot(312)
ax2.plot(q, intensity, '-', color='g')
ax2.set_xlabel("$q$ $(=4\pi\sin{\\theta}/\lambda$) / $\mathrm{\AA}^{-1}$"); ax2.set_ylabel("Intensity / counts")
ax2.set_ylim(-500, 7000); ax2.set_xlim(0.8, 5.8)
ax3 = fig.add_subplot(313)
ax3.plot(d, intensity, '-', color='b')
ax3.set_xlabel("$d$ / $\mathrm{\AA}$"); ax3.set_ylabel("Intensity / counts")
ax3.set_ylim(-500, 7000); ax3.set_xlim(1.05, 3.5)
fig.tight_layout(pad=0)
plt.show()
```
<img src="https://qiita-image-store.s3.amazonaws.com/0/169980/d37ff0b1-4ba0-dec3-3e1a-f2f3b19e5344.png" width=700>

特に問題なくグラフ化できていることが分かりますね。(XRD のデータで横軸d のグラフはあまり作らないと思いますが、それはさておき。)

## 2次元データの可視化（強度マップ）

実験をしていると、2つの独立なパラメータを変化させながら、何らかの量を計測したいことがあります。たとえば、試料表面上の位置X, Y を変えながら、各点で何らかの量Z を測定する場合です。このようなマッピング測定で得られるデータはX, Y, Z の3 次元になるので、直接可視化しようとすると等高線図のような2 次元強度マップにせざるを得ません。 

matplotlib ではいくつか2 次元強度マップを描くメソッドがありますが、ここではそのうちの一つであるpcolormesh について説明します。

このメソッドは、引数に2 次元ndarray を取ります。ここまで考えていたファイルの読み込み方だと、データは1 次元ndarray なので、2 次元ndarray への変換が必要になります。また、それだけでなく、**X, Y のサイズはZ よりも縦横方向に1 ずつだけ大きい必要があります。**

これに関しては、以下のサイトに分かりやすい図があります。
**[[python][matplotlib]matplotlibで等高線やら何やらのプロット](http://d.hatena.ne.jp/y_n_c/20091122/1258904025)**

pcolormesh はある点をプロットする際にその要素とその次の要素で囲まれた領域を塗り潰すため、X, Y はZ よりも縦横方向に1 ずつだけ大きい必要があるわけです。（植木算ですね。）**ただ、データ点数が多い場合は殆ど気にならないので、無視しても構いません。** 端点は欠けますが、X, Y, Z が同じサイズでもエラーは出ず、描画してくれます。

このようにpcolormesh の仕様は複雑なので、何段階かに分けて考えた方が分かりやすくなります。ここでは、以下の手順で進めます。

1. データを読み込む。
1. 読み込んだ1 次元ndarray を2 次元ndarray に変換する。
1. X, Y に相当するndarray のサイズを縦横方向に1 ずつ増やす。
1. 強度マップを描画する。

それでは、以下の内容のdata06.txt を読み込んで、2次元強度マップにしてみましょう。

```data06.txt
# X Y Z
-2 -2 -0.832
-1 -2 0.124
0 -2 0.584
1 -2 0.124
2 -2 -0.832
-2 -1 0.124
-1 -1 1.081
0 -1 1.540
1 -1 1.081
2 -1 0.124
-2 0 0.584
-1 0 1.540
0 0 2.000
1 0 1.540
2 0 0.584
-2 1 0.124
-1 1 1.081
0 1 1.540
1 1 1.081
2 1 0.124
-2 2 -0.832
-1 2 0.124
0 2 0.584
1 2 0.124
2 2 -0.832
```

まず、データを読み込みます。3 列をx, y, z にunpack します。

```Python
x, y, z = np.loadtxt("./data06.txt", unpack=True)
```

つぎに、x, y, z を2 次元配列X, Y, Z に変換します。このとき、変換後の2 次元配列の縦横のサイズが必要になります。今回の場合は、5×5 のサイズになります。

```
X = x.reshape(5, 5)
Y = y.reshape(5, 5)
Z = z.reshape(5, 5)
```

X, Y, Z をそれぞれ確認してみましょう。

```Python
# X
array([[-2., -1.,  0.,  1.,  2.],
       [-2., -1.,  0.,  1.,  2.],
       [-2., -1.,  0.,  1.,  2.],
       [-2., -1.,  0.,  1.,  2.],
       [-2., -1.,  0.,  1.,  2.]])
```

```Python
# Y
array([[-2., -2., -2., -2., -2.],
       [-1., -1., -1., -1., -1.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  1.,  1.],
       [ 2.,  2.,  2.,  2.,  2.]])
```


```Python
# Z
array([[-0.832,  0.124,  0.584,  0.124, -0.832],
       [ 0.124,  1.081,  1.54 ,  1.081,  0.124],
       [ 0.584,  1.54 ,  2.   ,  1.54 ,  0.584],
       [ 0.124,  1.081,  1.54 ,  1.081,  0.124],
       [-0.832,  0.124,  0.584,  0.124, -0.832]])
```

前述のように、pcolormesh でZ を正しく描画するには、X, Y はZ よりも縦横方向に1 ずつ大きいサイズを持っている必要がありました。これは、実験データを読み込んで描画しようと思ったときには大きな障害になります。なぜなら、読み込む時点ではX, Y, Z は同じサイズであることが殆どだからです。ただし、今回のような等間隔に並んだ簡単なデータであれば、ndarray の操作によって対応できます。

**2017/11/23追記: コメントで情報を頂きましたが、以下の処理はpandas を使うと簡単にできるようです。**

まず、X からサイズを増やしましょう。各点の間隔は1.0 で、追加すべきなのは上記の表示でみて下と右(index が大きい方)なので、

```Python
X = np.vstack((X, X[0]))
tmpX =  np.ones(len(X[:, 0])) * (X.max()+1.)
X = np.hstack((X, tmpY[:, np.newaxis]))
```

とします。すると、X は

```Python
array([[-2., -1.,  0.,  1.,  2.,  3.],
       [-2., -1.,  0.,  1.,  2.,  3.],
       [-2., -1.,  0.,  1.,  2.,  3.],
       [-2., -1.,  0.,  1.,  2.,  3.],
       [-2., -1.,  0.,  1.,  2.,  3.],
       [-2., -1.,  0.,  1.,  2.,  3.]])
```

となります。Y も同様に

```Python
Y = np.hstack((Y, Y[:, 0][:, np.newaxis]))
tmpY =  np.ones(len(Y[0])) * (Y.max()+1.)
Y = np.vstack((Y, tmpY))
```

とすると

```Python
array([[-2., -2., -2., -2., -2., -2.],
       [-1., -1., -1., -1., -1., -1.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  1.,  1.,  1.],
       [ 2.,  2.,  2.,  2.,  2.,  2.],
       [ 3.,  3.,  3.,  3.,  3.,  3.]])
```

となります。今回の例では簡単ですが、各点間の間隔が毎回違うなどの複雑な場合は、データ自体を手動で編集したほうが早いかもしれません。また、先にも書いたように、点数が多いときは無視しても大きな問題はありません。

さて、最後に強度マップを描画してみましょう。cmap にカラーマップ（色の付け方）の種類を指定可能です。好みの問題かと思いますが、ここではjet にしておきます。

```Python
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.pcolormesh(X, Y, Z, cmap='jet')
plt.show()
```
　<img src="https://qiita-image-store.s3.amazonaws.com/0/169980/f54f10b2-ceea-9797-a2ea-38c9d21b761c.png" width=350>

pcolormesh の仕様によって、X, Y がずれたように見えているのが気持ち悪いですね。X, Y が示す位置を中心として塗りつぶしたい場合は、X, Y それぞれからステップの半分を引いてやればOKです。それも含めて、コードと結果を以下にまとめて示します。

```Python
x, y, z = np.loadtxt("./data06.txt", unpack=True)

X = x.reshape(5, 5)
Y = y.reshape(5, 5)
Z = z.reshape(5, 5)

X = np.hstack((X, X[:, 0][:, np.newaxis]))
tmpX =  np.ones(len(X[0])) * (X.max()+1.)
X = np.vstack((X, tmpX))

Y = np.vstack((Y, Y[0]))
tmpY =  np.ones(len(Y[:, 0])) * (Y.max()+1.)
Y = np.hstack((Y, tmpY[:, np.newaxis]))

fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.pcolormesh(X-0.5, Y-0.5, Z, cmap='jet')
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-2.5, 2.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()
```

　<img src="https://qiita-image-store.s3.amazonaws.com/0/169980/28fbaa6f-14ef-b1cc-3cfb-d4ce668c3fba.png" width=350>

強度マップを示す際は強度を表すカラーバーを追加したくなります。しかし、matplotlib で強度マップと大きさが揃ったカラーバーを追加するのは意外と面倒です。これには、以下のようにmpl_toolkits.axes_grid1.make_axes_locatable を利用して、カラーバー用のaxes を追加する必要があります。

```Python
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
intensity_map = ax.pcolormesh(X-0.5, Y-0.5, Z, cmap='jet')

divider = make_axes_locatable(ax)
colorbar_ax = divider.append_axes("right", "5%", pad="3%")
fig.add_axes(colorbar_ax)
fig.colorbar(intensity_map, cax=colorbar_ax)

ax.set_xlabel("X"); ax.set_ylabel("Y"); ax.set_aspect('equal')

plt.show()
```

　<img src="https://qiita-image-store.s3.amazonaws.com/0/169980/c09a730c-516b-78fa-08ce-c370acf9d0ee.png" width=400>

二次元強度マップの縦サイズとぴったり合ったカラーバーが表示できました。

### 二次元X線回折像の例

ここでは、二次元X 線回折像[^4]を可視化する例を紹介したいと思います...が、ちょうど良いデータが無かったので適当にそれっぽいデータ（粉末を想定）を生成するところから始めます。

[^4]: またまたXRD の例で申し訳ないですが、他の例を思いつかなかったんです...

以下のようにしてデータを作ります。

```Python
# make X, Y, Z
x = np.arange(1, 501, 1, dtype=np.float64)
y = np.arange(1, 501, 1, dtype=np.float64)
X, Y = np.meshgrid(x, y)
Z = np.zeros(X.shape, dtype=np.float64)

# make polar coordinate axes
center = (50, 50)
R = np.sqrt( (X-center[0])**2 + (Y-center[1])**2 )
TH = np.arctan2( (Y-center[1]), (X-center[0]) )

# make diffraction rings
mask1 = np.logical_and( R > 60, R < 63 )
mask2 = np.logical_and( R > 120, R < 123 )
mask3 = np.logical_and( R > 300, R < 303 )
mask4 = np.logical_and( R > 420, R < 423 )
Z[mask1] = 500.0; Z[mask2] = 700.0; Z[mask3] = 1000.0; Z[mask4] = 500.0
from scipy.ndimage.filters import gaussian_filter
Z = gaussian_filter(Z, 1.0)

# add background
Z += 500.0*np.exp(-0.001*R)

# output text file
output = np.vstack( (X.flatten(), Y.flatten(), Z.flatten()) ).T
np.savetxt("data07.txt", output, fmt=["%.0f", "%.0f", "%.3f"])
```

これで、data07.txt が生成されます。縦横軸それぞれ500点ずつの、計250000点のデータが3列で入ったデータになります。

データを読み込んで描画してみましょう。データ点数が多いので、端点の処理は無視できるとして、スキップします。

```Python
x, y, z = np.loadtxt("data07.txt", unpack=True)
X = x.reshape(500, 500)
Y = y.reshape(500, 500)
Z = z.reshape(500, 500)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
intensity_map = ax.pcolormesh(X, Y, Z, cmap='viridis')

divider = make_axes_locatable(ax)
colorbar_ax = divider.append_axes("right", "5%", pad="3%")
fig.add_axes(colorbar_ax)
fig.colorbar(intensity_map, cax=colorbar_ax)

ax.set_xlabel("X"); ax.set_ylabel("Y")
ax.set_aspect('equal')

plt.show()
```
<img width="422" alt="download.png" src="https://qiita-image-store.s3.amazonaws.com/0/169980/b9322fa9-e515-824b-3b0f-b13c4abb9100.png">

それっぽい強度マップが見えていますね。

さて、なぜこんな例を出したかというと、**1 次元プロファイルの例でみたような座標変換は、2 次元強度マップでも有効**であることを示したかったからです。ここでは、極座標系に変換したマップを表示してみます。

中心の座標を$(cx, cy)$ とすると、極座標系の座標$R, \Theta$ は単に

$$
R = \sqrt{ (x - cx)^2 + (y - cy)^2 }
$$

$$
\Theta = \arctan{ \frac{y - cy}{x - cx} }
$$

で計算できます。直観的には面倒に思えますが、実は2 次元ndarray であるX, Y に対して、numpy メソッドを使うことによって、計算式ほぼそのままでコードが書けます。このようなことが可能なのは、縦・横軸の値が計算によって元のX, Y から変化しても、Z に対するインデックスの対応は保たれるからです。

```Python
center = (50, 50)
R = np.sqrt( (X-center[0])**2 + (Y-center[1])**2 )
TH = np.degrees(np.arctan2( (Y-center[1]), (X-center[0]) ))

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
intensity_map = ax.pcolormesh(R, TH, Z, cmap='viridis')

divider = make_axes_locatable(ax)
colorbar_ax = divider.append_axes("right", "5%", pad="3%")
fig.add_axes(colorbar_ax)
fig.colorbar(intensity_map, cax=colorbar_ax)

ax.set_xlabel("$R$ / pixel"); ax.set_ylabel("$\Theta$ / $^{\circ}$")

plt.show()
```

<img width="431" alt="download.png" src="https://qiita-image-store.s3.amazonaws.com/0/169980/0b51d524-9694-c81c-04bb-cf914964b507.png">

極座標系で表した結果が描画できました。想定される展開図になっていますね。リングが途中で途切れているので$\Theta$ 方向には値が無い領域がありますが、そのような部分には背景色が表示されます。

# まとめ

Python による実験データの可視化方法をざっくりまとめました。

numpy/scipy 等の科学計算ライブラリを実験データに適用すると、大体何でも出来ます。さらにJupyter Notebook を一緒に使えば、これが無料で良いのかというぐらい便利なデータ解析環境になるので、高価なグラフソフトが買いたくない人/わざわざ特殊な言語を覚えたくない人等にお勧めです。