コンテンツにスキップ

API Reference

波形データの読み込み

main

parse(file_name, file_name2=None, file_name3=None, fmt='vector', gain=None, bit=None, noheader=False)

強震動データを読み込み、適切な形式で読み込みます。

引数：

名前 タイプ デスクリプション デフォルト
file_name str

読み込むファイルのパス

 必須
file_name2 str

追加のファイルパス（SAC形式用など）．ほとんどの場合は不要

 None
file_name3 str

追加のファイルパス（SAC形式用など）．ほとんどの場合は不要

 None
fmt str

データのフォーマット ("vector", "nied", "jma", "win", "gns", "sac")．デフォルトはvecotrs形式（"vector"）

 'vector'
gain float

振幅倍率（win形式のみ）

 None
bit int

データビット数（win形式のみ）

 None
noheader bool

ヘッダーがあるかどうか（vectors形式のみ）

 False

戻り値：

名前 タイプ デスクリプション
vectors

読み込まれたデータの vectors オブジェクト
ソースコード位置： src/PySGM/main.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def parse(file_name,file_name2=None,file_name3=None,fmt="vector",
    gain=None,bit=None,noheader=False):
    """強震動データを読み込み、適切な形式で読み込みます。

    Args:
        file_name (str): 読み込むファイルのパス
        file_name2 (str, optional): 追加のファイルパス（SAC形式用など）．ほとんどの場合は不要
        file_name3 (str, optional): 追加のファイルパス（SAC形式用など）．ほとんどの場合は不要
        fmt (str, optional): データのフォーマット ("vector", "nied", "jma", "win", "gns", "sac")．デフォルトはvecotrs形式（"vector"）
        gain (float, optional): 振幅倍率（win形式のみ）
        bit (int, optional): データビット数（win形式のみ）
        noheader (bool, optional): ヘッダーがあるかどうか（vectors形式のみ）

    Returns:
        vectors: 読み込まれたデータの vectors オブジェクト
    """
    if fmt == "vector":
        v = vector.parse(file_name,noheader)
        return v

    elif fmt == "nied":
        file_basename,ext = os.path.splitext(file_name)
        v = nied.parse(file_basename,ext)
        return v

    elif fmt == "jma":
        file_basename,ext = os.path.splitext(file_name)
        if ext == ".csv":
            v = jma.csv_parse(file_name)
        else:
            v = jma.parse(file_name)
        return v

    elif fmt == "win":
        v = win.parse(file_name,gain,bit)
        return v

    elif fmt == "gns":
        v = gns.parse(file_name)
        return v

    elif fmt == "sac":
        v = sac.parse(file_name,file_name2,file_name3)
        return v

波形解析

このモジュールには、データの読み込み関数と、波形を扱うための主要クラスが含まれています．

Functions (関数)

parse(input_file, noheader=False) 

指定されたファイルを読み込み、3成分一括解析用の vectors オブジェクトを生成する

この関数は内部で vectors.input を呼び出す．主に PySGM 独自のテキスト形式（vector形式）を読み込む際に使用する．

引数：

名前 タイプ デスクリプション デフォルト
input_file str

読み込むファイルのパス

 必須
noheader bool

ファイルにヘッダー行がない場合は True を指定します

 False

戻り値：

名前 タイプ デスクリプション
vectors

読み込まれたデータを持つ vectors クラスのインスタンス
ソースコード位置： src/PySGM/vector.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def parse(input_file,noheader=False):
    """
        指定されたファイルを読み込み、3成分一括解析用の vectors オブジェクトを生成する

        この関数は内部で vectors.input を呼び出す．主に PySGM 独自のテキスト形式（vector形式）を読み込む際に使用する．

    Args:
        input_file (str): 読み込むファイルのパス
        noheader (bool): ファイルにヘッダー行がない場合は True を指定します

    Returns:
        vectors: 読み込まれたデータを持つ vectors クラスのインスタンス
    """
    v = vectors.input(input_file,noheader)
    return v

Classes (クラス)

vector 

単一成分の波形データを扱う基本クラス。

属性：

名前 タイプ デスクリプション
header dict

観測点情報や記録時間などを含むヘッダー
tim ndarray

時間軸データ（秒）
wave ndarray

波形データ（加速度など）
dt float

