Python / 1次元データの集合を2次元データに整形→カラーマップ作成

  • NumpyとMatplotlibを使って、1次元のデータの集合を2次元データに整形し、カラーマップ(カラープロットとかイメージプロットとも呼ばれる)を作ります。
  • Yを変えて取った沢山の1次元データ(X, Z)を2次元のデータ(X,Y,Z)に整形するところが本題です。

はじめに

実験をやっていると、測定量 $Z$ が2つのパラメータ $X$, $Y$ に依存するような状況がよくあります。こういう時、結果を $(X, Z)$ ないし $(Y, Z)$ のグラフで別々にまとめても良いですが、いっそ $(X,Y,Z)$ の2次元データにしてしまって、等高線プロットやカラーマップ一枚にまとめてしまうのも手です。

ただ、この「 $(X,Y,Z)$ の2次元データにする」というのが少し厄介で、普通実験では例え $Z$ が $Z=f(X, Y)$ という風に書けても、測定では $X$ と$Y$ のどちらかは固定しないといけないので、得られるのは $Z=f(X, Y_1), Z=f(X, Y_2), \cdots, Z=f(X, Y_N)$ のようなデータです。$Z=f(X,Y)$ の2次元(メッシュ状)データには自分で整形しなければいけません。

図1.png

ということで、この記事では上図のように測定の都合上分断された1次元データを2次元のメッシュ状データに整形して、カラーマップに起こしたいと思います。

サンプルデータ

ここからは下の dataset.py を実行してできるdatファイルを使って説明します。
(データの中身はべつに重要じゃないので、そういうデータだと思って、記録されているデータの形だけ見ておいてください。)

dataset.py
import numpy as np

#パラメータの定義
intensity = 50 #強度
HWHM = 10 #半値半幅
a = 10 #データのばらつきの大きさ

#データファイルの作成
for Y in np.arange (0, 10.1, 0.1):
    filename = 'Y=' + str(Y) + '.dat'
    X0 = (200*Y*Y + 2500)**0.5 - 50
    with open(filename, 'w') as file:
        file.writelines('X' + '\t' + 'Z' +'\n')
        for X in range (0, 101):
            Z = intensity * HWHM**2 / ((X-X0)**2 + HWHM**2) + 20 + a * np.random.rand()
            file.writelines(str(X) + '\t' + str(Z) + '\n')

名称未設定-1.png

Y=0.0~10.0まで、0.1刻みで101個のdatファイルができます。試しに Y=5.5.dat を開くと左図のようになっています。右図はZをXについてプロットしたグラフです。この山の位置はXだけでなくYにも依存しています。

[改訂版] 101個の1次元データを2次元データに整形

pandasを使うと圧倒的に簡単になることがわかったので、改訂版を書きます。

In[1]
import re, glob
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
In[2]
filelist = glob.glob('*.dat')
df = pd.DataFrame()
for file in filelist:
    match = re.match('Y=(.*).dat', file)
    df_sub = pd.read_table(file)
    df[float(match.group(1))] = df_sub.Z
df.columns.name = 'Y'
df.index = df_sub.X
df
Out[2]
Y   0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ... 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9
X                                                                                   
0   71.307065   71.624977   72.568247   70.073268   71.388264   71.429283   69.430455   66.104600   64.251044   61.960019   ... 20.892164   21.724984   20.259025   21.625291   22.658143   22.641024   20.799494   21.042593   20.667364   20.451245
1   66.347248   65.184597   66.907600   67.807422   67.879276   70.401552   72.100718   72.617697   72.195462   70.692071   ... 20.409888   22.230631   21.106551   22.801198   21.159110   20.973036   21.779757   20.625188   21.405971   21.577096
2   54.815612   54.960281   55.477640   57.619689   59.971637   60.601975   63.984228   65.729155   67.846441   69.637961   ... 21.854349   20.668861   20.172761   20.416828   21.374005   21.202518   21.688063   21.056256   22.637612   20.305400
3   46.311290   46.916455   47.512971   47.175870   48.731614   50.673641   52.572572   55.803255   59.562894   62.597950   ... 22.427942   20.156526   21.141887   22.187281   21.712688   22.921697   22.876228   22.972608   22.592168   21.185094
4   40.442910   38.820936   38.994950   41.859569   40.333883   42.995725   45.152994   47.650007   49.414120   53.309453   ... 21.873397   20.659303   21.022158   20.543980   23.023661   21.418374   22.771670   20.218522   22.349163   21.412955
5   35.598731   33.633262   35.782304   36.352184   35.778557   37.035520   38.947890   40.425551   41.307669   43.021355   ... 22.766781   20.876074   20.208458   21.890359   22.792392   22.499805   22.652404   22.497508   22.339281   21.668357
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

