pythonでトリパタイト（三軸図）を描く

地震波の特徴の把握や建物の応答を推定するためにトリパタイトがよく用いられますが、ふと考えてみると、加速度応答スペクトルや速度応答スペクトル、変位応答スペクトルはExcel等を利用して描いたことがありますが、トリパタイトは自分で描いたことが無いことに気付きました。そこで、今回プログラムを作成して実際に描いてみることにしました。プログラム言語には今回はpython3を使用することにします。

## 目次

 

## トリパタイトとは

ここでは、加速度応答スペクトル、速度応答スペクトル、変位応答スペクトルを１つの図に描いたものを指します。１つの図でこれらを表現するために３つの軸が存在することから三軸図、あるいはトリパタイト（Tripartite）と呼ばれます。一般的に、横軸は固有周期（あるいは固有振動数）、縦軸は速度応答、右上がりの軸が加速度応答、右下がりの軸が変位応答として表現されることが多いように思います。このとき、例えば速度応答スペクトルをSvとすると、実用上、変位応答スペクトルはSd≒Sv/ω、加速度応答スペクトルはSa≒Sv*ωとして定義され、それぞれ疑似変位応答スペクトル、疑似加速度応答スペクトルと呼ばれます（ここで、ω=2π/T）。１つの図で加速度、速度、変位の最大応答値を確認できるため、地震波に対する構造物の応答推定において非常に有用です。

 

## 描いてみる

では早速プログラムを作成していきたいと思います。

 

### 要件定義

早速プログラム作成、、、とは言っても、まずは要件定義が必要です。今から作成するプログラムはどのように運用されるものか、持つべき機能は何かなどといったことを考えます。今回は以下のように考えました。

* GUIは用意せず、コマンドプロンプト等から実行する
* 応答スペクトル計算は別プログラムに任せて、描画するだけのシンプルなプログラムとする
* いくつかのスペクトルを重ね描けるようにする
* 図の各種設定（表示範囲や線色など）をユーザ指定できるようにする
* 固有周期の線を重ね描けるようにする

以上の要件を基にプログラムを考えていきます。

 

### 構成
要件をふまえて今回は以下のような構成を考えました。

プログラムが読み込むデータは大きく２つです。１つは設定ファイルです。ユーザが各種設定を行えるように設定ファイル（config.json）を用意します。このファイルを用意することでユーザが指定できる項目を増やしたり、今後の拡張にも対応しやすくなります。json形式としたのは、人間にもプログラムにも読みやすいデータ形式であるためです。もう１つは応答スペクトルファイルです（ここでは速度応答スペクトルを想定します）。今回開発するプログラムでスペクトル計算まで含めてもよかったのですが、スペクトルの計算にも様々な条件や設定があり、そこまで含めると今回着目している「トリパタイトの描画」部分が分かりづらくなってしまう可能性があったため、ここではスペクトルデータは読み込むだけとすることにしました。なお、こうすることで、応答スペクトルをどのような形で用意しても、本プログラムを利用できるというメリットもあります。読み込むデータは以上の２つで、出力データはトリパタイトの画像ファイルとします。

 

### 完成例

解説は後ほど、ということで、まずは各ファイルの内容とプログラムの全コードおよびプログラム実行結果（出力画像）を示します。

#### 速度応答スペクトルファイル(response_spectrum.last.csv)

Version,2.0.5(2019/06/24)
Date,2019/10/21 16:43:07
Title,"Model20191021_001"
*** Response Spectrum Analysis ***,CASE000001
,Period       ,Vel-EL-CENTRO_NS,Vel-TAFT_EW  ,Vel-HACHINOHE_NS
,2.00000e-002,5.40086e-002,1.63185e-002,2.10775e-002
,2.12957e-002,9.60331e-002,3.08567e-002,3.54935e-002
,2.26754e-002,1.69266e-001,5.54261e-002,6.79695e-002
,2.41444e-002,2.33818e-001,6.98767e-002,8.04928e-002
,2.57087e-002,3.35093e-001,7.95677e-002,9.49895e-002
,2.73742e-002,4.13139e-001,8.43459e-002,9.08988e-002
,2.91477e-002,5.56141e-001,1.02678e-001,1.03497e-001
,3.10360e-002,6.23290e-001,1.08571e-001,1.36316e-001
（以下省略）