サンプリング間隔（秒）
ソースコード位置： src/PySGM/vector.py 
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
class vector:
    """
        単一成分の波形データを扱う基本クラス。

    Attributes:
        header (dict): 観測点情報や記録時間などを含むヘッダー
        tim (numpy.ndarray): 時間軸データ（秒）
        wave (numpy.ndarray): 波形データ（加速度など）
        dt (float): サンプリング間隔（秒）
    """
    def __init__(self,header,tim,wave):
        """
            vector クラスの初期化

        Args:
            header (dict): ヘッダー情報
            tim (numpy.ndarray): 時間軸データ
            wave (numpy.ndarray): 波形データ
        """
        self.header = header
        self.tim = tim
        self.wave = wave
        self.dt = tim[1] - tim[0]

    def copy(self):
        v2 = copy.deepcopy(self)
        return v2

    #----------------------------------------------#
    #  Wave Analysis
    #----------------------------------------------#
    def integration(wave,dt,low,high):
        """
            フーリエ変換を用いた数値積分

        Args:
            wave (numpy.ndarray): 積分対象の波形データ
            dt (float): サンプリング間隔（秒）
            low (float): 低域カットオフ周波数(Hz)
            high (float): 高域カットオフ周波数(Hz)

        Returns:
            numpy.ndarray: 積分後の波形データ
        """
        w = np.fft.fft(wave)
        freq = np.fft.fftfreq(len(wave),d=dt)
        df = freq[1] - freq[0]

        nt = 10
        low0  = max(0,low - nt*df)
        high0 = high + nt*df

        w2 = np.ones_like(w)
        for (i,f) in enumerate(freq):
            coef = (0.0+1.0j)*2.0*math.pi*f

            if abs(f) <= low0:
                w2[i] = 0.0 + 0.0j
            elif abs(f) < low:
                w2[i] = w[i] * (abs(f)-low0)/(low-low0) / coef
            elif abs(f) <= high:
                w2[i] = w[i] / coef
            elif abs(f) <= high0:
                w2[i] = w[i] * (abs(f)-high0)/(high-high0) / coef
            else:
                w2[i] = 0.0 + 0.0j

        wave_int = np.real(np.fft.ifft(w2))
        return wave_int

    def differentiation(wave,dt,low,high):
        w = np.fft.fft(wave)
        freq = np.fft.fftfreq(len(wave),d=dt)
        df = freq[1] - freq[0]

        nt = 10
        low0  = max(0,low - nt*df)
        high0 = high + nt*df

        w2 = np.ones_like(w)
        for (i,f) in enumerate(freq):
            coef = (0.0+1.0j)*2.0*math.pi*f

            if abs(f) < low0:
                w2[i] = 0.0 + 0.0j
            elif abs(f) < low:
                w2[i] = w[i] * (abs(f)-low0)/(low-low0) * coef
            elif abs(f) <= high:
                w2[i] = w[i] * coef
            elif abs(f) <= high0:
                w2[i] = w[i] * (abs(f)-high0)/(high-high0) * coef
            else:
                w2[i] = 0.0 + 0.0j

        wave_int = np.real(np.fft.ifft(w2))
        return wave_int

    def bandpass(wave,dt,low,high):
        w = np.fft.fft(wave)
        freq = np.fft.fftfreq(len(wave),d=dt)
        df = freq[1] - freq[0]

        nt = 10
        low0  = max(0,low - nt*df)
        high0 = high + nt*df

        w2 = np.ones_like(w)
        for (i,f) in enumerate(freq):
            if abs(f) < low0:
                w2[i] = 0.0 + 0.0j
            elif abs(f) < low:
                w2[i] = w[i] * (abs(f)-low0)/(low-low0)
            elif abs(f) <= high:
                w2[i] = w[i]
            elif abs(f) <= high0:
                w2[i] = w[i] * (abs(f)-high0)/(high-high0)
            else:
                w2[i] = 0.0 + 0.0j

        wave_int = np.real(np.fft.ifft(w2))
        return wave_int

    def bandpass_butter(wave,dt,low,high,order=5):
        fs = 1.0/dt
        nyq = 0.5*fs
        low_cut = low/nyq
        high_cut = high/nyq
        b,a = signal.butter(order,[low_cut,high_cut],btype='band')
        wave2 = signal.lfilter(b,a,wave)
        return wave2

    def direct_integration(wave,dt):
        wave_int = wave.cumsum()*dt
        return wave_int

    # ---------------------------------------------------    
    def boore_integration(acc,dt,t0):
        def boore_coeff(vel,n0,n1):
            ntim = len(vel)
            x = np.array(range(ntim))
            a0 = (vel[n0:n1].T @ (x[n0:n1]-n0)) / ((x[n0:n1]-n0).T @ (x[n0:n1]-n0))
            a1 = ((vel[n1:]-a0*(n1-n0)).T @ (x[n1:]-n1)) / ((x[n1:]-n1).T @ (x[n1:]-n1))
            return a0, a1

        def boore_residual(vel,n0,n1,a0,a1):
            x = range(len(vel))
            res = np.sum(vel[0:n0]**2)
            res += np.sum((vel[n0:n1] - a0*(range(n1-n0)))**2)
            res += np.sum((vel[n1:] - (a1*(range(len(vel)-n1)) + a0*(n1-n0)))**2)
            return res

        vel = vector.direct_integration(acc,dt)
        ntim = len(vel)
        n0_ini = int(t0/dt)

        for n0 in range(n0_ini,n0_ini+6000,100):
            for n in range(100,3000,10):
                n1 = n0 + n
                a0,a1 = boore_coeff(vel,n0,n1)
                res = boore_residual(vel,n0,n1,a0,a1)

                if ('res_min' not in locals()) or (res < res_min):
                    print(n0,n1,res)
                    res_min = res
                    n0_min = n0
                    n1_min = n1

        a0,a1 = boore_coeff(vel,n0_min,n1_min)
        x = np.array(range(ntim))
        y0 = a0*(x[n0_min:n1_min]-n0_min)
        y1 = a1*(x[n1_min:]-n1_min) + a0*(n1_min-n0_min)

        vel[n0_min:n1_min] -= y0
        vel[n1_min:] -= y1
        disp = vector.direct_integration(vel,dt)

        return disp

    # ---------------------------------------------------    
    def envelope_vector(wave):
        wave_ana = signal.hilbert(wave)
        wave_env = np.abs(wave_ana)
        return wave_env

    def envelope(self):
        v2 = copy.deepcopy(self)
        wave_ana = signal.hilbert(self.wave)
        v2.wave = np.abs(wave_ana)
        return v2

    def power_spectrum(wave,fs):
        freq, P = signal.periodogram(wave,fs)
        return freq, P

    def fourier_spectrum(wave,fs):
        def set_parameters(wave,dt):
            n = len(wave) / 2
            ntw = 2
            while ntw < n:
                ntw = ntw*2
            freq = abs(np.fft.fftfreq(ntw,d=dt))
            return freq, ntw

        dt = 1.0/float(fs)
        freq, ntw = set_parameters(wave,dt)
        w = np.fft.fft(wave[0:ntw])

        return freq, w

    def peak_ground(wave,print_result=True):
        PG = math.sqrt(np.max(wave**2))

        if print_result:
            print("PGX: ",PG)

        return PG

    def normalize(wave,scale=1.0):
        PG = math.sqrt(np.max(wave**2))
        return wave/PG*scale

    def roll(self,shift):
        nshift = int(shift/self.dt)
        v2 = copy.deepcopy(self)
        v2.wave = np.roll(self.wave,nshift)
        return v2

    def amplify(self,amp):
        v2 = copy.deepcopy(self)
        v2.wave = self.wave*amp
        return v2

    def sum(self,v1):
        v2 = copy.deepcopy(self)
        v2.wave += v1.wave
        return v2

    def extract(self,start,end):
        ns = int(start/self.dt)
        ne = int(end/self.dt)
        ntim = ne-ns
        v2 = copy.deepcopy(self)
        v2.wave = self.wave[ns:ne]
        v2.tim = self.tim[ns:ne]
        v2.header['ntim'] = ntim
        return v2

    #----------------------------------------------#
    def coherence(self,wave,window):
        freq,f0 = vector.fourier_spectrum(self.wave,1/self.dt)
        freq,f1 = vector.fourier_spectrum(wave.wave,1/self.dt)

        p0 = np.conjugate(f0)*f0
        p1 = np.conjugate(f1)*f1
        c01 = np.conjugate(f0)*f1

        nw = int(window/(freq[1]-freq[0]))
        if nw%2 == 0:
            nw += 1
        nw2 = int((nw-1)/2)
        w = signal.parzen(nw)

        a = np.r_[c01[nw2:0:-1],c01,c01[0],c01[-1:-nw2:-1]]
        c01_s = np.convolve(w/w.sum(),a,mode='valid')

        a = np.r_[p0[nw2:0:-1],p0,p0[0],p0[-1:-nw2:-1]]
        p0_s = np.convolve(w/w.sum(),a,mode='valid')
        a = np.r_[p1[nw2:0:-1],p1,p1[0],p1[-1:-nw2:-1]]
        p1_s = np.convolve(w/w.sum(),a,mode='valid')

        c01_abs = np.conjugate(c01_s)*c01_s

        return freq, c01_abs/(p0_s*p1_s)

    #----------------------------------------------#
    #  Residual functions
    #----------------------------------------------#
    def nonlinear_residual(self,v):
        v0 = np.sum(self.wave**2)
        v1 = np.sum(v.wave**2)
        v01 = np.sum((self.wave-v.wave)**2)
        return v01*v01/v0/v1

    def linear_residual(self,v):
        v0 = np.sum(self.wave**2)
        v01 = np.sum((self.wave-v.wave)**2)
        return v01/v0

    #----------------------------------------------#
    #  Plot functions
    #----------------------------------------------#
    def plot(self):

        start = self.tim[0]
        end = self.tim[-1]

        ns = int(start/self.dt)
        ne = int(end/self.dt)

        fig = plt.figure()

        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        ax1 = fig.add_subplot(111)

        if 'code' in self.header and 'record_time' in self.header:
            ax1.set_title(self.header['code']+" "+self.header['record_time'])

        ax1.plot(self.tim,self.wave,color='k',lw=1)

        amax = np.amax(np.abs(self.wave[ns:ne]))

        ax1.set_xlim(start,end)
        ax1.set_ylim(-1.2*amax,1.2*amax)
        ax1.set_xlabel("time [s]")

        plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0, hspace=0)
        plt.show()

    def plot_with(self,v):
        start = self.tim[0]
        end = min(self.tim[-1],v.tim[-1])

        ns = 0
        ne = int((end-start)/self.dt)

        fig = plt.figure()

        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        ax1 = fig.add_subplot(211)
        ax2 = fig.add_subplot(212)

        if 'code' in self.header and 'record_time' in self.header:
            ax1.set_title(self.header['code']+" "+self.header['record_time'])

        ax1.plot(self.tim,self.wave,color='b',lw=1)
        ax2.plot(v.tim,v.wave,color='k',lw=1)

        amax = np.amax(np.abs(self.wave[ns:ne]))

        ax1.set_xlim(start,end)
        ax2.set_xlim(start,end)

        ax1.set_ylim(-1.2*amax,1.2*amax)
        ax2.set_ylim(-1.2*amax,1.2*amax)

        ax2.set_xlabel("time [s]")

        plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0, hspace=0)
        plt.show()