95  21.804124   21.797596   23.010401   20.258773   20.073975   21.918238   22.169896   21.988170   20.120070   20.975194   ... 28.402855   28.041022   32.390451   36.689843   48.158836   58.781772   71.028003   67.065172   52.942427   42.138778
96  21.250080   21.630843   21.553191   22.952056   21.329605   20.270100   21.658320   20.191202   21.166837   20.145893   ... 27.438835   27.227526   31.282917   32.428334   40.363609   50.900390   63.092738   70.580351   63.513627   48.526807
97  22.389748   21.693057   20.886997   21.460203   22.610140   20.102447   23.021290   22.793081   22.306881   20.704143   ... 25.180590   26.366878   28.042743   30.121939   36.192960   40.735346   53.298041   67.425563   70.242088   60.144091
98  21.265201   21.367930   21.225976   20.466155   21.115541   20.294466   20.556839   22.789051   20.945778   21.343996   ... 23.949042   24.200023   26.902070   28.732446   32.388426   35.483425   44.613251   56.697203   68.448203   70.253491
99  22.688459   22.243006   22.604197   22.114754   22.967067   22.538572   21.954847   22.286714   22.779653   20.139557   ... 23.386930   24.568997   25.573001   27.852061   30.095048   32.057468   37.229132   47.773504   61.152210   71.645167
100 21.503768   20.480336   21.507903   21.943483   21.158995   20.880028   22.613661   21.468507   22.059082   20.855645   ... 24.196745   25.333946   24.965214   27.570366   27.059141   31.368592   34.169847   40.928392   48.591212   62.390900
101 rows × 101 columns

データフレームを作って、ファイルから一つ一つZ列(2列目)を読んで追加していきます。
その際、カラム名にはファイル名から取得したYの値を設定しています。

In[3]
plt.pcolor(df.index, df.columns, df.T)
plt.colorbar()
plt.axis('tight')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

download-24.png

場合によっては以下の[改定前]から、In[2]でやっているfilelistのソートが必要になるかもしれません。

[改訂前]

前バージョンもそのまま残しておきます。
pandasを使わない縛りプレイでは役に立つかも・・・?

ライブラリの読み込み

In[1]
import numpy as np
import re, glob
from functools import cmp_to_key
import matplotlib.pyplot as plt
%matplotlib inline

処理するファイルのリストを作成

In[2]
filelist = glob.glob('*.dat')

#filelistのソート
filelist2 = []
for file in filelist:
    match = re.match('Y=(.*).dat', file)
    tuple = (match.group(1), file)
    filelist2.append(tuple)

def cmp(a, b):
    if a == b: return 0
    return -1 if a < b else 1

def cmptuple(a, b):
    return cmp(float(a[0]), float(b[0]))

