iRICの計算結果をPythonでいじってみる：実河川の例

はじめに

河川に関する数値計算モデルプロジェクトiRICでは多くの数値計算モデルが計算条件設定の前処理と可視化の後処理を備えたGUIとともに使用可能です．

いろいろな機能があり，計算結果はサクサク確認できますが，細かな解析をしたかったり，こだわりの可視化をしたかったりとすると，やはり備え付けのGUIだけでは足りない部分も出てきます．

ここではiRICのモデルの一つであるNays2DHの計算結果をcsvで出力し，Pythonでいじくってみたというのを紹介していきます．

計算結果のcsvエクスポート

iRICではいくつかの形式で計算結果を外部ファイルに出力できます．

まず，上部のメニューバーから　ファイル→エクスポート→計算結果　と進みます．

出てきた画面で下記のように設定します．

形式はcsvファイルとします．ほかにもvtkなどでも出力できます．

時刻や領域は間引きができますが，ここでは行わないことにします．

その結果，プレフィックスで設定した名前の後ろに数字が付く（例：Result_1.csv）感じで連番ファイルが出力されます．

エクセルで開くと以下のようになっていて，一行目が時間，二行目が格子数，その下が計算結果になります．

Pythonによる読み込み

csvファイルを読み込んで解析したり，可視化したりしてみましょう．ここではPythonを使っていきます．

csvファイルを読み込んで，時間，x, y座標，河床高さ，水深を返す関数を作ってみます．

ifpが読み込むcsvファイルになります．

座標や河床変動，水深に相当する列のデータを指定します．どのデータがどのような順番で出力されているかはそれぞれのモデルによるので注意してrow[2]などの番号を設定します．この列番号は0から始めっていることに注意です．

import numpy as np 
import matplotlib.pyplot as plt 

def readResultfile(ifp):
    
    with open(ifp, 'r') as ifile:
            
        data    = csv.reader(ifile, delimiter=",")
        header  = str(next(data))
        number  = next(data)
        index   = next(data)
        
        wwt = header.replace('\n', '').split()
        wt  = str(wwt[-1])
        tt  = float(wt[0:-2])   # 時間（秒）
        
        nx = int(number[0])-1   # x方向格子数
        ny = int(number[1])-1   # y方向格子数
                
        x = []
        y = []
        z = []
        h = []
            
        X   = np.zeros((ny+1,nx+1))
        Y   = np.zeros((ny+1,nx+1))
        Z   = np.zeros((ny+1,nx+1))
        H   = np.zeros((ny+1,nx+1))
                            
            # CSVファイルからデータを読み込み
            
        for row in data:
            x.append(float(row[2]))   # x座標
            y.append(float(row[3]))   # y座標
            z.append(float(row[5]))   # 河床高さ
            h.append(float(row[4]))   # 水深
        
            # 1次元配列を二次元配列に変換
        
        ij = 0
        for j in range(ny+1):
            for i in range(nx+1):
                X[j][i] = x[ij]
                Y[j][i] = y[ij]
                Z[j][i] = z[ij]
                H[j][i] = h[ij]
                    
                ij = ij+1
                                            
    ifile.close()
    
    return tt, X, Y, Z, H

フォルダの指定

読み込みや出力先フォルダ名を設定します．

iRICからcsvファイルを出力したフォルダや，図を出力したいフォルダ名を指定します．

import csv
import os
import shutil

root = "E:/work/test/"                      # 親フォルダ
case = "1-550"                               # ケース名
csv_path = root+"csv/"+case+"/"               # csvが入っているフォルダ
fig_path = root+"fig/"+case+"/"               # ↑に対応する画像ファイルが入るフォルダ

if os.path.exists(fig_path):
    shutil.rmtree(fig_path)

os.mkdir(fig_path)

import matplotlib.cm as cm
import math
import re
import copy


wip = os.listdir(csv_path)

wwip = []

    #  番号順にリストを並び替える

for ip in wip:
    mmm = re.search("[0-9]+", ip)
    www = (mmm.group(), ip)
    wwip.append(www)

aaa = sorted(wwip, key=lambda wwip: int(wwip[0]))

inputfiles = []

for aa in aaa:
    inputfiles.append(aa[1])

地形と堤防位置の確認

ここでは北海道の音更川の計算結果を例にしていますが，この地形と堤防位置を可視化してみます．

まず，図のおさまりをよくするために，適当に回転させて，軸の数字が大きくなり過ぎないようにオフセットを取っておきます．

filename = csv_path+inputfiles[0]

  #  回転させる
    
d = -65.

d_rad = math.radians(d)

time, wX, wY, Z, H = readResultfile(filename)
    
X = wX * np.cos(d_rad)-wY*np.sin(d_rad)
Y = wX * np.sin(d_rad)+wY*np.cos(d_rad)

   # オフセットをとる

x = np.min(X)
y = np.min(Y)

X = X - x
Y = Y - y

堤防の部分では地形の勾配が大きくなっているはずなので，横断方向に曲率をとるラプラシアンフィルタを利用することを考えます．

ただ，これだけだとうまくいかなかったため，地形の高さも利用します．ところが，ある程度の区間だと勾配分で上流と下流で河床高さは大きく違うので，横断面内の平均河床高さを差し引くことで，比高のようなものに変換してみます．