__init__(header, tim, wave) 

vector クラスの初期化

引数：

名前 タイプ デスクリプション デフォルト
header dict

ヘッダー情報

 必須
tim ndarray

時間軸データ

 必須
wave ndarray

波形データ

 必須
ソースコード位置： src/PySGM/vector.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
def __init__(self,header,tim,wave):
    """
        vector クラスの初期化

    Args:
        header (dict): ヘッダー情報
        tim (numpy.ndarray): 時間軸データ
        wave (numpy.ndarray): 波形データ
    """
    self.header = header
    self.tim = tim
    self.wave = wave
    self.dt = tim[1] - tim[0]

integration(wave, dt, low, high) 

フーリエ変換を用いた数値積分

引数：

名前 タイプ デスクリプション デフォルト
wave ndarray

積分対象の波形データ

 必須
dt float

サンプリング間隔（秒）

 必須
low float

低域カットオフ周波数(Hz)

 必須
high float

高域カットオフ周波数(Hz)

 必須

戻り値：

タイプ デスクリプション

numpy.ndarray: 積分後の波形データ
ソースコード位置： src/PySGM/vector.py 
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def integration(wave,dt,low,high):
    """
        フーリエ変換を用いた数値積分

    Args:
        wave (numpy.ndarray): 積分対象の波形データ
        dt (float): サンプリング間隔（秒）
        low (float): 低域カットオフ周波数(Hz)
        high (float): 高域カットオフ周波数(Hz)

    Returns:
        numpy.ndarray: 積分後の波形データ
    """
    w = np.fft.fft(wave)
    freq = np.fft.fftfreq(len(wave),d=dt)
    df = freq[1] - freq[0]

    nt = 10
    low0  = max(0,low - nt*df)
    high0 = high + nt*df

    w2 = np.ones_like(w)
    for (i,f) in enumerate(freq):
        coef = (0.0+1.0j)*2.0*math.pi*f

        if abs(f) <= low0:
            w2[i] = 0.0 + 0.0j
        elif abs(f) < low:
            w2[i] = w[i] * (abs(f)-low0)/(low-low0) / coef
        elif abs(f) <= high:
            w2[i] = w[i] / coef
        elif abs(f) <= high0:
            w2[i] = w[i] * (abs(f)-high0)/(high-high0) / coef
        else:
            w2[i] = 0.0 + 0.0j

    wave_int = np.real(np.fft.ifft(w2))
    return wave_int

vectors

Bases: vector

3成分（EW, NS, UD）の強震動データを一括して管理・解析するためのクラス

属性：

名前 タイプ デスクリプション
header dict

観測点情報や記録時間などを含むヘッダー
tim ndarray

時間軸データ（秒）
ew ndarray

EW成分の波形データ
ns ndarray

NS成分の波形データ
ud ndarray

UD成分の波形データ
dt float

