はじめに

　この記事がちゃんとした(?)投稿としては初めてになります。至らない点がたくさんあるかと思いますが、温かい目で見ていただけると幸いです

概要

　内閣府での南海トラフや日本海溝・千島海溝沿いの巨大地震モデル検討会において想定されている津波データを用いて、津波が陸地を遡上する様子をアニメーションで描画してみます。

1. データを用意する

　描画する津波のデータには、内閣府での巨大地震モデル検討会で作成されたデータを用います。このデータはG空間情報センター[1]というサイトで公開されています。

　ユーザー登録をしないとデータをダウンロードできないので、ユーザー登録を行いましょう。

提供されているデータのうち、津波断層モデルに関わるものは

(1) 初期水位データ
(2) 地殻変動量データ
(3) 津波断層パラメータ
(4) 陸域における津波浸水深データ(最大包絡地)
(5) 海岸における津波高・津波到達時間データ
(6) 波形データ
(7) 液状化データ

の7つありますが、描画には津波の浸水深と到達時間のデータがあればよいので、

(4) 陸域における津波浸水深データ(最大包絡値)
(5) 海岸における津波高・津波到達時間データ

をダウンロードします。データはcsvファイルの形式でダウンロードします。

　データの中身は巨大地震モデルによって異なります。
　ここでは日本海溝モデルでのデータを使います。上記(4)のデータをダウンロードするとzipファイルに5つのcsvファイルが入っていました。

　ファイルの名前の最後に「第〇系」とあります。これは日本平面直交座標の系を示しています。
　都道府県と系の対応は以下のサイト[2]を参考にするとよいでしょう。

　ここでは、岩手県宮古市の地域での津波を描画することとします。岩手県は第10系に所属していますので、第10系のデータを使用します。
　開くと以下のようなデータが表示されました。

　データには緯度経度での位置情報と最大浸水深と最大浸水深となる時間(秒)、そして30cm、1m、3mの浸水深となる津波の到達時間(秒)情報が含まれていました。

　このcsvファイルを開いてみると、以下のようなメッセージが出ることがあります。

　Microsoft Excelでは扱うことができる行の最大行はおよそ100万行であり、それを超えるデータを扱う場合はデータのフィルタリングを行う必要があります。
　フィルタリングを行う方法は以下のサイト[3]を参考にするとよいでしょう。

　データには必ず緯度経度情報が含まれているので、上サイトでの方法を用いて自分が描画したい範囲の緯度経度のデータのみを取り出します(4章で、自分が描画したい範囲をどう得るのかについて記載します)。
　これでまず必要な津波データについては用意ができました。

2. 緯度経度を直交座標系に変換する

　津波を描画する際に描画先として指定する座標は、直交座標系となります。先ほど読み込んだデータの位置情報は緯度経度で表されていたので、緯度経度系⇒日本平面直角座標系への座標変換とOpenGLでの座標系における原点への平行移動を行う必要があります。
　
　そのためのコードを以下に示します。

temp_pre_lat_list = [] #変更前経度
temp_pre_lon_list = [] #変更前緯度
pre_lat_list = []
pre_lon_list = []

csv_transpose_list = [] #csv転置用

ox = const.EPSG_OX   #北西端x座標
oy = const.EPSG_OY   #北西端y座標

EPSG4326 = pyproj.Proj("+init=EPSG:4326")  #WGS84(GPSで得られる位置)
EPSG6678 = pyproj.Proj("+init=EPSG:6678")  #JGD2011

csv_file = open("座標変換前の津波データ.csv")

for row in csv.reader(csv_file):
    temp_pre_lon_list.append(row[0])
    temp_pre_lat_list.append(row[1])
    
del temp_pre_lon_list[0] #ヘッダーを読み飛ばす
del temp_pre_lat_list[0]

for item in temp_pre_lon_list:
    pre_lon_list.append(float(item))

for item in temp_pre_lat_list:
    pre_lat_list.append(float(item))
    