#### 設定ファイル(config.json)

{
  "graph":
    {
      "output_file_name":"Tripartite.png",
      "dpi":600,
      "graph_width":5.0,
      "graph_height":6.5,
      "min_period":0.01,
      "max_period":10.0,
      "min_vel":0.1,
      "max_vel":1000.0,
      "x_label":"T (S)",
      "y_label":"Sv (cm/s)",
      "acc_label":"Sa (cm/s/s)",
      "disp_label":"Sd (cm)",
      "eigen_period":[1.5, 0.7, 0.2]
    },
  "spectra":
    [
      {
        "name":"EL-CENTRO_NS",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":3,
        "color":"deepskyblue"
      },
      {
        "name":"TAFT_EW",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":4,
        "color":"black"
      },
      {
        "name":"HACHINOHE_NS",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":5,
        "color":"red"
      }
    ]
}

#### トリパタイト描画プログラム(Tripartite.py)

# -*- coding:utf8 -*-
import json, sys, csv, math
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# グラフ設定クラス
class GraphConfig:
    def __init__(self, config:dict):
        self.FileName = config['output_file_name']
        self.Dpi = config['dpi']
        self.Width = config['graph_width']
        self.Height = config['graph_height']
        self.MinPeriod = config['min_period']
        self.MaxPeriod = config['max_period']
        self.MinVel = config['min_vel']
        self.MaxVel = config['max_vel']
        self.XLabel = config['x_label']
        self.YLabel = config['y_label']
        self.AccLabel = config['acc_label']
        self.DispLabel = config['disp_label']
        self.EigenPeriod = config['eigen_period']

# スペクトルクラス
class Spectrum:
    def __init__(self, item):
        self.Name = item['name']
        self.Path = item['path']
        self.SkipRow = item['skip_row']
        self.PeriodCol = item['period_col'] - 1
        self.SpectrumCol = item['spectrum_col'] - 1
        self.Color = item['color']

# グラフ初期化
def graph_init():
    fig = plt.figure(figsize=(graph.Width, graph.Height))
    ax = fig.add_subplot(1,1,1)
    plt.subplots_adjust(left=0.15)
    return fig, ax

# グラフフォーマット
def graph_format(ax):
    ax.set_xlim(graph.MinPeriod, graph.MaxPeriod)
    ax.set_ylim(graph.MinVel, graph.MaxVel)
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.grid(which="both")
    ax.legend()
    ax.get_xaxis().set_major_formatter(plt.FormatStrFormatter('%.2f'))
    ax.get_yaxis().set_major_formatter(plt.FormatStrFormatter('%.2f'))
    ax.set_xlabel(graph.XLabel)
    ax.set_ylabel(graph.YLabel)

# 設定ファイルを読む
def read_config(path:str):
    with open(path) as f:
        config = json.load(f)
    return config

# 速度応答スペクトルを読む
def read_vel_resp_spectrum(path:str, skip_row:int, period_col:int, spectrum_col:int):
    df = pd.read_csv(path, header=None, skiprows=skip_row)
    return df[period_col].values, df[spectrum_col].values

# 加速度と変位のグリッド線を取得する
def get_grid_lines_of_acc_and_disp():
    x1 = graph.MinPeriod
    x2 = graph.MaxPeriod
    # acc
    acc_list = np.concatenate([np.linspace(int(10**i), int(10**(i+1)), 10) for i in range(4)])
    acc_lines = {}
    for acc in acc_list:
        # Sv = Sa / (2π/T)
        y1 = acc / (2.0 * np.pi / x1)
        y2 = acc / (2.0 * np.pi / x2)
        acc_lines[acc] = [x1, x2, y1, y2]
    # disp
    disp_list = np.concatenate([np.linspace(int(0.01*10**i), int(0.01*10**(i+1)), 10) for i in range(4)])
    disp_lines = {}
    for disp in disp_list:
        # Sv = Sd * (2π/T)
        y1 = disp * (2.0 * np.pi / x1)
        y2 = disp * (2.0 * np.pi / x2)
        disp_lines[disp] = [x1, x2, y1, y2]
    return acc_lines, disp_lines