filelist = [x[1] for x in sorted(filelist2, key=cmp_to_key(cmptuple))]
print(filelist)
In[2]の結果
['Y=0.0.dat', 'Y=0.1.dat', 'Y=0.2.dat', 'Y=0.3.dat', 'Y=0.4.dat', 'Y=0.5.dat', 'Y=0.6.dat', 'Y=0.7.dat', 'Y=0.8.dat', 'Y=0.9.dat', 'Y=1.0.dat', 'Y=1.1.dat', 'Y=1.2.dat', 'Y=1.3.dat', 'Y=1.4.dat', 'Y=1.5.dat', 'Y=1.6.dat', 'Y=1.7.dat', 'Y=1.8.dat', 'Y=1.9.dat', 'Y=2.0.dat', 'Y=2.1.dat', 'Y=2.2.dat', 'Y=2.3.dat', 'Y=2.4.dat', 'Y=2.5.dat', 'Y=2.6.dat', 'Y=2.7.dat', 'Y=2.8.dat', 'Y=2.9.dat', 'Y=3.0.dat', 'Y=3.1.dat', 'Y=3.2.dat', 'Y=3.3.dat', 'Y=3.4.dat', 'Y=3.5.dat', 'Y=3.6.dat', 'Y=3.7.dat', 'Y=3.8.dat', 'Y=3.9.dat', 'Y=4.0.dat', 'Y=4.1.dat', 'Y=4.2.dat', 'Y=4.3.dat', 'Y=4.4.dat', 'Y=4.5.dat', 'Y=4.6.dat', 'Y=4.7.dat', 'Y=4.8.dat', 'Y=4.9.dat', 'Y=5.0.dat', 'Y=5.1.dat', 'Y=5.2.dat', 'Y=5.3.dat', 'Y=5.4.dat', 'Y=5.5.dat', 'Y=5.6.dat', 'Y=5.7.dat', 'Y=5.8.dat', 'Y=5.9.dat', 'Y=6.0.dat', 'Y=6.1.dat', 'Y=6.2.dat', 'Y=6.3.dat', 'Y=6.4.dat', 'Y=6.5.dat', 'Y=6.6.dat', 'Y=6.7.dat', 'Y=6.8.dat', 'Y=6.9.dat', 'Y=7.0.dat', 'Y=7.1.dat', 'Y=7.2.dat', 'Y=7.3.dat', 'Y=7.4.dat', 'Y=7.5.dat', 'Y=7.6.dat', 'Y=7.7.dat', 'Y=7.8.dat', 'Y=7.9.dat', 'Y=8.0.dat', 'Y=8.1.dat', 'Y=8.2.dat', 'Y=8.3.dat', 'Y=8.4.dat', 'Y=8.5.dat', 'Y=8.6.dat', 'Y=8.7.dat', 'Y=8.8.dat', 'Y=8.9.dat', 'Y=9.0.dat', 'Y=9.1.dat', 'Y=9.2.dat', 'Y=9.3.dat', 'Y=9.4.dat', 'Y=9.5.dat', 'Y=9.6.dat', 'Y=9.7.dat', 'Y=9.8.dat', 'Y=9.9.dat', 'Y=10.0.dat']

イテレーターで処理するために、globモジュールを使って処理するdatファイルのリストを作成します。今回Yの値はファイル名に記録していますが、ファイルはこのYが小さい順に処理していきたいので、正規表現(と言っても、ワイルドカード.*を使うだけ)でYの値を一旦取り出して、リスト内のファイルをソートしています。ln[2]の結果を見ての通り、Yが小さい順に並びました。

Y軸のための一次元配列を作成

In[3]
Y = np.array([])
for file in filelist:
    match = re.match('Y=(.*).dat', file)
    Y = np.append(Y, float(match.group(1)))
print(Y)
In[3]の結果
[  0.    0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.    1.1
   1.2   1.3   1.4   1.5   1.6   1.7   1.8   1.9   2.    2.1   2.2   2.3
   2.4   2.5   2.6   2.7   2.8   2.9   3.    3.1   3.2   3.3   3.4   3.5
   3.6   3.7   3.8   3.9   4.    4.1   4.2   4.3   4.4   4.5   4.6   4.7
   4.8   4.9   5.    5.1   5.2   5.3   5.4   5.5   5.6   5.7   5.8   5.9
   6.    6.1   6.2   6.3   6.4   6.5   6.6   6.7   6.8   6.9   7.    7.1
   7.2   7.3   7.4   7.5   7.6   7.7   7.8   7.9   8.    8.1   8.2   8.3
   8.4   8.5   8.6   8.7   8.8   8.9   9.    9.1   9.2   9.3   9.4   9.5
   9.6   9.7   9.8   9.9  10. ]

リスト作成のところと一部重複しますが、後でY軸を設定するためにファイル名からYの数値だけ抽出して配列Yに納めます。

X軸のための一次元配列とカラーマップのための二次元配列を作成

