3934
Comment:
|
← Revision 12 as of 2024-03-11 20:31:28 ⇥
6847
|
Deletions are marked like this. | Additions are marked like this. |
Line 2: | Line 2: |
== frame fileに入っているチャンネル名を調べる == {{{#!python from lalframe.utils.frtools import get_channels print(get_channels(source)) }}} == frame fileを読む == |
|
Line 4: | Line 11: |
もしくは[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeriesDict.html#gwpy.timeseries.TimeSeriesDict|TimeSeriesDict]]を用いる。 単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。 |
もしくは[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeriesDict.html#gwpy.timeseries.TimeSeriesDict|TimeSeriesDict]]を用いる。<<BR>> 単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。<<BR>> |
Line 8: | Line 15: |
== frame fileを読む == 柏のKAGRAメインサーバー(以下 柏サーバー)では''/data'' 、神岡サーバーでは ''/frame0'' もしくは ''/frame1'' 以下にファイルが置いてある。 ディレクトリ構造の詳細は[[https://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=11052|JGW-G1911052]]を参照。 |
柏のKAGRAメインサーバー(以下 柏サーバー)では''/data'' 、神岡サーバーでは ''/frame0'' もしくは ''/frame1'' 以下にファイルが置いてある。<<BR>> ディレクトリ構造の詳細は[[https://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=11052|JGW-G1911052]]を参照。<<BR>> |
Line 42: | Line 47: |
時刻を指定してファイル一覧を取得するには、以下をコピペして使えば柏でも神岡でも使える。 | 時刻を指定してファイル一覧を取得するには、以下をコピペして使うとラク。<<BR>> |
Line 51: | Line 56: |
if os.path.isdir('/frame0'): path = '/frame0/full/' elif os.path.isdir('/data/full'): path = '/data/full/' else: return |
path = '/data/KAGRA/raw/full/' |
Line 62: | Line 65: |
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) |
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)] |
Line 75: | Line 72: |
=== 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 }}} == 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|ddtxml]]というパッケージを使えばpythonで直接.xmlファイルを読むことができる。<<BR>> ただしddtxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。<<BR>> ファイルへアクセスする方法は {{{#!python acc = dttxml.DiagAccess('fname.xml') }}} とのことなので、適当なファイルを開いて {{{#!python print(vars(acc)) }}} や {{{#!python acc.__dict__.keys() }}} で調べていくのでも良い。 === "Power Spectrum" を読んでgwpyのFrequencySeriesにする === Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。<<BR>> 次の関数を定義して引数にファイル名を入れれば、ASDのdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える {{{#!python 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 }}} |
データの読み書き
frame fileに入っているチャンネル名を調べる
frame fileを読む
時系列データを読む場合は、基本的にgwpy.timeseriesのTimeSeries もしくはTimeSeriesDictを用いる。
単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。
(gwpy1.0.0では、バグにより単一チャンネルであってもTimeSeriesで複数ファイルをまとめ読みできないとの報告がある。)
柏のKAGRAメインサーバー(以下 柏サーバー)では/data 、神岡サーバーでは /frame0 もしくは /frame1 以下にファイルが置いてある。
ディレクトリ構造の詳細はJGW-G1911052を参照。
GPS時刻を調べるには、GPS Time Converterなどを使うのが手っ取り早い。
例えば、2020/3/9 14:00 UTCあたりの干渉計信号(strain)とREFL定盤の加速度計の信号を読むには、
1 from gwpy.timeseries import TimeSeriesDict
2 import numpy as np
3 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']
4 channels = ['K1:CAL-CS_PROC_C00_STRAIN_DBL_DQ', 'K1:PEM-ACC_MCF_TABLE_REFL_Z_OUT_DQ']
5 data = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', pad=np.nan, nproc=4)
6 print(data[channels[0]])
7 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)
時刻を指定してファイル一覧を取得するには、以下をコピペして使うとラク。
GPS時刻でも、sources = get_sources('Mar 9 2020 14:00:00', 'Mar 9 2020 14:05:00')みたいにUTC時刻でもイケる。
1 import os
2 import numpy as np
3 from gwpy.timeseries import TimeSeriesDict
4 from gwpy.time import to_gps
5
6 def get_sources(start, end):
7 path = '/data/KAGRA/raw/full/'
8
9 if type(start) in (int, float):
10 ts = int(start)
11 te = int(end)
12 else:
13 ts = int(to_gps(start))
14 te = int(to_gps(end))
15
16 gps = np.arange(ts-ts%32, te+1, 32, dtype=int)
17 sources = [path + str(t//100000) + '/K-K1_C-{}-32.gwf'.format(t) for t in gps]
18 sources = [s for s in sources if os.path.isfile(s)]
19
20 return sources
もしくは、Kozapyのmylibにある GetFilelist(gpsstart,gpsend)を使うのも良い。
nds2
https://dac.icrr.u-tokyo.ac.jp/KAGRA/DAWG/DMG/Manuals/nds2/nds2-enduser
数時間に渡る連続データを読む際のメモリ管理
超高分解能なスペクトルを作りたい時など、あまりにも長時間のデータを一度に読み込もうとすると、メモリ不足で失敗する。
そのような場合は、ある程度のファイル数で区切って.resample()しながら継接ぎすればよい。
重いオブジェクトはその都度delしよう。
1 ts = float(開始GPS時刻)
2 te = ts + 2**10 #~17分
3
4 for i in tqdm(range(26)): #~7時間半弱
5 sources = get_sources(ts, te)
6 read = TimeSeriesDict.read(sources, channels, format='gwf.lalframe', pad=np.nan, nproc=32)
7 read = read[channels[0]].resample(1024)
8 if i==0 : data = read.copy()
9 else: data.append(read)
10
11 del read, sources
12 ts = ts + 2**10
13 te = te + 2**10
FrequencySeriesの保存
HDF5がおすすめ。 https://gwpy.github.io/docs/stable/spectrum/io.html#gwpy-frequencyseries-io-hdf5
Diagguiの.xmlファイルを読む
オンサイトで仕事をすると、diagguiでスペクトルやコヒーレンスを測定することが多い。
diagguiで保存したデータ(.xmlファイル)を解析するには、.txtファイルに書き出してそれを読む方法が最も簡単だが、いちいいち手でexportしなければならないため面倒だしヒューマンエラーも起きやすい。
ddtxmlというパッケージを使えばpythonで直接.xmlファイルを読むことができる。
ただしddtxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点。
ファイルへアクセスする方法は
1 acc = dttxml.DiagAccess('fname.xml')
とのことなので、適当なファイルを開いて
1 print(vars(acc))
や
1 acc.__dict__.keys()
で調べていくのでも良い。
"Power Spectrum" を読んでgwpyのFrequencySeriesにする
Diagguiの"Power Spectrum"というのは、用語の間違えがそのまま放置されていて、実際にはASDであることに注意。
次の関数を定義して引数にファイル名を入れれば、ASDのdictができる。channel名を指定すればgwpyのFrequencySeriesとして扱える