= データの読み書き =


== frame fileを読む ==

時系列データを読む場合は、基本的にgwpy.timeseriesの[[https://gwpy.github.io/docs/stable/timeseries/|TimeSeries]]
もしくは[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeriesDict/|TimeSeriesDict]]を用いる。<<BR>>
単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。<<BR>>

[[https://gwdoc.icrr.u-tokyo.ac.jp/DocDB/0165/T2516508/001/read_gwf.html|Jupyterのサンプルコード]] (要JGWdocアカウント)

frame fileに入っているチャンネル名を調べるには、コマンドラインで
{{{
$ FrChannels {frame fileのpath} | grep {探したいチャンネル名に含まれる文字列}
}}}
を使うか、pythonで
{{{#!python
from lalframe.utils.frtools import get_channels
print(get_channels(source))
}}}
とすればよい。前者の方が高速なのでおススメ。



== FrequencySeriesの保存 ==
HDF5がおすすめ。
[[https://gwpy.github.io/docs/stable/spectrum/io.html#gwpy-frequencyseries-io-hdf5]]


== Diagguiの.xmlファイルを読む ==
オンサイトで仕事をすると、diagguiでスペクトルやコヒーレンスを測定することが多い。<<BR>>
diagguiで保存したデータ(.xmlファイル)を解析するには、.txtファイルに書き出してそれを読む方法が最も簡単だが、いちいいち手でexportしなければならないため面倒だしヒューマンエラーも起きやすい。<<BR>>
[[https://git.ligo.org/cds/dttxml|dttxml]]というパッケージを使えばpythonで直接.xmlファイルを読むことができる。<<BR>>
ただしdttxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。<<BR>>
ファイルへアクセスする方法は
{{{#!python
acc = dttxml.DiagAccess('fname.xml')
}}}
とのことなので、適当なファイルを開いて
{{{#!python
print(vars(acc))
}}}
や
{{{#!python
acc.__dict__.keys()
}}}
で調べていくのでも良い。

=== "Power Spectrum" 等を読んでgwpyのFrequencySeriesにする ===
Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。また、gwpyのcoherenceはいわゆる2乗コヒーレンスだが、Diagguiのcoherenceは2乗ではないものである。<<BR>>
次の関数を定義して引数にファイル名を入れれば、ASD等のdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える。<<BR>>

{{{#!python
import dttxml
from gwpy.frequencyseries import FrequencySeries

def get_ASDs(file):
    diaggui = dttxml.DiagAccess(file).results.PSD
    return {ch: FrequencySeries(diaggui[ch].PSD[0], frequencies=diaggui[ch].FHz, 
                                epoch=diaggui[ch].gps_second,
                                channel=diaggui[ch].channelA,
                                name   =diaggui[ch].channelA) 
            for ch in diaggui.keys()}

def get_COHs(file):
    diaggui = dttxml.DiagAccess(file).results.COH      
    return {chA: {diaggui[chA].channelB[i]:
                  FrequencySeries(diaggui[chA].coherence[i], frequencies=diaggui[chA].FHz, 
                                  epoch=diaggui[chA].gps_second,
                                  channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
                                  name   =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
                  for i in range(len(diaggui[chA].channelB))}
            for chA in diaggui.keys()}

def get_CSDs(file):
    diaggui = dttxml.DiagAccess(file).results.CSD
    return {chA: {diaggui[chA].channelB[i]:
                  FrequencySeries(diaggui[chA].CSD[i], frequencies=diaggui[chA].FHz, 
                                  epoch=diaggui[chA].gps_second,
                                  channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
                                  name   =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
                  for i in range(len(diaggui[chA].channelB))}
            for chA in diaggui.keys()}


def get_TFs(file):
    diaggui = dttxml.DiagAccess(file).results._mydict['TF']
    return {chA: {diaggui[chA].channelB[i]:
                  FrequencySeries(diaggui[chA].xfer[i], frequencies=diaggui[chA].FHz, 
                                  epoch=diaggui[chA].gps_second,
                                  channel=diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA,
                                  name   =diaggui[chA].channelB[i]+' / '+diaggui[chA].channelA)
                  for i in range(len(diaggui[chA].channelB))}
            for chA in diaggui.keys()}
}}}