In[4]
X = np.array([])
Z = np.empty((0,101), float)

#カラーマップのための配列を作成
for file in filelist:
    Z_sub = np.array([])
    with open(file, 'r') as read:
            for line in read.readlines()[1:]: #1行目をスキップ
                items = line[:-1].split('\t')
                Z_sub = np.append(Z_sub, float(items[1]))
    Z = np.vstack((Z, Z_sub))

#X軸のための配列を作成
with open(file, 'r') as read:
            for line in read.readlines()[1:]: #1行目をスキップ
                items = line[:-1].split('\t')
                X = np.append(X, float(items[0]))

print(X)
print(Z)
In[4]の結果
[   0.    1.    2.    3.    4.    5.    6.    7.    8.    9.   10.   11.
   12.   13.   14.   15.   16.   17.   18.   19.   20.   21.   22.   23.
   24.   25.   26.   27.   28.   29.   30.   31.   32.   33.   34.   35.
   36.   37.   38.   39.   40.   41.   42.   43.   44.   45.   46.   47.
   48.   49.   50.   51.   52.   53.   54.   55.   56.   57.   58.   59.
   60.   61.   62.   63.   64.   65.   66.   67.   68.   69.   70.   71.
   72.   73.   74.   75.   76.   77.   78.   79.   80.   81.   82.   83.
   84.   85.   86.   87.   88.   89.   90.   91.   92.   93.   94.   95.
   96.   97.   98.   99.  100.]
[[ 77.80746092  70.89995942  75.85870509 ...,  23.75101033  26.63014195
   27.41822976]
 [ 70.22928442  77.20944503  68.4186785  ...,  23.11059562  25.32787551
   23.83316096]
 [ 78.82482384  75.82282473  69.04286475 ...,  28.4786889   28.88143063
   29.02818251]
 ..., 
 [ 28.05311936  25.86981344  27.63047642 ...,  79.70742496  74.3232517
   72.44646855]
 [ 26.02612835  22.579358    26.90219008 ...,  72.76800112  76.14112188
   74.62499125]
 [ 29.6708667   22.40732954  29.63804179 ...,  71.55372667  76.75416415
   75.76960793]]

ln[2]で作成したfilelistでイテレーションを行い、datファイルの2列目(Zの値)を逐一読み出して2次元配列 Z に納めます。少し詳しく説明します。まず、この配列Zは扱うデータによって初期化が必要で、今回はdatファイルのデータ行数がX=0~100までの101行なので

Z = np.empty((0,101), float)

というように(0,101) で初期化しています(一番初めの図で示したように、datファイルで行だったものが列に転置されます)。filelistでのイテレーションが始まって、

for line in read.readlines()[1:]: #1行目をスキップ
    items = line[:-1].split('\t')
    Z_sub = np.append(Z_sub, float(items[1]))

のところですが、ファイルを一つ開いて一行ずつ一列目の要素をZ_subに格納していきます。#1行目をスキップと書いたところは、今回datファイルの一行目に「X Z」という文字列をインデックスとして書き込んでおり、それを読み込ませないためです。最終行まで行ったら、このZ_subは1×101の配列になるわけですが、今度はこれをZに縦に積んでいきます。それをfilelistの最後のファイルまで繰り返します。

Zを作り終えたところで最後に開いたファイルをもう一度開き、X軸作成のために一列目の値(X)を配列Xに格納します。

カラーマップの作成

In[5]
plt.pcolor(X, Y, Z) #カラーマップを作成
plt.colorbar() #カラーバーを表示
plt.axis('tight')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

2次元のZができたので、plt.pcolor(X, Y, Z)でカラーマップを作ります。改めて書きますが、X はX軸のための一次元配列、Y はY軸のための一次元配列、Z はカラーマップのために作った二次元配列です。plt.pcolor(Z)だけでもカラーマップはできるんですが、スケールがわからなくなるので、軸のための配列X, Yを用意しています。plt.axis('tight')は、カラーマップの余白を消すための記述です。

ということで、下図がIn[5]の結果です。Jupyter上ではインラインで表示されます。

figure_1.png

めでたくカラーマップができました。