「HyperSpy」は、Pythonで記述された多次元データ解析のためのライブラリです。次のような特徴があります。
HyperSpyは非常に強力なライブラリですが、日本語の情報は、書籍もWebサイトもまったく見つかりません。英語の公式ドキュメントを見るのが確実という状況です。
https://hyperspy.org/hyperspy-doc/current/index.html
そこで本講義では、このHyperSpyについて入門レベルのご紹介をします。けっこう貴重な講義かも!?
はじめに、本講義で使用するファイルを皆さんの環境にダウンロードするため、次のコードを実行してください。
次のソースコードを、実行してください。Google Colaboratoryでソースコードを実行するには、ソースコードの任意の場所をクリックした後、次のどちらかの操作を行います。
!curl -L -o Fei_HAADF-DE_location.dm3 https://github.com/tendo-sms/python_intermediate_2023/raw/main/file_2/Fei_HAADF-DE_location.dm3
!curl -L -o core_shell.hdf5 https://github.com/tendo-sms/python_intermediate_2023/raw/main/file_2/core_shell.hdf5
HyperSpyの公式ページでは何種類かのインストール方法が紹介されていますが今回は、
という理由で、pipコマンドでのインストール方法を使用します。
!pip install hyperspy
pipコマンドの終了を待つ間に・・
HyperSpyでは、マウス操作で操作できるインタラクティブなグラフを描画できます。
Jupyter NotebookやJupyterLabを使用する場合、デフォルトではインタラクティブなグラフではなく、単なる画像に変換されて表示されます。
次のコマンドを実行することで、インタラクティブなグラフを描画できるようになります。
%matplotlib nbagg
しかし・・・実は、Google Colaboratoryでは、インタラクティブなグラフを描画する方法がありません。単なる画像の表示のみ対応しています。
本講義ではGoogle Colaboratoryを利用するため、皆さんのお手元でインタラクティブな操作を実感することはできません。そのため、講師のPCでJupyter Notebookを使ったインタラクティブな操作をご覧いただきます。
HyperSpyの機能によっては、GUI操作が求められるものがもあります。
この講義でもGUIの表示が必要なサンプルプログラムがあるので、次のパッケージをpipコマンドでインストールします。
!pip install hyperspy_gui_traitsui==1.5.3
事前準備はここまでです。
いよいよ、HyperSpyを使ってプログラムを記述していきます。
import hyperspy.api as hs
HyperSpyにはいくつかのサンプルデータが元から用意されています。とりあえず、次の2つのサンプルデータをプロットしてみましょう。
グラフのプロットは、次の構文です。
オブジェクト.plot()
次のコードを実行してください。
import hyperspy.api as hs
# SEMによるエネルギー分散型X線分析(EDS)データ
hs.datasets.example_signals.EDS_SEM_Spectrum().plot()
# 鉄針のホログラムデータ
hs.datasets.example_signals.object_hologram().plot()
メタデータの確認方法などは、このあと別の例でご説明します。
HyperSpyは、大別すると以下のフォーマットに対応しています。
バイナリファイル (Gatan’s dm3/dm4、EDAX spc and spd、MRCなど)
テキストファイル (ASCII、CSV、TSV、等間隔行列など)
HDF5形式ファイル
【ご参考】全ての対応フォーマットは公式ドキュメントの以下に記載があります。
https://hyperspy.org/hyperspy-doc/current/user_guide/io.html#supported-formats
Hyperspyは入力ファイルの拡張子やその内容から、ファイルの種類を推定し、適切なIOプラグインを選択してくれます。
今回の講義では、Gatanのdm3ファイルと、HDF5ファイルを題材としてご紹介します。
【ご参考】
hs.print_known_signal_types()
最初にダウンロードしたファイルの中に、STEMで計測された16 x 4の2次元データ「Fei_HAADF-DE_location.dm3」があります。
これを読み込んでプロットしてみましょう。
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
s.plot()
元の入力ファイル(HyperSpyの形式に変換される前のもの)に含まれていたメタデータは、次の変数に格納されています。
オブジェクト.original_metadata
先ほど読み込んだファイルに元から含まれていたメタデータを確認してみます。
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
s.original_metadata
【ワンポイント】
データを画面表示するにはprint関数を使うのが基本ですが、ここではGoogle Colaboratoryの独自機能を活用して、print関数を使わずにデータを画面表示しました。この独自機能は、単に変数名を記述するだけで、データを見やすく画面表示することができます。他にもJupyter Notebook等、Pythonが「対話モード」で動作している環境はこのような機能を持ち、処理方法を試行錯誤する際には非常に便利です。
一方で、処理方法確定時に作ることが多い、拡張子が.pyとなっている「普通のPythonプログラム」は「スクリプトモード」で動作し、変数名だけでは画面表示をしません。
Python初級者セミナーでは、「対話モード」特有機能を多用して変なクセを付けてしまわないように、データの画面表示は一貫してprint関数を使ってご説明しました。中級者編にご参加されている皆様は心配ないと思いますので、便利さ重視でソースコードを記述しています。
HyperSpyが元の入力ファイルから判断して自動的に付与したメタデータは、次の変数に格納されています。
オブジェクト.metadata
先ほど読み込んだファイルに対して、HyperSpyが自動的に付与したメタデータを確認してみます。
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
s.metadata
これらのメタデータを更新したり、新たに追加したりすることもできます。
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
# メタデータの更新
s.metadata.General.title = "Update title"
# メタデータの追加
s.metadata.New_value = "New value"
s.metadata
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
s.data
プロットされる実データは、NumPy配列です。
そのため、データに対してNumPyの様々な機能を使うことができます。
import numpy as np
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
# 配列形状の取得
print(s.data.shape)
# 平均値の取得
print(np.mean(s.data))
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
# 色々と試してみましょう
作成・編集したデータを保存するには、「saveメソッド」を使用します。
引数にファイル名を指定し、ファイル名の拡張子を参照しファイルフォーマットを決定します。
保存できるフォーマットは、下記のWriteがYesのものです。
https://hyperspy.org/hyperspy-doc/current/user_guide/io.html#supported-formats
拡張子を指定しない場合は、デフォルトではhspy(HyperSpy標準のHDF5形式)で保存します。
import hyperspy.api as hs
s = hs.load("Fei_HAADF-DE_location.dm3")
# メタデータの更新
s.metadata.General.title = "Update title"
# 拡張子を省略してhspy形式で保存
s.save("output")
保存したファイルを開いて、タイトルの変更が反映されていることを確認してみます。
import hyperspy.api as hs
s = hs.load("output.hspy")
s.metadata
ここからは、HyperSpyで取り扱うデータについて、さらに深堀りしていきます。
HyperSpyのデータには、「シグナル」と「ナビゲーション」という概念があります。
まずは、これらの概念についてご説明します。
シグナル
ナビゲーション
hs.load()を利用してロードしたデータは、全て「BaseSigalクラス」を継承したサブクラスに保管されます。
BaseSignalクラスには、様々なサブクラスがあります。
サブクラスは、階層構造になっています。
サブクラスは、データの持つ以下の要素によって自動的に選択されます。
HyperSpyのデータをprint関数で出力すると、シグナルやナビゲーションの情報が分かります。
次のコードは、データ「core_shell.hdf5」の情報を表示します。まずは実行してみてください。
(非推奨のメタデータを使っていることを示すWarningが出ることもありますが、ここでは無視します。)
import hyperspy.api as hs
s = hs.load("core_shell.hdf5")
print(s)
出力結果は、次の意味となります。
<シグナルクラス, title: データタイトル, dimensions: (ナビゲーションの次元数|シグナルの次元数)>
シグナルクラスは、次の情報から「EDSTEMSpectrum」と判定されています。
BaseSignal subclass | signal_dimension | signal_type | dtype |
---|---|---|---|
BaseSignal | - | - | real |
Signal1D | 1 | - | real |
EELSSpectrum | 1 | EELS | real |
EDSSEMSpectrum | 1 | EDS_SEM | real |
EDSTEM | 1 | EDS_TEM | real |
Signal2D | 2 | - | real |
HologramImage | 2 | hologram | real |
DielecticFuntion | 1 | - | complex |
ComplexSignal | - | - | complex |
ComplexSignal1D | 1 | - | complex |
Complex2D | 2 | - | complex |
ナビゲーションの次元数は2次元であることが分かります。
使用可能なサブクラスは、hs.print_known_signal_types()で確認できます。
import hyperspy.api as hs
hs.print_known_signal_types()
なお、HyperSpyの拡張パッケージをインストールすることで、サブクラスを追加することもできます。
先ほど情報を確認したEDS TEMのスペクトルデータを、plotでグラフ表示してみます。
import hyperspy.api as hs
s = hs.load("core_shell.hdf5")
s.plot()
ナビゲーションの右上をよく見てください。赤い四角があるのが分かるでしょうか。
この赤い四角は、マウスでドラッグすることで場所を変更できます。赤い四角がある場所に対応するシグナルデータが、下半分にグラフで描画されます。
・・・なのですが、最初にご説明したとおり、残念ながらGoogle Colaboratoryではこのようなインタラクティブな操作はできません。
講師の画面でお見せしますので、そちらをご覧ください。
今度は、1次元のナビゲーションを持つデータを見てみたいと思います。
HyperSpyの公式サンプルには該当するデータがないため、ここでは、適当なデータをNumPy配列から作成してみます。
次の構文でNumPy配列からシグナルクラスののオブジェクトを作成することができます。
hs.signals.クラス名(NumPy配列)
次のコードは、ランダムな値で作成したNumPy配列から、Signal1Dサブクラスのオブジェクトを作成します。
import numpy as np
import hyperspy.api as hs
s = hs.signals.Signal2D(np.random.random_sample((10, 20, 100)))
print(s)
ナビゲーションが1次元、シグナルが2次元のデータが作成できました。
これをプロットしてみましょう。
import numpy as np
import hyperspy.api as hs
s = hs.signals.Signal2D(np.random.random_sample((10, 20, 100)))
print(s)
s.plot()
ナビゲーションの左端に、赤い線があるのが分かるでしょうか。
この赤い線をマウスでドラッグすることで、その場所に対応するシグナルデータが、下半分に描画されます。
こちらも残念ながらGoogle Colaboratoryでは操作できませんので、講師の画面でお見せします。
ここまでは、単にグラフを表示・保存するだけでした。
ここからは応用として、スムージングやモデルフィッティングなどの、データの活用例をご紹介します。
ここからのご説明は、ソースコードを上から順番に実行する必要があるので、動かしてみる際はご注意ください。
最初に、必要なモジュールとしてHyperSpyとMatplotlibをimportします。
import hyperspy.api as hs
import matplotlib.pyplot as plt
サンプルデータとして、HyperSpyライブラリから提供されているEDS-TEM(エネルギー分散型X線分光法-透過型電子顕微鏡)を使ってみます。
# テスト用のEDS-TEMデータを読み込み
s = hs.datasets.example_signals.EDS_TEM_Spectrum()
# プロット
s.plot()
HyperSpyでは、スペクトルなどの1次元データの処理ツールが充実しています。
別のツールと組み合わせて、自分で実装しなくてもHyperSpyだけで色々と調査・解析することができます。
バックグラウンド(本来の信号以外に計測されるノイズや不要な信号)の除去モデルとして、HyperSpyでは以下の形状をサポートしています。
GUIがあれば、範囲を指定できますがGoogle Colaboratoryでは利用できないのでCUIモードで実行します。
# バックグラウンドを除去
sb = s.remove_background(signal_range=(0, 10), zero_fill=False, background_type="PowerLaw")
# 同時に描画したいときはSignal1Dのリストにする
hs.plot.plot_spectra([s, sb], legend=["org", "background remove"])
簡易的にやってみたい場合はfastオプションを使うと速く結果をみることもできます。
# Lowess data smoothing
sb.plot()
sb.smooth_lowess()
# Total variation data smoothing
sb.plot()
sb.smooth_tv()
# Savitzky-Golay filter
sb.plot()
sb.smooth_savitzky_golay()
【ワンポイント】
スムージングのデータそのものを取得したいときは、次のように指定することで取得できます。
snpm = sb.smooth_lowess(display=False)["traitsui"].get()["single_spectrum"]
GUIが使える環境だと、動的な操作でスパイクを確認しながら除去できます。 Google Colaboratoryでは表示できないので紹介だけです。
ただし、元の信号は失われてしまうのでコピーしてからスパイクの除去を実施します。
# スムーシングした信号をコピー
sbp = sb.deepcopy()
# スパイクノイズの除去
sbp.spikes_removal_tool(interactive=False, threshold=1)
hs.plot.plot_spectra([sb, sbp], legend=["before", "remove spike"])
# ピークの検出
peaks = sbp.find_peaks1D_ohaver()
print(peaks[0]["position"])
sbp.plot()
for p in peaks[0]["position"]:
plt.axvline(x=p, color="blue", linestyle="--")
components1D
components1D(特定のシグナルタイプ対象)
まずは、モデルを定義してみましょう。
# モデルの作成
m = sbp.create_model()
中身を確認してみます。
# モデルの一覧
m.components
色々なモデルがすでに定義されているようですね。
実は、HyperSpyではEDS(エネルギー分散型X線分光)スペクトルに基づいて複数のモデルを自動的に作成してくれています。
ここでは、わかりやすくするために全部のモデルを消去してから作り直してみましょう。
# モデルを全消去
m.clear()
m.components
全部消去されましたね。
では、ガウシアンモデルをひとつ追加してみます。
# ガウシアンモデルを追加
m.append(hs.model.components1D.Gaussian())
m.components
データとモデル結果を描画してみます。
plot_componentsオプションは複数のモデルがある場合に、それぞれ表示させるかどうかを決めることができます。
# プロット
m.plot(plot_components=False)
まだガウシアンモデルはパラメータをフィッティングできていないことがわかります。
モデルのパラメータを確認してみましょう。
# パラメータの表示
m.print_current_values()
では、フィッティングさせてみます。
# フィッティング
m.multifit()
m.print_current_values()
パラメータが変化しました。うまくフィッティングできたか描画してみます。
# プロット
m.plot(plot_components=False)
ひとつの山はフィッティングできたようですが、複数あるので不十分に感じます。
ここからが、HyperSpyの見せ所です。
モデルをどんどん追加していくことができます。
# ガウシアンモデルを2つ追加
for _ in range(2):
m.append(hs.model.components1D.Gaussian())
m.multifit()
m.plot(plot_components=False)
思った通りになっているでしょうか。
パラメータは直接、編集することもできますので、ある程度は自動的にしてあとは手動でも解析することもできます。
これらの方法はあとから課題に挑戦してみてください。
今回は詳しくは説明できませんでしたが、HyperSpyには機械学習アルゴリズムも豊富に提供されています。
HyperSpyのパワフルさはいかがでしたでしょうか?
これらの一連の作業が、HyperSpyだけで可能になります。
複数ライブラリを組み合わせずとも、HyperSpyだけでコードを完結させるとメンテナンス性も向上して、研究者間でのコード共有化もはかどるのではないでしょうか?