Differences between revisions 1 and 6 (spanning 5 versions)

データの読み書き

時系列データを読む場合は、基本的にgwpy.timeseriesのTimeSeries もしくはTimeSeriesDictを用いる。
単一ファイルの単一チャンネルを読むだけなら前者で十分だが、結局のところ後者の方が上位互換なので、常にそちらを使う方がおススメ。
(gwpy1.0.0では、バグにより単一チャンネルであってもTimeSeriesで複数ファイルをまとめ読みできないとの報告がある。)

frame fileを読む

柏の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     if os.path.isdir('/frame0'):        path = '/frame0/full/'
   8     elif os.path.isdir('/data/full'):   path = '/data/full/'
   9     else:        return
  10 
  11     if type(start) in (int, float):
  12         ts = int(start)
  13         te = int(end)
  14     else:
  15         ts = int(to_gps(start))
  16         te = int(to_gps(end))
  17 
  18     t0 = ts - ts%32
  19     td = te + (32 - te%32) - ts
  20     N  = int(td/32) +1
  21     sources = []
  22     for i in range(N):
  23         dname = path + str(int(t0 + i*32)//100000)
  24         fname = '/K-K1_C-' + str(int(t0 + i*32)) + '-32.gwf'
  25         source = dname +fname
  26         if os.path.isfile(source):   sources.append(source)
  27 
  28     return sources

もしくは、Kozapyのmylibにある GetFilelist(gpsstart,gpsend)を使うのも良い。

数時間に渡る連続データを読む際のメモリ管理

超高分解能なスペクトルを作りたい時など、あまりにも長時間のデータを一度に読み込もうとすると、メモリ不足で失敗する。
そのような場合は、ある程度のファイル数で区切って.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の保存

Diagguiの.xmlファイルを読む

オンサイトで仕事をすると、diagguiでスペクトルやコヒーレンスを測定することが多い。
diagguiで保存したデータ(.xmlファイル)を解析するには、.txtファイルに書き出してそれを読む方法が最も簡単だが、いちいいち手でexportしなければならないため面倒だしヒューマンエラーも起きやすい。
ddtxmlというパッケージを使えばpythonで直接.xmlファイルを読むことができる。
ただしddtxmlはドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点

KAGRA/Subgroups/PEM/PythonMemoJP/ReadFile (last edited 2024-03-11 20:31:28 by tatsuki.washimi)

-  ⇤ ← Revision 1 as of 2020-03-18 15:03:37 → 
  Size: 2642
  Editor: tatsuki.washimi
  Comment:
+   ← Revision 6 as of 2022-04-08 10:42:12 → ⇥
  Size: 5894
  Editor: tatsuki.washimi
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 4:
-もしくは[[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 10:
-柏の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 16:
-from gwpy.timeseries import TimeSeriesDict
-Line 42:
+Line 41:
+時刻を指定してファイル一覧を取得するには、以下をコピペして使えば柏でも神岡でも使える。<<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):
    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)''を使うのも良い。


=== 数時間に渡る連続データを読む際のメモリ管理 ===
超高分解能なスペクトルを作りたい時など、あまりにも長時間のデータを一度に読み込もうとすると、メモリ不足で失敗する。<<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はドキュメントが全く整備されていないので、自力でソースコード読むなりしないと使えないのが難点