サンプリング間隔（秒）
ソースコード位置： src/PySGM/vector.py 
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
class vectors(vector):
    """
        3成分（EW, NS, UD）の強震動データを一括して管理・解析するためのクラス

    Attributes:
        header (dict): 観測点情報や記録時間などを含むヘッダー
        tim (numpy.ndarray): 時間軸データ（秒）
        ew (numpy.ndarray): EW成分の波形データ
        ns (numpy.ndarray): NS成分の波形データ
        ud (numpy.ndarray): UD成分の波形データ
        dt (float): サンプリング間隔（秒）
    """
    def __init__(self,header,tim,ew,ns,ud):
        """vectors クラスの初期化

        Args:
            header (dict): ヘッダー情報
            tim (numpy.ndarray): 時間軸データ
            ew (numpy.ndarray): EW成分のデータ
            ns (numpy.ndarray): NS成分のデータ
            ud (numpy.ndarray): UD成分のデータ
        """
        self.header = header
        self.tim = tim
        self.ew = ew
        self.ns = ns
        self.ud = ud
        self.dt = tim[1] - tim[0]

    def copy(self):
        v2 = copy.deepcopy(self)
        return v2

    #----------------------------------------------#
    #  Construction functions
    #----------------------------------------------#
    def append(self,v):
        """vectors オブジェクトを後ろに繋げる

        Args:
            v (vectors): 後ろに追加する vectors オブジェクト

        Returns:
            vectors: 後ろにvが追加された新しい vectors インスタンス
        """
        v2 = self
        v2.header = self.header
        v2.ew = np.hstack((self.ew,v.ew))
        v2.ns = np.hstack((self.ns,v.ns))
        v2.ud = np.hstack((self.ud,v.ud))

        dt = self.tim[1] - self.tim[0]
        shift = np.full_like(v.tim,self.tim[-1]+dt)
        tim = v.tim + shift

        v2.tim = np.hstack((self.tim,tim))
        v2.header['ntim'] = len(v2.tim)

        return v2

    def to_vector(self,comp):
        if comp == "ew":
            v = vector(self.header,self.tim,self.ew)
        elif comp == "ns":
            v = vector(self.header,self.tim,self.ns)
        elif comp == "ud":
            v = vector(self.header,self.tim,self.ud)
        return v

    #----------------------------------------------#
    #  Correction functions
    #----------------------------------------------#
    def trend_removal(self,sec=None):
        """
            3成分（EW, NS, UD）すべてに対して、線形トレンドの除去（基線補正）を一括適用する

        Args:
            なし

        Returns:
            なし（インスタンス内の各成分が補正後の波形に更新される）
        """

        if sec is None:
            self.ew = self.ew - np.average(self.ew)
            self.ns = self.ns - np.average(self.ns)
            self.ud = self.ud - np.average(self.ud)
        else:
            n = int(sec/self.dt)
            self.ew = self.ew - np.average(self.ew[:n])
            self.ns = self.ns - np.average(self.ns[:n])
            self.ud = self.ud - np.average(self.ud[:n])

    def rotation(self,rot,deg='deg'):
        """
            水平2成分（EW, NS）を指定した角度だけ回転させる

        Args:
            rot (float): 回転角
            deg (str): 角度の単位。'deg'（度）または 'rad'（ラジアン）

        Returns:
            vectors: 回転後の成分を持つ新しい vectors インスタンス
        """
        if deg == 'deg':
            rad = rot/180.0 * math.pi
        else:
            rad = rot

        p = self.ew * math.sin(rad) + self.ns * math.cos(rad)
        n = self.ew * math.cos(rad) - self.ns * math.sin(rad)

        v2 = copy.deepcopy(self)
        v2.ew = n
        v2.ns = p

        return v2

    def zero_padding(self,ntim):
        self.trend_removal()
        ntim0 = len(self.tim)

        if ntim0 < ntim:
            zero = np.zeros(ntim-ntim0,dtype=float)

            v2 = self
            v2.header = self.header
            v2.ew = np.hstack((self.ew,zero))
            v2.ns = np.hstack((self.ns,zero))
            v2.ud = np.hstack((self.ud,zero))

            dt = self.tim[1] - self.tim[0]
            v2.tim = np.linspace(0,ntim*dt,ntim,endpoint=False)
            v2.header['ntim'] = len(v2.tim)

        else:
            v2 = copy.deepcopy(self)

        return v2

    def time_shift(self,shift):
        nshift = int(shift/self.dt)
        v2 = copy.deepcopy(self)
        v2.ew = np.roll(self.ew,-nshift)
        v2.ns = np.roll(self.ns,-nshift)
        v2.ud = np.roll(self.ud,-nshift)
        return v2

    def normalize(self,scale=1.0):
        v2 = copy.deepcopy(self)
        v2.ew = vector.normalize(self.ew,scale)
        v2.ns = vector.normalize(self.ns,scale)
        v2.ud = vector.normalize(self.ud,scale)
        return v2

    def level_shift(self,level):
        v2 = copy.deepcopy(self)
        v2.ew = self.ew + level
        v2.ns = self.ns + level
        v2.ud = self.ud + level
        return v2

    def trim(self,ntim):
        v2 = copy.deepcopy(self)
        v2.ew = self.ew[0:ntim]
        v2.ns = self.ns[0:ntim]
        v2.ud = self.ud[0:ntim]
        v2.tim = self.tim[0:ntim]
        v2.header['ntim'] = ntim

        return v2

    def adjust(self,v):
        ntim = v.header['ntim']

        if self.header['ntim'] > ntim:
            v2 = self.trim(ntim)
        else:
            v2 = self.zero_padding(ntim)

        return v2

    def every(self,n):
        v2 = copy.deepcopy(self)
        v2.ew = self.ew[0::n]
        v2.ns = self.ns[0::n]
        v2.ud = self.ud[0::n]
        v2.tim = self.tim[0::n]
        v2.dt = v2.tim[1] - v2.tim[0]

        return v2

    def resampling(self,fs):
        f = int(1.0/self.dt)
        if f > fs:
            n = int(f/fs)
            v2 = self.every(n)
        else:
            v2 = self

        return v2

    def amplify(self,amp):
        v2 = copy.deepcopy(self)
        v2.ew = self.ew*amp
        v2.ns = self.ns*amp
        v2.ud = self.ud*amp
        return v2

    def extract(self,start,end=60,to_end=True):
        ns = int(start/self.dt)
        if to_end:
            ne = self.header['ntim']
        else:
            ne = int(end/self.dt)
        ntim = ne-ns

        v2 = copy.deepcopy(self)
        v2.ew = self.ew[ns:ne]
        v2.ns = self.ns[ns:ne]
        v2.ud = self.ud[ns:ne]
        v2.tim = self.tim[ns:ne]
        v2.header['ntim'] = ntim
        return v2

    def reset_tim(self):
        t0 = self.tim[0]
        self.tim = self.tim - t0
        self.header['ntim'] = len(self.tim)

    #----------------------------------------------#
    def set_origin_time(self,origin_time_str):
        origin_time = datetime.datetime.strptime(origin_time_str,"%Y/%m/%d %H:%M:%S")
        record_time = datetime.datetime.strptime(self.header['record_time'],"%Y/%m/%d %H:%M:%S")
        dt = record_time - origin_time
        self.tim = self.tim + dt.total_seconds()

    #----------------------------------------------#
    #  Wave Analysis
    #----------------------------------------------#
    def integration(self,low=0.2,high=50):
        """
            3成分（EW, NS, UD）すべてに対して、フーリエ変換を用いて数値積分する．指定した周波数帯域（low - high）以外はカットされる．

        Args:
            low (float): 低域カットオフ周波数（Hz）．デフォルトは 0.2
            high (float): 高域カットオフ周波数（Hz）．デフォルトは 50

        Returns:
            なし（インスタンス内の self.ew, self.ns, self.ud が積分後の波形に更新される）
        """
        v2 = copy.deepcopy(self)
        v2.ew = vector.integration(self.ew,self.dt,low,high)
        v2.ns = vector.integration(self.ns,self.dt,low,high)
        v2.ud = vector.integration(self.ud,self.dt,low,high)

        return v2

    def differentiation(self,low=0.05,high=50):
        v2 = copy.deepcopy(self)
        v2.ew = vector.differentiation(self.ew,self.dt,low,high)
        v2.ns = vector.differentiation(self.ns,self.dt,low,high)
        v2.ud = vector.differentiation(self.ud,self.dt,low,high)

        return v2

    def bandpass(self,low=0.05,high=10):
        """
            3成分（EW, NS, UD）すべてに対して、バンドパスフィルタを一括適用する

        Args:
            low (float): 低域カットオフ周波数（Hz）
            high (float): 高域カットオフ周波数（Hz）

        Returns:
            なし（インスタンス内の各成分がフィルタ適用後の波形に更新される）
        """
        v2 = copy.deepcopy(self)
        v2.ew = vector.bandpass(self.ew,self.dt,low,high)
        v2.ns = vector.bandpass(self.ns,self.dt,low,high)
        v2.ud = vector.bandpass(self.ud,self.dt,low,high)

        return v2

    def bandpass_butter(self,low=0.05,high=10,order=5):
        v2 = copy.deepcopy(self)
        v2.ew = vector.bandpass_butter(self.ew,self.dt,low,high,order)
        v2.ns = vector.bandpass_butter(self.ns,self.dt,low,high,order)
        v2.ud = vector.bandpass_butter(self.ud,self.dt,low,high,order)
        return v2

    def direct_integration(self):
        v2 = copy.deepcopy(self)
        v2.ew = vector.direct_integration(self.ew,self.dt)
        v2.ns = vector.direct_integration(self.ns,self.dt)
        v2.ud = vector.direct_integration(self.ud,self.dt)
        return v2

    def boore_integration(self,t0=15.0):
        v2 = copy.deepcopy(self)
        v2.trend_removal(sec=t0)
        v2.ew = vector.boore_integration(self.ew,self.dt,t0)
        v2.ns = vector.boore_integration(self.ns,self.dt,t0)
        v2.ud = vector.boore_integration(self.ud,self.dt,t0)
        return v2

    def envelope(self):
        v2 = copy.deepcopy(self)
        v2.ew = vector.envelope_vector(self.ew)
        v2.ns = vector.envelope_vector(self.ns)
        v2.ud = vector.envelope_vector(self.ud)
        return v2

    def power_spectrum(self):
        fs = 1.0/self.dt
        freq, Pew = vector.power_spectrum(self.ew,fs)
        freq, Pns = vector.power_spectrum(self.ns,fs)
        freq, Pud = vector.power_spectrum(self.ud,fs)

        self.ps = spectrum.spectrum(self.header,freq,Pew,Pns,Pud)
        return self.ps

    def fourier_spectrum(self):
        fs = 1.0/self.dt
        freq, Few = vector.fourier_spectrum(self.ew,fs)
        freq, Fns = vector.fourier_spectrum(self.ns,fs)
        freq, Fud = vector.fourier_spectrum(self.ud,fs)

        self.fs = spectrum.spectrum(self.header,freq,Few,Fns,Fud)
        return self.fs

    def fourier_amp_spectrum(self):
        fs = 1.0/self.dt
        freq, Few = vector.fourier_spectrum(self.ew,fs)
        freq, Fns = vector.fourier_spectrum(self.ns,fs)
        freq, Fud = vector.fourier_spectrum(self.ud,fs)

        self.fs = spectrum.spectrum(self.header,freq,np.abs(Few),np.abs(Fns),np.abs(Fud))
        return self.fs

    def spectrum_ratio(self,ps0):
        fs = 1.0/self.dt
        freq, Pew = vector.power_spectrum(self.ew,fs)
        freq, Pns = vector.power_spectrum(self.ns,fs)
        freq, Pud = vector.power_spectrum(self.ud,fs)

        nf = min(len(freq),len(ps0.freq))

        sr_ew = np.sqrt(Pew[0:nf] / ps0.ew[0:nf])
        sr_ns = np.sqrt(Pns[0:nf] / ps0.ns[0:nf])
        sr_ud = np.sqrt(Pud[0:nf] / ps0.ud[0:nf])

        sr = spectrum.spectrum(self.header,freq,sr_ew,sr_ns,sr_ud)
        return sr

    def response_spectrum(self):
        self.period = np.logspace(-1,1,100)

        self.Sa_ew, self.Sv_ew, self.Sd_ew = response.response_spectrum(self.ew,self.period,self.dt)
        self.Sa_ns, self.Sv_ns, self.Sd_ns = response.response_spectrum(self.ns,self.period,self.dt)
        self.Sa_ud, self.Sv_ud, self.Sd_ud = response.response_spectrum(self.ud,self.period,self.dt)

    def response_spectrum_max(self):
        """
            水平成分の応答スペクトルをRotD100で算出する．実行するとvectorsオブジェクトに以下のAttributesが追加される．

        Returns:
            None: このメソッドは値を返さない．代わりに以下の属性が追加される．
                - **Sa**: 加速度応答スペクトル（h=5%)
                - **Sv**: 速度応答スペクトル（h=5%)
                - **Sd**: 変位応答スペクトル（h=5%)
                - **Sa_rot**: 加速度応答スペクトルが最大となる方位（ラジアン）
                - **Sv_rot**: 速度応答スペクトルが最大となる方位（ラジアン）
                - **Sd_rot**: 変位応答スペクトルが最大となる方位（ラジアン）
        """
        self.period = np.logspace(-1,1,100)

        self.Sa, self.Sv, self.Sd, self.Sa_rot, self.Sv_rot, self.Sd_rot = response.response_spectrum_max(self.ew,self.ns,self.period,self.dt)

    def response_spectrum_peak_period(self):
        self.period = np.logspace(-1,1,100)
        pSv,peak_pSv,peak_period = response.response_spectrum_pseudo_FD(self.ew,self.ns,self.period,self.dt)

        return peak_pSv, peak_period

    def response_spectrum_peak_period_psv(self):
        self.period = np.logspace(-1,1,100)
        pSv,peak_pSv,peak_period = response.response_spectrum_pseudo(self.ew,self.ns,self.period,self.dt)

        return peak_pSv, peak_period

    def calc_SI(self):
        SI = response.calc_SI_FD(self.ew,self.ns,self.dt)
        return SI

    def peak_ground_3d(self,print_result=True):
        """
            3成分合成波の最大値（PGAやPGV）を算出する

        Args:
            print_result (bool): 結果を表示するかどうか

        Returns:
            float: 最大値
        """
        PG = math.sqrt(np.max(self.ew**2 + self.ns**2 + self.ud**2))

        if print_result:
            print("PGX (3d): ",PG)

        return PG

    def peak_ground_2d(self,print_result=True):
        """
            水平2成分合成波の最大値（PGAやPGV）を算出する

        Args:
            print_result (bool): 結果を表示するかどうか

        Returns:
            float: 最大値
        """
        PG = math.sqrt(np.max(self.ew**2 + self.ns**2))

        if print_result:
            print("PGX (2d): ",PG)

        return PG

    def peak_ground_all(self,print_result=True):
        PG_ew = math.sqrt(np.max(self.ew**2))
        PG_ns = math.sqrt(np.max(self.ns**2))
        PG_ud = math.sqrt(np.max(self.ud**2))

        if print_result:
            print("PGX (ew): ",PG_ew)
            print("PGX (ns): ",PG_ns)
            print("PGX (ud): ",PG_ud)

        return PG_ew, PG_ns, PG_ud

    def jma_seismic_intensity(self,print_result=True):
        """
            気象庁計測震度を計算

        Args:
            print_result (bool,optional): 結果をコンソールに表示する場合はTrue

        Returns:
            float: 計算された計測震度値
        """
        si = jsi.jsi(self.ew,self.ns,self.ud,self.dt)

        if print_result:
            print("JMA Seismic Intensity: ", si)

        return si

    def GM2025_modified_intensity(self,k=0.4,print_result=True):
        si = jsi.GM2025(self.ew,self.ns,self.ud,self.dt,k)

        if print_result:
            print("GM2025 Intensity: ", si)

        return si

    def realtime_seismic_intensity(self):
        rsi = realtime_jsi.realtime_jsi(self.ew,self.ns,self.ud,self.dt)
        return rsi

    def peak_period_Sv(self):
        self.period = np.logspace(-1,1,100)
        peak_Sv,peak_period = response.response_spectrum_peak_period_sv(self.ew,self.ns,self.period,self.dt)

        return peak_Sv, peak_period

    def demand_curve(self,ductility=4.0):
        self.period = np.logspace(-1,np.log10(5.0),25)

        self.Dc_ew = demand_curve.demand_curve(self.ew,self.period,ductility,self.dt)
        self.Dc_ns = demand_curve.demand_curve(self.ns,self.period,ductility,self.dt)

    def cumulative_square(self):
        return np.sum(self.ns**2 + self.ew**2) * self.dt

    def hv_spectrum(self,nt=4096,start=60,ncut=10,window=0.2):

        fsample = 1.0/self.dt
        nseg = int(((len(self.tim) - 30) - start)/(nt/2))

        index = start
        for i in range(0,nseg):
            ew = self.ew[index:index+nt]
            ns = self.ns[index:index+nt]
            ud = self.ud[index:index+nt]
            rms = math.sqrt(np.var(ew) + np.var(ns) + np.var(ud))

            try:
                rms_list += [rms]
                index_list += [index]
            except:
                rms_list = [rms]
                index_list = [index]

            index += int(nt/2)

        r = np.array(rms_list)
        sorted_index = np.argsort(r)

        for i in sorted_index[0:10]:
            index = index_list[i]

            ew = self.ew[index:index+nt]
            ns = self.ns[index:index+nt]
            ud = self.ud[index:index+nt]

            ew -= np.average(ew)
            ns -= np.average(ns)
            ud -= np.average(ud)

            freq, Few = vector.fourier_spectrum(ew,fsample)
            freq, Fns = vector.fourier_spectrum(ns,fsample)
            freq, Fud = vector.fourier_spectrum(ud,fsample)

            fs = spectrum.spectrum(self.header,freq,Few,Fns,Fud)
            fs_smooth = fs.smoothing(window)

            try:
                fs_list.append(fs_smooth)
            except:
                fs_list = spectrum.spectrum_list(fs_smooth)

        fs_list.average(window)

        hv_list = []
        for s in fs_list.sr:
            hv_tmp = np.sqrt(np.abs(s.ew)**2 + np.abs(s.ns)**2) / np.abs(s.ud)
            hv_list += [hv_tmp]

        hv = np.sqrt(fs_list.ave.ew**2 + fs_list.ave.ns**2) / fs_list.ave.ud

        return fs_list.freq, hv, hv_list

    #----------------------------------------------#
    #  I/O functions
    #----------------------------------------------#
    def input(input_file,noheader=False):
        with open(input_file,'r') as file:
            lines = list(file)

        if noheader:
            header = ""
            tim,ew,ns,ud = np.loadtxt(lines[0:], usecols=(0,1,2,3), unpack=True)
        else:
            h = lines[0].strip('#').rstrip()
            header = eval(h)
            tim,ew,ns,ud = np.loadtxt(lines[1:], usecols=(0,1,2,3), unpack=True)

        v = vectors(header,tim,ew,ns,ud)
        return v

    def input_pickle(input_file):
        with open(input_file,'rb') as f:
            v = pickle.load(f)
        return v

    #----------------------------------------------#
    def output(self,output_file,fmt="%15.7f"):
        """
            3成分（EW, NS, UD）の時刻歴波形をファイルに書き出すする

        Args:
            output_ffile (str): 出力先ファイル名
            fmt (str): 出力フォーマット
        """
        header = str(self.header)
        output = np.c_[self.tim,self.ew,self.ns,self.ud]
        np.savetxt(output_file,output,fmt=fmt,header=header,comments="#")

    def output_csv(self,output_file,fmt="%15.7f"):
        header = ["{}: {}".format(key,self.header[key]) for key in self.header.keys()]
        output = np.c_[self.tim,self.ew,self.ns,self.ud]
        with open(output_file,'w') as f:
            writer = csv.writer(f)
            writer.writerow(header)
            writer.writerow(["tim","ew","ns","ud"])
            writer.writerows(output)

    def output_pickle(self,output_file):
        with open(output_file,'wb') as f:
            pickle.dump(self,f)

    def output_rs(self,output_file,fmt="%15.7f"):
        header = str(self.header)

        output = np.c_[self.period,self.Sa_ew,self.Sa_ns,self.Sa_ud]
        np.savetxt(output_file+".sa",output,fmt=fmt,header=header,comments="#")

        output = np.c_[self.period,self.Sv_ew,self.Sv_ns,self.Sv_ud]
        np.savetxt(output_file+".sv",output,fmt=fmt,header=header,comments="#")

        output = np.c_[self.period,self.Sd_ew,self.Sd_ns,self.Sd_ud]
        np.savetxt(output_file+".sd",output,fmt=fmt,header=header,comments="#")

    def output_rs_max(self,output_file,fmt="%15.7f"):
        header = str(self.header)

        output = np.c_[self.period,self.Sa,self.Sa_rot]
        np.savetxt(output_file+".samax",output,fmt=fmt,header=header,comments="#")

        output = np.c_[self.period,self.Sv,self.Sv_rot]
        np.savetxt(output_file+".svmax",output,fmt=fmt,header=header,comments="#")

        output = np.c_[self.period,self.Sd,self.Sd_rot]
        np.savetxt(output_file+".sdmax",output,fmt=fmt,header=header,comments="#")


    def output_ps(self,output_file,fmt="%15.7f"):
        self.ps.output(output_file,fmt=fmt)


    def output_dc(self,output_file,fmt="%15.7f"):
        header = str(self.header)

        output = np.c_[self.period,self.Dc_ew,self.Dc_ns]
        np.savetxt(output_file+".dc",output,fmt=fmt,header=header,comments="#")

    #----------------------------------------------#
    #  Plot functions
    #----------------------------------------------#
    def plot_all(self,start=0,end=60,to_end=False,output_file=None):
        """
            3成分（EW, NS, UD）の時刻歴波形を一括してプロットします。

        Args:
            start (float): 表示の開始時刻（秒）
            end (float): 表示の終了時刻（秒）
            to_end (bool): 波形データの最後まで表示する場合はTrue
            output_file (str): 外部ファイルに画像として出力する場合は，ファイル名を指定する．ファイル名の拡張子に応じた形式で保存される．
        """
        if to_end:
            start = self.tim[0]
            end = self.tim[-1]
            amax = np.amax(np.r_[np.abs(self.ew),np.abs(self.ns),np.abs(self.ud)])
        else:
            ns = int(start/self.dt)
            ne = int(end/self.dt)
            amax = np.amax(np.r_[np.abs(self.ew[ns:ne]),np.abs(self.ns[ns:ne]),np.abs(self.ud[ns:ne])])

        fig = plt.figure()

        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        ax1 = fig.add_subplot(311)
        ax2 = fig.add_subplot(312)
        ax3 = fig.add_subplot(313)

        if 'code' in self.header and 'record_time' in self.header:
            ax1.set_title(self.header['code']+" "+self.header['record_time'])


        ax1.plot(self.tim,self.ew,color='k',lw=1)
        ax2.plot(self.tim,self.ns,color='k',lw=1)
        ax3.plot(self.tim,self.ud,color='k',lw=1)



        ax1.set_xlim(start,end)
        ax2.set_xlim(start,end)
        ax3.set_xlim(start,end)

        ax1.set_ylim(-1.2*amax,1.2*amax)
        ax2.set_ylim(-1.2*amax,1.2*amax)
        ax3.set_ylim(-1.2*amax,1.2*amax)

        ax1.set_xticklabels([])
        ax2.set_xticklabels([])

        ax3.set_xlabel("time [s]")

        ax1.text(0.02,0.85,"EW",transform=ax1.transAxes)
        ax2.text(0.02,0.85,"NS",transform=ax2.transAxes)
        ax3.text(0.02,0.85,"UD",transform=ax3.transAxes)

        plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0, hspace=0)
        if output_file is None:
            plt.show()
        else:
            plt.savefig(output_file)

    def plot_rs_all(self):
        fig = plt.figure()

        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        plt.subplot(331)
        plt.title("Sa")
        plt.plot(self.period,self.Sa_ew,color='maroon',lw=1.5)
        plt.ylabel("EW")

        plt.subplot(332)
        plt.title("Sv")
        plt.plot(self.period,self.Sv_ew,color='darkorange',lw=1.5)

        plt.subplot(333)
        plt.title("Sd")
        plt.plot(self.period,self.Sd_ew,color='darkgreen',lw=1.5)


        plt.subplot(334)
        plt.plot(self.period,self.Sa_ns,color='maroon',lw=1.5)
        plt.ylabel("NS")

        plt.subplot(335)
        plt.plot(self.period,self.Sv_ns,color='darkorange',lw=1.5)

        plt.subplot(336)
        plt.plot(self.period,self.Sd_ns,color='darkgreen',lw=1.5)


        plt.subplot(337)
        plt.plot(self.period,self.Sa_ud,color='maroon',lw=1.5)
        plt.ylabel("UD")
        plt.xlabel("period [s]")

        plt.subplot(338)
        plt.plot(self.period,self.Sv_ud,color='darkorange',lw=1.5)
        plt.xlabel("period [s]")

        plt.subplot(339)
        plt.plot(self.period,self.Sd_ud,color='darkgreen',lw=1.5)
        plt.xlabel("period [s]")

        axs = plt.gcf().get_axes()
        for ax in axs:
            plt.axes(ax)
            plt.xscale("log")
            plt.yscale("log")
            plt.xlim(0.1,10)
            ax.xaxis.set_major_formatter(FormatStrFormatter('%.1f'))
            ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))

        plt.show()

    def plot_rs_max(self):
        fig = plt.figure()

        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        plt.subplot(231)
        plt.xscale("log")
        plt.yscale("log")
        plt.title("Sa")
        plt.plot(self.period,self.Sa,color='maroon',lw=1.5)
        plt.ylabel("Sa")

        plt.subplot(232)
        plt.xscale("log")
        plt.yscale("log")
        plt.title("Sv")
        plt.plot(self.period,self.Sv,color='darkorange',lw=1.5)
        plt.ylabel("Sv")

        plt.subplot(233)
        plt.xscale("log")
        plt.yscale("log")
        plt.title("Sd")
        plt.plot(self.period,self.Sd,color='darkgreen',lw=1.5)
        plt.ylabel("Sd")


        plt.subplot(234)
        plt.xscale("log")
        plt.ylim(-90,90)
        plt.plot(self.period,self.Sa_rot,color='maroon',lw=1.5)
        plt.ylabel("azimth")
        plt.xlabel("period [s]")

        plt.subplot(235)
        plt.xscale("log")
        plt.ylim(-90,90)
        plt.plot(self.period,self.Sv_rot,color='darkorange',lw=1.5)
        plt.xlabel("period [s]")

        plt.subplot(236)
        plt.xscale("log")
        plt.ylim(-90,90)
        plt.plot(self.period,self.Sd_rot,color='darkgreen',lw=1.5)
        plt.xlabel("period [s]")

        plt.show()

    def plot_ps_all(self):
        self.ps.plot_ps_all()