# フィルタ計算する関数　→　堤防位置の検知に利用してみる
def filter2D(src, kernel):
    # カーネルサイズ
    m, n = kernel.shape
    # 畳み込み演算をしない領域の幅
    ignore = int((m - 1) / 2)
    height, width = src.shape
    
    img_edge = np.zeros((height, width))
    
    for y in range(ignore, height - ignore):
        for x in range(ignore, width - ignore):
            # 畳み込み演算
            img_edge[y][x] = np.sum(src[y - ignore: y + ignore + 1, x - ignore: x + ignore + 1] * kernel)
    
    return img_edge

ny, nx = np.shape(Z)

nor_z = copy.copy(Z)

ave_z = np.mean(Z,axis=0)

  # 横断面内の平均河床高さを引いて勾配分を取り除いた形の地形を作る

for i in range(nx):
    for j in range(ny):
        nor_z[j][i] = Z[j][i]-ave_z[i]

    # 上記作成した地形に対して横断方向にラプラシアンフィルタ－をかけて堤防位置特定に利用

kernel_y = np.array([
    [ 0,  1,  0],
    [ 0, -2,  0],
    [ 0,  1,  0],
])

fil_z = filter2D(nor_z, kernel_y)

堤防位置を横断方向の河床高さ曲率と比高を基準としてトライアンドエラーで抽出していきます．いろいろ手作業過ぎるコードですが．．

xl = []
yl = []
xr = []
yr = []
il = []
jl = []
ir = []
jr = []

nyc = int(0.5*ny)
zt = 0.7
dz = -0.4

  #  横断面平均河床高さとラプラシアンの組み合わせによる堤防位置の特定

for i in range(nx):
    nl = 0
    nr = 0
    for j in reversed(range(0,nyc)):
        if nor_z[j][i]>zt and fil_z[j][i]<dz and nl==0:
            xl.append(X[j][i])
            yl.append(Y[j][i])
            il.append(i)
            jl.append(j)
            nl = nl+1
            
    for j in range(nyc,ny):
        if nor_z[j][i]>zt and fil_z[j][i]<dz and nr==0:
            xr.append(X[j][i])
            yr.append(Y[j][i])
            ir.append(i)
            jr.append(j)
            nr = nr+1

比高差と左右岸の堤防位置をmatplotlibで確認してみましょう．

plt.figure(figsize=(10,2.5),tight_layout=True)
plt.axes().set_aspect('equal')
plt.xlabel('X(m)')
plt.ylabel('Y(m)')
plt.subplots_adjust(bottom=0.5)

plt.contourf(X, Y, nor_z)
plt.colorbar(orientation='horizontal', pad=0.3)

plt.plot(X[jl,il],Y[jl,il],'o',color='r',markersize=0.3)
plt.plot(xr,yr,'o',color='b',markersize=0.3)
plt.show()
plt.close()

勾配分の河床高さを差し引いたことで形が見やすくなり，堤防位置も特定できているようです．

計算結果の可視化

csvの連番ファイルを読んで結果をjpgファイルで出力してみます．ここでは水深を出力していきます．上記で特定した堤防位置も同時に書き込むことで，堤防まで流路が移動しているかどうか，確認することができます．

for inputfile in inputfiles[::1]:
    filename = csv_path+inputfile
    figname = fig_path+inputfile+".jpg"
    
    print(filename)
      
    tt, wX, wY, Z, H = readResultfile(filename)
    
    X = wX * np.cos(d_rad)-wY*np.sin(d_rad)
    Y = wX * np.sin(d_rad)+wY*np.cos(d_rad)
    
    x = np.min(X)
    y = np.min(Y)
    
    X = X - x
    Y = Y - y
    
            
    plt.figure(figsize=(10,2.5),tight_layout=True)
    plt.axes().set_aspect('equal')
    plt.xlabel('X(m)')
    plt.ylabel('Y(m)')
    plt.subplots_adjust(bottom=0.2)
    
    levels = np.linspace(0.1,5,51)
    labels = np.linspace(0.1,5,6)

    plt.title("Time= {0:.0f} min".format(tt/60)) 

    plt.contourf(X, Y, H, levels, cmap=cm.Blues, extend="max")
    plt.colorbar(ticks=labels, label='Depth (m)', orientation='horizontal', fraction=0.08, pad=0.3)
    plt.plot(xl,yl,'o',color='k',markersize=0.2)
    plt.plot(xr,yr,'o',color='k',markersize=0.2)
    plt.savefig(figname, dpi=600)
    plt.close()

計算終了時の結果は以下のようになっていて，一部の流路が堤防付近まで及んでいることがわかります．

連番動画を結合して動画にする

opencvを使ってmp4の動画ファイルを作ってみましょう．先ほど出力したjpgファイルのリストを作成して動画を作っていきます．コーデックなどはこの記事を参考にしました．

import cv2
import glob

from natsort import natsorted

img_array = []

filenames = natsorted(glob.glob(fig_path+"*.jpg"))

for filename in filenames:
    img = cv2.imread(filename)
    height, width, layers = img.shape
    size = (width, height)
    img_array.append(img)
    
    
name = root+case+'.mp4'
codec = cv2.VideoWriter_fourcc(*'X264')
framerate = 10.
out = cv2.VideoWriter(name, codec, framerate, size)

for i in range(len(img_array)):
    out.write(img_array[i])
out.release()

これで問題なくできると思います．コーデックの調整が必要な場合もあります．