# グリッドを描画する
def draw_grid(grid_type, fig, ax, label, lines, text_pos, angle, vertalalign):
    for k, v in lines.items():
        x1, x2, y1, y2 = tuple(v)
        ax.plot([x1, x2], [y1, y2], color='darkgray', linewidth='0.5')
        text = '{:.0f}'.format(k)
        if (grid_type == 'acc'):
            if (np.log10(k)).is_integer():
                ax.text(x1, y1, text, rotation=angle, verticalalignment=vertalalign)
        elif (grid_type == 'disp'):
            if k >= 1.0 and (np.log10(k)).is_integer():
                ax.text(x2, y2, text, rotation=angle, verticalalignment=vertalalign)
    fig.text(text_pos[0], text_pos[1], label, rotation=angle)

# 加速度および変位のグリッドを描画する
def draw_grids(fig, ax):
    acc_lines, disp_lines = get_grid_lines_of_acc_and_disp()
    draw_grid('acc', fig, ax, graph.AccLabel, acc_lines, [0.3, 0.65], 45, 'bottom')
    draw_grid('disp', fig, ax, graph.DispLabel, disp_lines, [0.7, 0.75], -45, 'top')

# 速度応答スペクトルを描画する
def draw_vel_spectrum(ax):
    for spec in spectra:
        x, y = read_vel_resp_spectrum(spec.Path, spec.SkipRow, spec.PeriodCol, spec.SpectrumCol)
        ax.plot(x, y, label=spec.Name, color=spec.Color)

# 固有周期の線を描画する
def draw_period_line(ax):
    for period in graph.EigenPeriod:
        ax.plot([period, period], [graph.MinVel, graph.MaxVel], color='black', linestyle='dashed')

# トリパタイトを描画する
def draw_tripartite():
    fig, ax = graph_init()
    draw_grids(fig, ax)
    draw_vel_spectrum(ax)
    draw_period_line(ax)
    graph_format(ax)
    plt.savefig(graph.FileName, dpi=graph.Dpi)
    plt.show()

if __name__ == '__main__':
    config = read_config('./config.json')
    graph = GraphConfig(config['graph'])
    spectra = []
    for item in config['spectra']:
        spectra.append(Spectrum(item))
    draw_tripartite()

#### 出力画像(Tripartite.png)

 

### 解説

以前のプログラムTipsで紹介した事柄は割愛したいと思います。
```if __name__ == "__main__"```の意味やクラスの呼び出し方などは過去の記事をご覧ください。

【プログラムTIPS】pythonでCSVファイルを固定書式に変換する

#### 速度応答スペクトルデータ

速度応答スペクトルは固有周期と減衰定数をパラメータとして、ある地震動に対する１質点系モデルの最大応答速度をプロットして線で繋ぐことで得られます。実際に応答スペクトルを得るためには、Nigam・Jennings法などを用いて応答計算を行い最大値をプロットするも良し、定式化された計算式を用いてプロットするも良し、既存プログラムを利用するも良し、だと思います。私が学生だったときは最初に挙げた方法でスペクトルを計算していましたが、今回は既存プログラムを利用しました。RESP-F3Tであれば、以下のような簡単な記述でスペクトルを取得することができます。

InputUnit N cm
<
< 波形
<
Wave EL-CENTRO_NS 2688 0.0200 Format "(7F10.1)" File "C:\Program Files (x86)\KKE\RESP-D\SampleWave2006s.dat" Skip 1 Magnify 1
Wave TAFT_EW 2714 0.0200 Format "(7F10.1)" File "C:\Program Files (x86)\KKE\RESP-D\SampleWave2006s.dat" Skip 1429 Magnify 1
Wave HACHINOHE_NS 2550 0.0200 Format "(10F7.2)" File "C:\Program Files (x86)\KKE\RESP-D\SampleWave2006s.dat" Skip 2488 Magnify 1
<
< 解析コマンド
<
<
Analysis  ResponseSpectrum
    Wave EL-CENTRO_NS TAFT_EW HACHINOHE_NS
    Period Begin 0.02 End 10.0 LogScale
    OutputResponse Velocity

#### 設定ファイル

json形式のデータです。データ形式の詳細は「json」でググると多くの記事がありますのでここでは割愛しますが、人にもある程度読みやすい形になっているかと思います。大きく分けて「graph」、「spectra」という２つの項目を用意して、それぞれで、描画する図自体に関する設定と、スペクトルに関する設定を行えるようにしています。例えば、「graph」の中にある「output_file_name」では出力画像ファイル名を、「spectra」の中にある「path」ではスペクトルデータのファイルパスを設定しています。