lon_list = tuple(pre_lon_list)
lat_list = tuple(pre_lat_list)

lon_lat = list(pyproj.transform(EPSG4326, EPSG6678, lon_list, lat_list)) #この時点ではcsvにデータが行で格納されてしまう

#平行移動込み
for i in range(len(lon_list)):
    csv_transpose_list.append((round(lon_lat[0][i]-ox), round(lon_lat[1][i]-oy))) #ここで列データに格納できるようにする

with open("座標変換後の津波データ.csv", mode = "w", newline = '') as f:
    writer = csv.writer(f)
    writer.writerows(csv_transpose_list)

3. データを読み込む

　データはcsvファイルなので、csvモジュールを用いてデータを読み込みます。
　先ほど示したようにデータは列ごとに種類が分かれているので、列ごとにデータをリストに格納していきます。

　以下にコードを示します。

WaveReader.py

class WaveReader(object):
    def __init__(self, filename):
        #津波データの読み込み
        self.water_points = []
        self.wave_levels = [0.3, 1.0, 3.0]
        self.x_list = []
        self.y_list = []
        self.tstep_30cm_list = []
        self.tstep_1m_list = []
        self.tstep_3m_list = []
        
        print("Loading Wave data...")

        csv_file = open(filename) #座標変換後のデータを使う
        reader = csv.reader(csv_file)
        header = next(reader)  # ヘッダーを読み飛ばす

        for row in reader:
            tstep_x = float(row[0]) * X_SCALE_ADJUSTMENT_RATIO #X_SCALE_ADJUSTMENT_RATIO:テクスチャ画像に合うように変換
            tstep_y = float(row[1]) * Y_SCALE_ADJUSTMENT_RATIO
            self.x_list.append(tstep_x)
            self.y_list.append(tstep_y)  
            # 秒データを計算ステップに変換
            tstep_30cm = float(row[2]) / DT_SEC #DT_SEC:津波描画の時間ステップ(今回は10秒=1ステップと設定し、DT_SEC=10とする)
            tstep_1m = float(row[3]) / DT_SEC
            tstep_3m = float(row[4]) / DT_SEC
            self.tstep_30cm_list.append(tstep_30cm)
            self.tstep_1m_list.append(tstep_1m)
            self.tstep_3m_list.append(tstep_3m)

        print("Completed Loading Wave data.")

　ここではcsvにあるデータを列ごとに格納しているのですが、各データをちょこっといじります。
　まず、緯度経度データについては津波を描画する際にそのウィンドウサイズの座標比をそろえる必要があるので、その比をかけます。
　また、浸水時刻データは秒数で表されていますが、計算コスト削減のため、ここでは時間ステップで表現することとします。

4. 津波の描画

　津波の描画を行うためのデータの準備はできました。ここからは津波の描画について説明します。この記事では、日本海溝巨大地震による津波の襲来が懸念されている岩手県宮古市河南地区において想定されている津波を描画します。

4-1 背景のテクスチャ

　津波を描画する際に、背景としてその地域の航空画像のテクスチャを貼り付けます。
　航空画像の取得には様々な方法がありますが、一番簡単な方法としては「Google Earth Pro」を用いる方法がよいでしょう。Google Earth Pro には目印を追加することで、下図のように地図上に画鋲が表示されます。この画鋲を描画範囲の緯度経度と合わせて刺して、スクリーンショットを撮影し、画鋲に合わせて画像をトリミングすれば背景に使う航空画像が得られます。

　ちなみに描画範囲の決定にはOpenStreetMap[4]を使用しました。今回はOpenStreetMapのデータを直接使用する訳ではありませんが、データを取得する「エクスポート」の機能では、表示している地図の東西南北の経緯がわかるので便利かと思います。

