|
Size: 5112
Comment:
|
← Revision 16 as of 2025-02-14 16:25:54 ⇥
Size: 4963
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 2: | Line 2: |
時系列データを読む場合は、基本的に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]]を用いる。 単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。 (gwpy1.0.0では、バグにより単一チャンネルであっても[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries|TimeSeries]]で複数ファイルをまとめ読みできないとの報告がある。) |
|
| Line 10: | Line 5: |
| 柏のKAGRAメインサーバー(以下 柏サーバー)では''/data'' 、神岡サーバーでは ''/frame0'' もしくは ''/frame1'' 以下にファイルが置いてある。 ディレクトリ構造の詳細は[[https://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=11052|JGW-G1911052]]を参照。 GPS時刻を調べるには、[[https://www.andrews.edu/~tzs/timeconv/timeconvert.php?|GPS Time Converter]]などを使うのが手っ取り早い。 |
|
| Line 14: | Line 6: |
| 例えば、2020/3/9 14:00 UTCあたりの干渉計信号(strain)とREFL定盤の加速度計の信号を読むには、 | 時系列データを読む場合は、基本的に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 16: | Line 18: |
| 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]]) |
from lalframe.utils.frtools import get_channels print(get_channels(source)) |
| Line 24: | Line 21: |
| {{{ 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) }}} 時刻を指定してファイル一覧を取得するには、以下をコピペして使えば柏でも神岡でも使える。 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): if os.path.isdir('/frame0'): path = '/frame0/full/' elif os.path.isdir('/data/full'): path = '/data/full/' else: return if type(start) in (int, float): ts = int(start) te = int(end) else: ts = int(to_gps(start)) te = int(to_gps(end)) t0 = ts - ts%32 td = te + (32 - te%32) - ts N = int(td/32) +1 sources = [] for i in range(N): dname = path + str(int(t0 + i*32)//100000) fname = '/K-K1_C-' + str(int(t0 + i*32)) + '-32.gwf' source = dname +fname if os.path.isfile(source): sources.append(source) return sources }}} もしくは、[[https://github.com/gw-detchar/Kozapy/blob/master/samples/mylib/mylib.py|Kozapyのmylib]]にある'' GetFilelist(gpsstart,gpsend)''を使うのも良い。 |
とすればよい。前者の方が高速なのでおススメ。 |
| Line 77: | Line 24: |
| === 数時間に渡る連続データを読む際のメモリ管理 === 超高分解能なスペクトルを作りたい時など、あまりにも長時間のデータを一度に読み込もうとすると、メモリ不足で失敗する。<<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 101: | Line 28: |
== 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()} }}} |
データの読み書き
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()}