__init__(header, tim, ew, ns, ud)

vectors クラスの初期化

引数：

名前 タイプ デスクリプション デフォルト
header dict

ヘッダー情報

 必須
tim ndarray

時間軸データ

 必須
ew ndarray

EW成分のデータ

 必須
ns ndarray

NS成分のデータ

 必須
ud ndarray

UD成分のデータ

 必須
ソースコード位置： src/PySGM/vector.py 
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
def __init__(self,header,tim,ew,ns,ud):
    """vectors クラスの初期化

    Args:
        header (dict): ヘッダー情報
        tim (numpy.ndarray): 時間軸データ
        ew (numpy.ndarray): EW成分のデータ
        ns (numpy.ndarray): NS成分のデータ
        ud (numpy.ndarray): UD成分のデータ
    """
    self.header = header
    self.tim = tim
    self.ew = ew
    self.ns = ns
    self.ud = ud
    self.dt = tim[1] - tim[0]

append(v)

vectors オブジェクトを後ろに繋げる

引数：

名前 タイプ デスクリプション デフォルト
v vectors

後ろに追加する vectors オブジェクト

 必須

戻り値：

名前 タイプ デスクリプション
vectors

後ろにvが追加された新しい vectors インスタンス
ソースコード位置： src/PySGM/vector.py 
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
def append(self,v):
    """vectors オブジェクトを後ろに繋げる

    Args:
        v (vectors): 後ろに追加する vectors オブジェクト

    Returns:
        vectors: 後ろにvが追加された新しい vectors インスタンス
    """
    v2 = self
    v2.header = self.header
    v2.ew = np.hstack((self.ew,v.ew))
    v2.ns = np.hstack((self.ns,v.ns))
    v2.ud = np.hstack((self.ud,v.ud))

    dt = self.tim[1] - self.tim[0]
    shift = np.full_like(v.tim,self.tim[-1]+dt)
    tim = v.tim + shift

    v2.tim = np.hstack((self.tim,tim))
    v2.header['ntim'] = len(v2.tim)

    return v2