{
  "graph":
    {
      "output_file_name":"Tripartite.png",
      "dpi":600,
      "graph_width":5.0,
      "graph_height":6.5,
      "min_period":0.01,
      "max_period":10.0,
      "min_vel":0.1,
      "max_vel":1000.0,
      "x_label":"T (S)",
      "y_label":"Sv (cm/s)",
      "acc_label":"Sa (cm/s/s)",
      "disp_label":"Sd (cm)",
      "eigen_period":[1.5, 0.7, 0.2]
    },
  "spectra":
    [
      {
        "name":"EL-CENTRO_NS",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":3,
        "color":"deepskyblue"
      },
      {
        "name":"TAFT_EW",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":4,
        "color":"black"
      },
      {
        "name":"HACHINOHE_NS",
        "path":"./spectrum_data/response_spectrum.last.csv",
        "skip_row":5,
        "period_col":2,
        "spectrum_col":5,
        "color":"red"
      }
    ]
}

#### pythonコード

ようやく、、、プログラムの核であるpythonのソースコードについて触れます。
まず、プログラムを実行すると以下の部分が呼ばれます。

if __name__ == '__main__':
    config = read_config('./config.json')
    graph = GraphConfig(config['graph'])
    spectra = []
    for item in config['spectra']:
        spectra.append(Spectrum(item))
    draw_tripartite()



ここでは、

１．設定ファイル読み込み

２．読み込んだ設定ファイルに従って分かりやすく簡単に各種設定や描画処理が行えるように、オブジェクト生成（描画のための前処理）

３．トリパタイト描画処理実行命令

という流れになっています。「graph」はGraphConfigクラスのオブジェクト、「spectra」はSpectrumクラスのオブジェクトの配列であり、設定ファイルの内容とリンクするように生成されています。オブジェクト指向プログラミングでは、こうしたオブジェクトを生成することでプログラムを構成しやすくしています。一方で、この規模のプログラムであれば、わざわざクラスを用意せずにjsonを読み込んで作成した辞書オブジェクトをそのまま使用する方法も考えられます。


続いて以下が呼ばれます。


def draw_tripartite():
    fig, ax = graph_init()
    draw_grids(fig, ax)
    draw_vel_spectrum(ax)
    draw_period_line(ax)
    graph_format(ax)
    plt.savefig(graph.FileName, dpi=graph.Dpi)
    plt.show()



ここでは、

１．グラフの初期設定

２．グリッド描画（加速度と変位のグリッド）

３．速度応答スペクトル描画

４．固有周期の線描画

５．グラフ調整

６．グラフ保存

７．グラフ表示

という流れです。以降では、トリパタイトを描くうえで主要な処理である２、３の項目について中身を見ていきます。


* グリッド描画（加速度と変位のグリッド）

該当部分を再掲します。


def get_grid_lines_of_acc_and_disp():
    x1 = graph.MinPeriod
    x2 = graph.MaxPeriod
    # acc
    acc_list = np.concatenate([np.linspace(int(10**i), int(10**(i+1)), 10) for i in range(4)])
    acc_lines = {}
    for acc in acc_list:
        # Sv = Sa / (2π/T)
        y1 = acc / (2.0 * np.pi / x1)
        y2 = acc / (2.0 * np.pi / x2)
        acc_lines[acc] = [x1, x2, y1, y2]
    # disp
    disp_list = np.concatenate([np.linspace(int(0.01*10**i), int(0.01*10**(i+1)), 10) for i in range(4)])
    disp_lines = {}
    for disp in disp_list:
        # Sv = Sd * (2π/T)
        y1 = disp * (2.0 * np.pi / x1)
        y2 = disp * (2.0 * np.pi / x2)
        disp_lines[disp] = [x1, x2, y1, y2]
    return acc_lines, disp_lines

