データの読み書き
frame fileを読む
時系列データを読む場合は、基本的にgwpy.timeseriesのTimeSeries もしくはTimeSeriesDictを用いる。
単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。
Jupyterのサンプルコード (要JGWdocアカウント)
frame fileに入っているチャンネル名を調べるには、コマンドラインで
$ FrChannels {frame fileのpath} | grep {探したいチャンネル名に含まれる文字列}を使うか、pythonで
とすればよい。前者の方が高速なのでおススメ。
FrequencySeriesの保存
HDF5がおすすめ。 https://gwpy.github.io/docs/stable/spectrum/io.html#gwpy-frequencyseries-io-hdf5
Diagguiの.xmlファイルを読む
オンサイトで仕事をすると、diagguiでスペクトルやコヒーレンスを測定することが多い。
diagguiで保存したデータ(.xmlファイル)を解析するには、.txtファイルに書き出してそれを読む方法が最も簡単だが、いちいいち手でexportしなければならないため面倒だしヒューマンエラーも起きやすい。
dttxmlというパッケージを使えばpythonで直接.xmlファイルを読むことができる。
ただしdttxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。
ファイルへアクセスする方法は
1 acc = dttxml.DiagAccess('fname.xml')
とのことなので、適当なファイルを開いて
1 print(vars(acc))
や
1 acc.__dict__.keys()
で調べていくのでも良い。
"Power Spectrum" 等を読んでgwpyのFrequencySeriesにする
Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。また、gwpyのcoherenceはいわゆる2乗コヒーレンスだが、Diagguiのcoherenceは2乗ではないものである。
次の関数を定義して引数にファイル名を入れれば、ASD等のdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える。
1 import dttxml
2 from gwpy.frequencyseries import FrequencySeries
3
4 def get_ASDs(file):
5 diaggui = dttxml.DiagAccess(file).results.PSD
6 return {ch: FrequencySeries(diaggui[ch].PSD[0], frequencies=diaggui[ch].FHz,
7 epoch=diaggui[ch].gps_second,
8 channel=diaggui[ch].channelA,
9 name =diaggui[ch].channelA)
10 for ch in diaggui.keys()}
11
12 def get_COHs(file):
13 diaggui = dttxml.DiagAccess(file).results.COH
14 return {chA: {diaggui[chA].channelB[i]:
15 FrequencySeries(diaggui[chA].coherence[i], frequencies=diaggui[chA].FHz,
16 epoch=diaggui[chA].gps_second,
17 channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
18 name =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
19 for i in range(len(diaggui[chA].channelB))}
20 for chA in diaggui.keys()}
21
22 def get_CSDs(file):
23 diaggui = dttxml.DiagAccess(file).results.CSD
24 return {chA: {diaggui[chA].channelB[i]:
25 FrequencySeries(diaggui[chA].CSD[i], frequencies=diaggui[chA].FHz,
26 epoch=diaggui[chA].gps_second,
27 channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
28 name =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
29 for i in range(len(diaggui[chA].channelB))}
30 for chA in diaggui.keys()}
31
32
33 def get_TFs(file):
34 diaggui = dttxml.DiagAccess(file).results._mydict['TF']
35 return {chA: {diaggui[chA].channelB[i]:
36 FrequencySeries(diaggui[chA].xfer[i], frequencies=diaggui[chA].FHz,
37 epoch=diaggui[chA].gps_second,
38 channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
39 name =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
40 for i in range(len(diaggui[chA].channelB))}
41 for chA in diaggui.keys()}