bandpass(low=0.05, high=10) 

3成分（EW, NS, UD）すべてに対して、バンドパスフィルタを一括適用する

引数：

名前 タイプ デスクリプション デフォルト
low float

低域カットオフ周波数（Hz）

 0.05
high float

高域カットオフ周波数（Hz）

 10

戻り値：

タイプ デスクリプション

なし（インスタンス内の各成分がフィルタ適用後の波形に更新される）
ソースコード位置： src/PySGM/vector.py 
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
def bandpass(self,low=0.05,high=10):
    """
        3成分（EW, NS, UD）すべてに対して、バンドパスフィルタを一括適用する

    Args:
        low (float): 低域カットオフ周波数（Hz）
        high (float): 高域カットオフ周波数（Hz）

    Returns:
        なし（インスタンス内の各成分がフィルタ適用後の波形に更新される）
    """
    v2 = copy.deepcopy(self)
    v2.ew = vector.bandpass(self.ew,self.dt,low,high)
    v2.ns = vector.bandpass(self.ns,self.dt,low,high)
    v2.ud = vector.bandpass(self.ud,self.dt,low,high)

    return v2

integration(low=0.2, high=50) 

3成分（EW, NS, UD）すべてに対して、フーリエ変換を用いて数値積分する．指定した周波数帯域（low - high）以外はカットされる．