def draw_grid(grid_type, fig, ax, label, lines, text_pos, angle, vertalalign):
    for k, v in lines.items():
        x1, x2, y1, y2 = tuple(v)
        ax.plot([x1, x2], [y1, y2], color='darkgray', linewidth='0.5')
        text = '{:.0f}'.format(k)
        if (grid_type == 'acc'):
            if (np.log10(k)).is_integer():
                ax.text(x1, y1, text, rotation=angle, verticalalignment=vertalalign)
        elif (grid_type == 'disp'):
            if k &amp;amp;amp;amp;amp;gt;= 1.0 and (np.log10(k)).is_integer():
                ax.text(x2, y2, text, rotation=angle, verticalalignment=vertalalign)
    fig.text(text_pos[0], text_pos[1], label, rotation=angle)

def draw_grids(fig, ax):
    acc_lines, disp_lines = get_grid_lines_of_acc_and_disp()
    draw_grid('acc', fig, ax, graph.AccLabel, acc_lines, [0.3, 0.65], 45, 'bottom')
    draw_grid('disp', fig, ax, graph.DispLabel, disp_lines, [0.7, 0.75], -45, 'top')



加速度グリッドと変位グリッドはまずそれぞれどのような刻み（目盛り）にしたいかを考えて、その刻み毎の加速度あるいは変位を速度に変換することで速度応答スペクトル図上に加速度と変位の線を描く、という考え方になるかと思います。５行目では加速度グリッド用の目盛りリストを作成しています。一行に収めたため少し分かりづらいですが、ログ表示に対応した目盛りとするため、1～10、10～100、100～1000、1000～10000のそれぞれの区間を10分割した配列を作成しています。そして６行目からSv=Sa/ωの関係に基づいて加速度を速度に変換して速度応答スペクトル図上に線を描画するためのデータ（始点と終点の値がセットとなったデータ）にしています。変位グリッドも同様です。変位の場合はSv=Sd*ωで変換します。２２行目から実際にグリッドとテキスト（目盛りの値など）を描画する処理を記述しています。ちなみに、２行目にあるgraph.MinPeriodは設定ファイルで定義した最小固有周期（グラフ横軸の最小値）を意味しますが、最初の描画前処理で「graph」というオブジェクトを作成していたため、このように簡単に設定データを参照できるようになっています。

※ここでは、numpyというライブラリも使用しています（コード中の「np」）。


* 速度応答スペクトル描画

こちらも該当部分のみを再掲します。


def read_vel_resp_spectrum(path:str, skip_row:int, period_col:int, spectrum_col:int):
    df = pd.read_csv(path, header=None, skiprows=skip_row)
    return df[period_col].values, df[spectrum_col].values

def draw_vel_spectrum(ax):
    for spec in spectra:
        x, y = read_vel_resp_spectrum(spec.Path, spec.SkipRow, spec.PeriodCol, spec.SpectrumCol)
        ax.plot(x, y, label=spec.Name, color=spec.Color)



変換も何も無く、ただ描画するだけであるためグリッド描画よりシンプルです。スペクトルデータを設定ファイルで定義した仕様（ファイルパスや読み飛ばし行数、データ列など）に基づいて読み込み（７行目）、読み込んだデータをmatplotlibという描画ライブラリのplotという関数に渡すだけです（８行目）。ここで同時に、設定ファイルで定義した凡例の文字列や色を指定しています。

※ここでは、pandasというライブラリも使用しています（コード中の「pd」）。




 


## まとめ


今回は各種細かな設定ができるようなものを作成したため、少し複雑になった部分もありますが、機能を限定すればもっと簡潔に短いコード量で作成することも可能です。ライブラリの仕様や細かな処理についてはここでは触れられませんでしたが、本記事がpythonによる簡単なプログラム作成のための参考になれば幸いに思います。


完成したプログラム全体はGitHubに公開しています。下記よりご覧いただけます。

※本プログラムによるいかなるトラブル、損失、損害について弊社は一切責任を負いません。本記事に掲載している図やプログラムコードはあくまでも一例です。

トリパタイト（三軸図）描画プログラム


関連する製品 RESP-F3T


RESP-F3T

3次元フレーム汎用解析プログラム



詳細はこちらから


採用情報


構造計画研究所 RESPチームでは、いっしょに働いていただけるエンジニアを募集しています。

構造設計・構造解析だけでなく、プログラミング技術を活かして新しいものを生み出したいと思っている方、ぜひご応募ください。


採用HPはこちら→https://www.kke.co.jp/recruit/

返信を残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です