Differences between revisions 12 and 16 (spanning 4 versions)

データの読み書き

frame fileを読む

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

Jupyterのサンプルコード (要JGWdocアカウント)

frame fileに入っているチャンネル名を調べるには、コマンドラインで

$ FrChannels {frame fileのpath} | grep {探したいチャンネル名に含まれる文字列}

を使うか、pythonで

   1 from lalframe.utils.frtools import get_channels
   2 print(get_channels(source))

とすればよい。前者の方が高速なのでおススメ。

FrequencySeriesの保存

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()}

KAGRA/Subgroups/PEM/PythonMemoJP/ReadFile (last edited 2025-02-14 16:25:54 by tatsuki.washimi)

-  ⇤ ← Revision 12 as of 2024-03-11 20:31:28 → 
  Size: 6847
  Editor: tatsuki.washimi
  Comment:
+   ← Revision 16 as of 2025-02-14 16:25:54 → ⇥
  Size: 4963
  Editor: tatsuki.washimi
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 2:
-== frame fileに入っているチャンネル名を調べる ==
+== 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で
-Line 7:
+Line 21:
+とすればよい。前者の方が高速なのでおススメ。
-Line 8:
+Line 23:
-== frame fileを読む ==
-Line 10:
+Line 24:
-時系列データを読む場合は、基本的にgwpy.timeseriesの[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries|TimeSeries]]
もしくは[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeriesDict.html#gwpy.timeseries.TimeSeriesDict|TimeSeriesDict]]を用いる。<<BR>>
単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。<<BR>>
(gwpy1.0.0では、バグにより単一チャンネルであっても[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries|TimeSeries]]で複数ファイルをまとめ読みできないとの報告がある。)

柏のKAGRAメインサーバー(以下 柏サーバー)では''/data'' 、神岡サーバーでは ''/frame0'' もしくは ''/frame1'' 以下にファイルが置いてある。<<BR>>
ディレクトリ構造の詳細は[[https://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=11052|JGW-G1911052]]を参照。<<BR>>
GPS時刻を調べるには、[[https://www.andrews.edu/~tzs/timeconv/timeconvert.php?|GPS Time Converter]]などを使うのが手っ取り早い。

例えば、2020/3/9 14:00 UTCあたりの干渉計信号(strain)とREFL定盤の加速度計の信号を読むには、
{{{#!python
from gwpy.timeseries import TimeSeriesDict
import numpy as np
sources = ['/data/full/12677/K-K1_C-1267797600-32.gwf', '/data/full/12677/K-K1_C-1267797632-32.gwf', '/data/full/12677/K-K1_C-1267797664-32.gwf']
channels = ['K1:CAL-CS_PROC_C00_STRAIN_DBL_DQ', 'K1:PEM-ACC_MCF_TABLE_REFL_Z_OUT_DQ']
data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', pad=np.nan, nproc=4)
print(data[channels[0]])
print(data[channels[1]])
}}}
{{{
TimeSeries([ -2.13041590e-12,  -2.12980671e-12,  -2.12920164e-12,
            ...,   3.02297020e-13,   3.03052655e-13,
              3.03810103e-13]
           unit: dimensionless,
           t0: 1267797600.0 s,
           dt: 6.103515625e-05 s,
           name: K1:CAL-CS_PROC_C00_STRAIN_DBL_DQ,
           channel: K1:CAL-CS_PROC_C00_STRAIN_DBL_DQ)
TimeSeries([-0.00049336,  0.00021882, -0.00037257, ...,
             0.0007291 , -0.00051094,  0.00017448]
           unit: dimensionless,
           t0: 1267797600.0 s,
           dt: 0.00048828125 s,
           name: K1:PEM-ACC_MCF_TABLE_REFL_Z_OUT_DQ,
           channel: K1:PEM-ACC_MCF_TABLE_REFL_Z_OUT_DQ)
}}}

時刻を指定してファイル一覧を取得するには、以下をコピペして使うとラク。<<BR>>
GPS時刻でも、''sources = get_sources('Mar 9 2020 14:00:00', 'Mar 9 2020 14:05:00')''みたいにUTC時刻でもイケる。
{{{#!python
import os
import numpy as np
from gwpy.timeseries import TimeSeriesDict
from gwpy.time import to_gps
                                                                                  
def get_sources(start, end):
    path = '/data/KAGRA/raw/full/'
    
    if type(start) in (int, float):
        ts = int(start)
        te = int(end)
    else:
        ts = int(to_gps(start))
        te = int(to_gps(end))

    gps = np.arange(ts-ts%32, te+1, 32, dtype=int)
    sources = [path + str(t//100000) + '/K-K1_C-{}-32.gwf'.format(t)  for t in gps]
    sources = [s for s in sources if os.path.isfile(s)]

    return sources
}}}
もしくは、[[https://github.com/gw-detchar/Kozapy/blob/master/samples/mylib/mylib.py|Kozapyのmylib]]にある'' GetFilelist(gpsstart,gpsend)''を使うのも良い。

=== nds2 ===
[[https://dac.icrr.u-tokyo.ac.jp/KAGRA/DAWG/DMG/Manuals/nds2/nds2-enduser]]
=== 数時間に渡る連続データを読む際のメモリ管理 ===
超高分解能なスペクトルを作りたい時など、あまりにも長時間のデータを一度に読み込もうとすると、メモリ不足で失敗する。<<BR>>
そのような場合は、ある程度のファイル数で区切って[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.resample|.resample()]]しながら継接ぎすればよい。<<BR>>
重いオブジェクトはその都度[[https://algorithm.joho.info/programming/python/del-py/|del]]しよう。

{{{#!python
ts = float(開始GPS時刻)
te = ts + 2**10 #~17分

for i in tqdm(range(26)):  #~7時間半弱
   sources = get_sources(ts, te)
   read = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', pad=np.nan, nproc=32)
   read = read[channels[0]].resample(1024)
   if i==0 : data = read.copy()
   else:     data.append(read)

   del read, sources
   ts = ts + 2**10
   te = te + 2**10
}}}
-Line 104:
+Line 33:
-[[https://git.ligo.org/cds/dttxml|ddtxml]]というパッケージを使えばpythonで直接.xmlファイルを読むことができる。<<BR>>
ただしddtxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。<<BR>>
+[[https://git.ligo.org/cds/dttxml|dttxml]]というパッケージを使えばpythonで直接.xmlファイルを読むことができる。<<BR>>
ただしdttxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。<<BR>>
-Line 120:
+Line 49:
-=== "Power Spectrum" を読んでgwpyのFrequencySeriesにする ===
Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。<<BR>>
次の関数を定義して引数にファイル名を入れれば、ASDのdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える
+=== "Power Spectrum" 等を読んでgwpyのFrequencySeriesにする ===
Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。また、gwpyのcoherenceはいわゆる2乗コヒーレンスだが、Diagguiのcoherenceは2乗ではないものである。<<BR>>
次の関数を定義して引数にファイル名を入れれば、ASD等のdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える。<<BR>>
-Line 125:
+Line 54:
-def get_asd(xml_file):
    data = dttxml.DiagAccess(xml_file).results.PSD
    channels = data.keys()
    
    ASD ={}
    for ch in channels:
        ASD[ch] = FrequencySeries(data[ch].PSD[0], frequencies=data[ch].FHz, channel=data[ch].channelA) 
    
    return ASD
+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()}