引数：

名前 タイプ デスクリプション デフォルト
low float

低域カットオフ周波数（Hz）．デフォルトは 0.2

 0.2
high float

高域カットオフ周波数（Hz）．デフォルトは 50

 50

戻り値：

タイプ デスクリプション

なし（インスタンス内の self.ew, self.ns, self.ud が積分後の波形に更新される）
ソースコード位置： src/PySGM/vector.py 
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
def integration(self,low=0.2,high=50):
    """
        3成分（EW, NS, UD）すべてに対して、フーリエ変換を用いて数値積分する．指定した周波数帯域（low - high）以外はカットされる．

    Args:
        low (float): 低域カットオフ周波数（Hz）．デフォルトは 0.2
        high (float): 高域カットオフ周波数（Hz）．デフォルトは 50

    Returns:
        なし（インスタンス内の self.ew, self.ns, self.ud が積分後の波形に更新される）
    """
    v2 = copy.deepcopy(self)
    v2.ew = vector.integration(self.ew,self.dt,low,high)
    v2.ns = vector.integration(self.ns,self.dt,low,high)
    v2.ud = vector.integration(self.ud,self.dt,low,high)

    return v2

jma_seismic_intensity(print_result=True) 

気象庁計測震度を計算

引数：

名前 タイプ デスクリプション デフォルト
print_result (bool, optional)

結果をコンソールに表示する場合はTrue

 True

戻り値：

名前 タイプ デスクリプション
float

計算された計測震度値
ソースコード位置： src/PySGM/vector.py 
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
def jma_seismic_intensity(self,print_result=True):
    """
        気象庁計測震度を計算

    Args:
        print_result (bool,optional): 結果をコンソールに表示する場合はTrue

    Returns:
        float: 計算された計測震度値
    """
    si = jsi.jsi(self.ew,self.ns,self.ud,self.dt)

    if print_result:
        print("JMA Seismic Intensity: ", si)

    return si

output(output_file, fmt='%15.7f') 

3成分（EW, NS, UD）の時刻歴波形をファイルに書き出すする

引数：

名前 タイプ デスクリプション デフォルト
output_ffile str

出力先ファイル名

 必須
fmt str

出力フォーマット

 '%15.7f'
ソースコード位置： src/PySGM/vector.py 
961
962
963
964
965
966
967
968
969
970
971
def output(self,output_file,fmt="%15.7f"):
    """
        3成分（EW, NS, UD）の時刻歴波形をファイルに書き出すする

    Args:
        output_ffile (str): 出力先ファイル名
        fmt (str): 出力フォーマット
    """
    header = str(self.header)
    output = np.c_[self.tim,self.ew,self.ns,self.ud]
    np.savetxt(output_file,output,fmt=fmt,header=header,comments="#")

peak_ground_2d(print_result=True) 

水平2成分合成波の最大値（PGAやPGV）を算出する

引数：

名前 タイプ デスクリプション デフォルト
print_result bool

結果を表示するかどうか

 True

戻り値：

名前 タイプ デスクリプション
float