　このようにして津波を表示させる範囲を指定します。この緯度経度の範囲をメモしておいて、1章においてExcelデータのフィルタリングを行いましょう。
　
　ここからはテクスチャを読み込む方法についてです。まず、ライブラリをインポートします。他に使用するライブラリも合わせてインポートしてしまいましょう。

main.py

from pyglet import image
from pyglet.gl import *
from pyglet.window import *
from WaveReader import WaveReader
from time import sleep
import const

　pyglet[5]はゲームやマルチメディアアプリケーションの作成に用いられるライブラリで、画面の設定やテクスチャの貼り付けも行います。constは定数を格納した自作ライブラリです。
　テクスチャの話に戻ります。出力画面の設定とそこに貼り付けるテクスチャの設定は以下のようにします。

main.py

#出力画面の設定
window = pyglet.window.Window(const.WINDOW_WIDTH_DEFAULT, const.WINDOW_HEIGHT_DEFAULT)
#const.WINDOW_WIDTH_DEFAULT=テクスチャの横方向ピクセル数
#const.WINDOW_HEIGHT_DEFAULT=テクスチャの縦方向ピクセル数

#テクスチャ読み込み
pic = image.load('aerial_photo/kanan.png')
sprite = pyglet.sprite.Sprite(img = pic)

def setup():
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glClearColor(1, 1, 1, 1)
    glLoadIdentity()

　また、出力画面に時間経過を表示させるといいでしょう。

main.py

#経過時間表示
label_t = pyglet.text.Label("",
                            font_name='Meiryo',
                            font_size=30,
                            #color=(0, 0, 0, 255),
                            bold=True,
                            x=const.WINDOW_WIDTH_DEFAULT * 0.03, y=const.WINDOW_HEIGHT_DEFAULT * 0.95,
                            anchor_x='left', anchor_y='center')
@window.event 
def on_draw():
    global Mesh_list
    window.clear()
    sprite.draw()

    #1clock=10sで表示
    now_seconds = clock * 10
    h = '{0:02d}'.format(int(now_seconds / 3600))
    m = '{0:02d}'.format(int(int(now_seconds / 60) % 60))
    s = '{0:02d}'.format(int(now_seconds % 60))
    
    hms = h + "h " + m + "m " + s + "s"
    
    label_t.text = hms
    label_t.draw()

    for mesh in Mesh_list: #オブジェクトの描画
        mesh.draw_water()

　このように設定すると以下のように経過時間が表示されました。一番最後の2行については後ほど説明します。

4-2 津波メッシュの作成

　まず、3章で読み込んだ津波メッシュデータからメッシュ1つずつ取り出し、draw_meshのインスタンスを生成しリストに格納します。draw_meshではメッシュサイズの設定や描画時の色の設定を行います。ここで、格納するデータは

・データ地点情報(x, y)
・30cm,1m,3m津波の到達時間(時間ステップ)

となります。津波の到達時間に関しては、もともと秒でのデータでしたがWaveReaderにて時間ステップに変換しています。先ほども示したように、時間ステップは1ステップあたり10秒としています。

main.py

class instance_draw_mesh(object):
    def __init__(self):
        global Mesh_list
        self.wavereader = WaveReader("wave_csv/宮古市河南津波データ.csv")
        self.x_list = self.wavereader.x_list
        self.y_list = self.wavereader.y_list
        self.tstep_30cm_list = self.wavereader.tstep_30cm_list
        self.tstep_1m_list = self.wavereader.tstep_1m_list
        self.tstep_3m_list = self.wavereader.tstep_3m_list

        for i in range(len(self.x_list)):
            self.x = self.x_list[i]
            self.y = self.y_list[i]
            self.t_30cm = self.tstep_30cm_list[i]
            self.t_1m = self.tstep_1m_list[i]
            self.t_3m = self.tstep_3m_list[i]
            Mesh_list.append(draw_mesh(self.x, self.y, self.t_30cm, self.t_1m, self.t_3m))