最大値
ソースコード位置： src/PySGM/vector.py 
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
def peak_ground_2d(self,print_result=True):
    """
        水平2成分合成波の最大値（PGAやPGV）を算出する

    Args:
        print_result (bool): 結果を表示するかどうか

    Returns:
        float: 最大値
    """
    PG = math.sqrt(np.max(self.ew**2 + self.ns**2))

    if print_result:
        print("PGX (2d): ",PG)

    return PG

peak_ground_3d(print_result=True) 

3成分合成波の最大値（PGAやPGV）を算出する

引数：

名前 タイプ デスクリプション デフォルト
print_result bool

結果を表示するかどうか

 True

戻り値：

名前 タイプ デスクリプション
float

最大値
ソースコード位置： src/PySGM/vector.py 
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
def peak_ground_3d(self,print_result=True):
    """
        3成分合成波の最大値（PGAやPGV）を算出する

    Args:
        print_result (bool): 結果を表示するかどうか

    Returns:
        float: 最大値
    """
    PG = math.sqrt(np.max(self.ew**2 + self.ns**2 + self.ud**2))

    if print_result:
        print("PGX (3d): ",PG)

    return PG

plot_all(start=0, end=60, to_end=False, output_file=None) 

3成分（EW, NS, UD）の時刻歴波形を一括してプロットします。

引数：

名前 タイプ デスクリプション デフォルト
start float

表示の開始時刻（秒）

 0
end float

表示の終了時刻（秒）

 60
to_end bool

波形データの最後まで表示する場合はTrue

 False
output_file str

外部ファイルに画像として出力する場合は，ファイル名を指定する．ファイル名の拡張子に応じた形式で保存される．

 None
ソースコード位置： src/PySGM/vector.py 
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
def plot_all(self,start=0,end=60,to_end=False,output_file=None):
    """
        3成分（EW, NS, UD）の時刻歴波形を一括してプロットします。

    Args:
        start (float): 表示の開始時刻（秒）
        end (float): 表示の終了時刻（秒）
        to_end (bool): 波形データの最後まで表示する場合はTrue
        output_file (str): 外部ファイルに画像として出力する場合は，ファイル名を指定する．ファイル名の拡張子に応じた形式で保存される．
    """
    if to_end:
        start = self.tim[0]
        end = self.tim[-1]
        amax = np.amax(np.r_[np.abs(self.ew),np.abs(self.ns),np.abs(self.ud)])
    else:
        ns = int(start/self.dt)
        ne = int(end/self.dt)
        amax = np.amax(np.r_[np.abs(self.ew[ns:ne]),np.abs(self.ns[ns:ne]),np.abs(self.ud[ns:ne])])

    fig = plt.figure()

    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    ax1 = fig.add_subplot(311)
    ax2 = fig.add_subplot(312)
    ax3 = fig.add_subplot(313)

    if 'code' in self.header and 'record_time' in self.header:
        ax1.set_title(self.header['code']+" "+self.header['record_time'])


    ax1.plot(self.tim,self.ew,color='k',lw=1)
    ax2.plot(self.tim,self.ns,color='k',lw=1)
    ax3.plot(self.tim,self.ud,color='k',lw=1)



    ax1.set_xlim(start,end)
    ax2.set_xlim(start,end)
    ax3.set_xlim(start,end)

    ax1.set_ylim(-1.2*amax,1.2*amax)
    ax2.set_ylim(-1.2*amax,1.2*amax)
    ax3.set_ylim(-1.2*amax,1.2*amax)

    ax1.set_xticklabels([])
    ax2.set_xticklabels([])

    ax3.set_xlabel("time [s]")

    ax1.text(0.02,0.85,"EW",transform=ax1.transAxes)
    ax2.text(0.02,0.85,"NS",transform=ax2.transAxes)
    ax3.text(0.02,0.85,"UD",transform=ax3.transAxes)

    plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0, hspace=0)
    if output_file is None:
        plt.show()
    else:
        plt.savefig(output_file)

response_spectrum_max() 

水平成分の応答スペクトルをRotD100で算出する．実行するとvectorsオブジェクトに以下のAttributesが追加される．

戻り値：

名前 タイプ デスクリプション
None

このメソッドは値を返さない．代わりに以下の属性が追加される． - Sa: 加速度応答スペクトル（h=5%) - Sv: 速度応答スペクトル（h=5%) - Sd: 変位応答スペクトル（h=5%) - Sa_rot: 加速度応答スペクトルが最大となる方位（ラジアン） - Sv_rot: 速度応答スペクトルが最大となる方位（ラジアン） - Sd_rot: 変位応答スペクトルが最大となる方位（ラジアン）
ソースコード位置： src/PySGM/vector.py 
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
def response_spectrum_max(self):
    """
        水平成分の応答スペクトルをRotD100で算出する．実行するとvectorsオブジェクトに以下のAttributesが追加される．

    Returns:
        None: このメソッドは値を返さない．代わりに以下の属性が追加される．
            - **Sa**: 加速度応答スペクトル（h=5%)
            - **Sv**: 速度応答スペクトル（h=5%)
            - **Sd**: 変位応答スペクトル（h=5%)
            - **Sa_rot**: 加速度応答スペクトルが最大となる方位（ラジアン）
            - **Sv_rot**: 速度応答スペクトルが最大となる方位（ラジアン）
            - **Sd_rot**: 変位応答スペクトルが最大となる方位（ラジアン）
    """
    self.period = np.logspace(-1,1,100)

    self.Sa, self.Sv, self.Sd, self.Sa_rot, self.Sv_rot, self.Sd_rot = response.response_spectrum_max(self.ew,self.ns,self.period,self.dt)

rotation(rot, deg='deg') 

水平2成分（EW, NS）を指定した角度だけ回転させる

引数：

名前 タイプ デスクリプション デフォルト
rot float

回転角

 必須
deg str

角度の単位。'deg'（度）または 'rad'（ラジアン）

 'deg'

戻り値：

名前 タイプ デスクリプション
vectors

回転後の成分を持つ新しい vectors インスタンス
ソースコード位置： src/PySGM/vector.py 
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
def rotation(self,rot,deg='deg'):
    """
        水平2成分（EW, NS）を指定した角度だけ回転させる

    Args:
        rot (float): 回転角
        deg (str): 角度の単位。'deg'（度）または 'rad'（ラジアン）

    Returns:
        vectors: 回転後の成分を持つ新しい vectors インスタンス
    """
    if deg == 'deg':
        rad = rot/180.0 * math.pi
    else:
        rad = rot

    p = self.ew * math.sin(rad) + self.ns * math.cos(rad)
    n = self.ew * math.cos(rad) - self.ns * math.sin(rad)

    v2 = copy.deepcopy(self)
    v2.ew = n
    v2.ns = p

    return v2

trend_removal(sec=None) 

3成分（EW, NS, UD）すべてに対して、線形トレンドの除去（基線補正）を一括適用する

戻り値：

タイプ デスクリプション

なし（インスタンス内の各成分が補正後の波形に更新される）
ソースコード位置： src/PySGM/vector.py 
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
def trend_removal(self,sec=None):
    """
        3成分（EW, NS, UD）すべてに対して、線形トレンドの除去（基線補正）を一括適用する

    Args:
        なし

    Returns:
        なし（インスタンス内の各成分が補正後の波形に更新される）
    """

    if sec is None:
        self.ew = self.ew - np.average(self.ew)
        self.ns = self.ns - np.average(self.ns)
        self.ud = self.ud - np.average(self.ud)
    else:
        n = int(sec/self.dt)
        self.ew = self.ew - np.average(self.ew[:n])
        self.ns = self.ns - np.average(self.ns[:n])
        self.ud = self.ud - np.average(self.ud[:n])