class draw_mesh(object):
    def __init__(self, x, y, t_30cm, t_1m, t_3m):
        self.mesh_size = 10.0
        self.half_mesh = self.mesh_size / 2.0
        self.color_list = [[0.4, 0.4, 1.0], [0.1, 0.1, 1.0], [0.0, 0.0, 1.0]]
        self.x = x
        self.y = y
        self.t_30cm = t_30cm
        self.t_1m = t_1m
        self.t_3m = t_3m

4-3 津波メッシュの描画

　津波を描画するメソッドをdraw_meshに書き込みます。ここでは以下のような2つのメソッドを考えました。

main.py

#class draw_mesh内のメソッド

#津波の描画
    def draw_water(self):
        global Mesh_list, clock
        if clock > self.t_30cm:
            self.color = self.color_list[0]
            if clock > self.t_1m:
                self.color = self.color_list[1]
                if clock > self.t_3m:
                    self.color = self.color_list[2]
                    self.draw_square(self.color, self.x, self.y)
                else:
                   self.draw_square(self.color, self.x, self.y)
            else:
                self.draw_square(self.color, self.x, self.y)
        else:
            return

    def draw_square(self, color, x, y):
        glColor3f(color[0], color[1], color[2])
        glBegin(GL_QUADS) 
        glVertex2i(int(x - self.half_mesh), int(y + self.half_mesh))
        glVertex2i(int(x + self.half_mesh), int(y + self.half_mesh))
        glVertex2i(int(x - self.half_mesh), int(y - self.half_mesh))
        glVertex2i(int(x + self.half_mesh), int(y - self.half_mesh))
        glEnd()

　draw_waterでは、clock(後述)と各高さの津波の到達時間と比較し、浸水判定を行います。
　draw_squareでは、メッシュを描画する際の四角形を定義しています。

　これらのメソッドを時間経過に応じて繰り返す必要があります。pygletでは同じ処理を繰り返す際にpyglet.clock.schedule(<繰り返す処理>)というメソッドを用います。ここでの<繰り返す処理>としてupdateという関数を定義します。

main.py

def update(dt):
    global Mesh_list, activeFlag, clock
    if activeFlag:
        for mesh in Mesh_list:
            mesh.draw_water()
        clock += 1
    sleep(0.2)
    return True

　ここでclockを0から1つずつ増加させていき、draw_waterにて浸水判定をします。

5 結果

　ここまでの手順をまとめたmainコードを以下に示す。

main.py

from pyglet import image
from pyglet.gl import *
from pyglet.window import *
from WaveReader import WaveReader
from time import sleep
import const

window = pyglet.window.Window(const.WINDOW_WIDTH_DEFAULT, const.WINDOW_HEIGHT_DEFAULT)

Mesh_list = []
activeFlag = True
clock = 0

#テクスチャ読み込み
pic = image.load('aerial_photo/kanan.png')
sprite = pyglet.sprite.Sprite(img = pic)

#経過時間表示
label_t = pyglet.text.Label("",
                            font_name='Meiryo',
                            font_size=30,
                            #color=(0, 0, 0, 255),
                            bold=True,
                            x=const.WINDOW_WIDTH_DEFAULT * 0.03, y=const.WINDOW_HEIGHT_DEFAULT * 0.95,
                            anchor_x='left', anchor_y='center')

def setup(): #背景の色の設定
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glClearColor(1, 1, 1, 1)
    glLoadIdentity()

@window.event
def on_draw():
    global Mesh_list
    window.clear()
    sprite.draw()

    #1clock=10sで表示
    now_seconds = clock * 10
    h = '{0:02d}'.format(int(now_seconds / 3600))
    m = '{0:02d}'.format(int(int(now_seconds / 60) % 60))
    s = '{0:02d}'.format(int(now_seconds % 60))

    hms = h + "h " + m + "m " + s + "s"

    label_t.text = hms
    label_t.draw()

    for mesh in Mesh_list: #オブジェクトの描画
        mesh.draw_water() 

def update(dt):
    global Mesh_list, activeFlag, clock
    if activeFlag:
        for mesh in Mesh_list:
            mesh.draw_water()
        clock += 1
    sleep(0.2)
    return True

class instance_draw_mesh(object):
    def __init__(self):
        global Mesh_list
        self.wavereader = WaveReader("wave_csv/宮古市河南津波データ.csv")
        self.x_list = self.wavereader.x_list
        self.y_list = self.wavereader.y_list
        self.tstep_30cm_list = self.wavereader.tstep_30cm_list
        self.tstep_1m_list = self.wavereader.tstep_1m_list
        self.tstep_3m_list = self.wavereader.tstep_3m_list

        for i in range(len(self.x_list)):
            self.x = self.x_list[i]
            self.y = self.y_list[i]
            self.t_30cm = self.tstep_30cm_list[i]
            self.t_1m = self.tstep_1m_list[i]
            self.t_3m = self.tstep_3m_list[i]
            Mesh_list.append(draw_mesh(self.x, self.y, self.t_30cm, self.t_1m, self.t_3m))

class draw_mesh(object):
    def __init__(self, x, y, t_30cm, t_1m, t_3m):
        self.mesh_size = 10.0
        self.half_mesh = self.mesh_size / 2.0
        self.color_list = [[0.4, 0.4, 1.0], [0.1, 0.1, 1.0], [0.0, 0.0, 1.0]]
        self.x = x
        self.y = y
        self.t_30cm = t_30cm
        self.t_1m = t_1m
        self.t_3m = t_3m

    #津波の描画
    def draw_water(self):
        global Mesh_list, clock
        if clock > self.t_30cm:
            self.color = self.color_list[0]
            if clock > self.t_1m:
                self.color = self.color_list[1]
                if clock > self.t_3m:
                    self.color = self.color_list[2]
                    self.draw_square(self.color, self.x, self.y)
                else:
                   self.draw_square(self.color, self.x, self.y)
            else:
                self.draw_square(self.color, self.x, self.y)
        else:
            return

    def draw_square(self, color, x, y):
        glColor3f(color[0], color[1], color[2])
        glBegin(GL_QUADS) 
        glVertex2i(int(x - self.half_mesh), int(y + self.half_mesh))
        glVertex2i(int(x + self.half_mesh), int(y + self.half_mesh))
        glVertex2i(int(x - self.half_mesh), int(y - self.half_mesh))
        glVertex2i(int(x + self.half_mesh), int(y - self.half_mesh))
        glEnd()

def main():
    instance_draw_mesh() #インスタンスの生成
    pyglet.clock.schedule(update) #指定した時間(dt)のたびに関数(update)を呼び出す
    pyglet.app.run() #実行

if __name__=='__main__':
    main()

　これを実行すると、以下のようになりました。

途中、津波の描画が止まる部分があると思いますが、これは当地域の陸閘の影響があるようです。
また、今回出力したgifでは津波の描画が進むにつれて、描画が遅くなってしまう問題があります(スナップショットを残しておき、パラパラ漫画のようにつなげていけば一定速度で描画できているようには見えますが…)。

6 まとめ

　今回は、内閣府巨大地震モデル検討会の津波断層モデルを用いて津波の遡上の様子を描画しました。このデータはとてもサイズが大きく、メッシュデータに落とし込むのにも使い勝手があまりよろしくなかったですが、描画範囲をしぼれば割と苦労なく描画ができたと思います。

参考文献

[1] G空間情報センター https://front.geospatial.jp/
[2] 日本平面直角座標の系 https://www.sinfonica.or.jp/faq/gis/minf_hzahyo.html
[3] Excel最大行（100万行）の壁を軽々超える方法 https://atmarkit.itmedia.co.jp/ait/articles/2109/16/news024.html
[4] OpenStreetMap https://www.openstreetmap.org/
[5] pyglet https://pypi.org/project/pyglet/