趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

pythonの関数で複数の戻り値を設けた場合の型による速度比較

勉強会で話が出たので整理しました。



速度比較

 複数の戻り値をタプルとディクショナリで戻す場合の速度を比較してみました。

一応私のPCのスペックです。

ノートPCのスペック : VAIO Pro PG - 趣味で計算流砂水理

  • タプル
def returntuple(num):
    A, B, C, D, E = float(0), float(0), float(0), float(0), float(0)
    
    for _ in range(num):
        A += 1.1
        B += 1.11
        C += 1.111
        D += 1.1111
        E += 1.11111
    
    return A, B, C, D, E
%%timeit
returntuple(10000)

1.27 ms ± 33.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

  • ディクショナリ
def returndic(num):
    A, B, C, D, E = float(0), float(0), float(0), float(0), float(0)
    
    for _ in range(num):
        A += 1.1
        B += 1.11
        C += 1.111
        D += 1.1111
        E += 1.11111
    
    return {'A':A, 'B':B, 'C':C, 'D':D, 'E':E}
%%timeit
returndic(10000)

1.32 ms ± 29.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

両者の速度はほとんど変わらない。

numba.jit最適化環境での速度比較

次にnumba.jit最適化時の速度を比較しました。

  • 参考:numba.jitとは

computational-sediment-hyd.hatenablog.jp

computational-sediment-hyd.hatenablog.jp

  • タプル
from numba import jit
@jit(nopython=True)
def returntuple(num):
    A, B, C, D, E = float(0), float(0), float(0), float(0), float(0)
    
    for _ in range(num):
        A += 1.1
        B += 1.11
        C += 1.111
        D += 1.1111
        E += 1.11111
    
    return A, B, C, D, E
%%timeit
returntuple(10000)

16.2 µs ± 5.84 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

  • ディクショナリ
@jit(nopython=True)
def returndic(num):
    A, B, C, D, E = float(0), float(0), float(0), float(0), float(0)
    
    for _ in range(num):
        A += 1.1
        B += 1.11
        C += 1.111
        D += 1.1111
        E += 1.11111
    
    return {'A':A, 'B':B, 'C':C, 'D':D, 'E':E}
%%timeit
returndic(10000)

42.6 µs ± 23.2 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

タプルの方が2.6倍早くなります。 大規模な配列だともっと差がでるようです。

まとめ

  • 戻り値の型は何でも大丈夫。気にしなくても良い。
  • しかし、numba.jitを使う場合等、速度を重視する場合タプルの方が良い。

  • 例の話はpix2pixでした。

lemファイルをgeotiffファイルに変換するpythonスクリプト

ユーザーがかなり限定されますが...

※ テストが十分でないのでエラー等指摘いただけると助かります。


lemファイルとは

  • lemファイルとは、航空LiDARの数値標高データの形式の一つです。

参照:https://www.harrisgeospatial.co.jp/Portals/74/VIS_Japan/documents/3_ENVI2017_CaseStudy_Toda.pdf?ver=2017-11-09-141636-000#page=24

  • lemファイルとその属性情報を示すcsvファイルがセット。
  • フォーマットは次を参照

https://gigaplus.makeshop.jp/jmcnet/download/Form2m.txt

pythonスクリプト

  • lem2geotiff.pyの関数translateが本体です。
  • 処理の流れは以下のとおりです。

    1. csvファイルより必要情報を読み取る。
    2. lemファイルを読み、ndarrayに変換する。
    3. ndarrayからxarrayを作成する。
    4. xarrayをgeotiffとして出力
  • 引数は次の4つのみ。filepath以外はオプショナルです。

    • filepath : str - lemファイルのフルパス
    • valm1111 : float, default 60000 - データ作成範囲外の場合のフラグ-1111の値
    • valm9999 : float, default 50000 - 海部及び陸水部の場合のフラグ-9999の値
    • dtype : type, default np.uint16 - geotiffファイルのデータ型
def translate(filepath, valm1111=60000, valm9999=50000, dtype=np.uint16 ):
    fpath = Path(filepath)
    
    # reading csv file
    df = pd.read_csv(fpath.with_suffix('.csv'), encoding='ShiftJIS',header=None, index_col=0)
    
    nx = int(df.loc['東西方向の点数'].values[0])
    ny = int(df.loc['南北方向の点数'].values[0])
    y1 = float(df.loc['区画左下X座標'].values[0])
    x0 = float(df.loc['区画左下Y座標'].values[0])
    y0 = float(df.loc['区画右上X座標'].values[0])
    x1 = float(df.loc['区画右上Y座標'].values[0])
    dx = float(df.loc['東西方向のデータ間隔'].values[0])
    dy = float(df.loc['南北方向のデータ間隔'].values[0])
    nepsg = int(df.loc['平面直角座標系番号'].values[0])
    
    # reading lem file
    with open(filepath, "r") as f:
        lines = f.readlines()
    
    # Important: The dimensions of the array must be in the order from y to x.
    arr = np.empty((ny, nx))
    for i, l in enumerate(lines) :
        ll = l.replace('\n','')
        s = ll[10:]
        v = [s[i: i+5] for i in range(0, len(s), 5)]
        arr[i,:] = np.array(v)
    
    # make xarray
    val = np.array(arr, np.float)
    
    # Changing the flag -9999 for the sea and land and water sections and the flag -1111 outside the data creation area into any given value.
    # dafault : -9999 to 50000, -1111 to 60000
    val = np.where(val == -1111, valm1111, val)
    val = np.where(val == -9999, valm9999, val)
    
    # changing data type(default uint16)
    val = val.astype(dtype)
    
    dx *= 100
    dy *= 100
    xarr = np.arange(x0,x1,dx)
    xarr = xarr + dx/2
    yarr = np.arange(y0,y1,-dy)
    yarr = yarr - dy/2
    
    epsg = str(6668 + nepsg)
    
    ds = xr.Dataset({'z': (['y','x'], val) }, coords={'x': xarr/100, 'y': yarr/100}) 
    ds = ds.rio.write_crs('EPSG:' + epsg, inplace = True)
    # ds.rio.reproject("epsg:****")
    
    # export geotiff file
    ds['z'].rio.to_raster( fpath.stem + '.tif')

Usage

Exmaple 1

単一ファイルの場合

import lem2geotiff as l2g

l2g.translate('lem/****.lem')

Exmaple 2

複数ファイルの一括処理

import lem2geotiff as l2g
import glob

for ff in glob.glob( 'lem/*.lem'):
    l2g.translate(ff)
おまけ

lemファイルは膨大なファイル数になることが多いので、ソフトウェアで取り扱う場合はVRTファイルを使うことが多い。

from osgeo import gdal 

my_vrt = gdal.BuildVRT('output.vrt', glob.glob( '*.tif'))
my_vrt = None
  • gdalの場合
gdalbuildvrt out.vrt *.tif

Environment

rasterとxarrayの変換にはrioxarrayモジュールを使用している。

詳細を省略するが、gdalを使用するよりもコード量が大幅に少なくなるため、非常に便利である。

Anaconda

仮想環境lemを作成する。

conda create -y -n lem python=3.7
conda activate lem 
conda install -y -c conda-forge rioxarray
conda install -y -c conda-forge gdal
conda install -y -c conda-forge jupyter

Google Colaboratory

  • rioxarrayをインストールする。
  • ドライブをマウントして作業ディレクトリに移動する。
!pip install rioxarray
from google.colab import drive
drive.mount('/content/drive')
%cd "/content/drive/My Drive/*"

github

github.com

参考サイト

raster - Python equivalent of gdalbuildvrt - Geographic Information Systems Stack Exchange

[Python]可読性を上げるための、docstringの書き方を学ぶ(NumPyスタイル) - Qiita

docstringのstyle3種の例 - Qiita

numpydoc docstring guide — numpydoc v1.2.dev0 Manual

numba.jitによる高速化:クラス編

  • numba.jitによる高速化のクラス編です。公式にもほとんど情報が無いので結構苦労しました。

※2024/2/25追記:以下も分かりやすいです。

qiita.com

  • numba.jitを知らない方はまず関連記事を読んだほうがスムーズだと思います。

関連記事

computational-sediment-hyd.hatenablog.jp

computational-sediment-hyd.hatenablog.jp

過去記事 と同様に半径1の円を簡素化した64角形の面積を求める問題を例とします。


computational-sediment-hyd.hatenablog.jp


from shapely.geometry import Point
import numpy as np
p = Point(0.0, 0.0)
x = p.buffer(1.0)
s = x.simplify(0.001, preserve_topology=False)
polygon = np.array( [ p for p in s.exterior.coords])
len(polygon)
# 65
s

def Area(p):
    s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
    return 0.5*np.abs(np.sum(s))

Area(np.array(polygon)) 
# 3.1365484905459393

numba version

numbaのバージョンを確認する。この手の話はまだまだ仕様変更があると思われます。

import numba
numba.__version__

'0.46.0'

PCのスペック

一応速度比較を行っているので。

ノートPCのスペック : VAIO Pro PG - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

単純なクラスの高速化

  • 面積を求めるメソッドAreaを持つクラスshapeを定義する。
  • クラス変数の型を指定して、次のようにクラスの前に@jitclass(クラス変数の型)をつける。
from numba import jitclass   
from numba import float64

spec = [
     ('polygon', float64[:,:])
]
@jitclass(spec)
class shape(object):
    def __init__(self, polygon):
        self.polygon = polygon 
        
    def Area(self):
        p = self.polygon
        s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
        return 0.5*np.abs(np.sum(s))
%%timeit
c = shape(polygon)
c.Area()
  • 参考サイト

Compiling Python classes with @jitclass — Numba 0.50.1 documentation

速度比較:%%timeで計測

numbaにより約10倍高速化する。

  • numbaあり

    13.2 µs ± 229 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

  • numbaなし

    119 µs ± 614 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

ネストされたクラスの高速化

  • 先程のクラスshape内に新たなクラスmodelを設ける。
  • クラス変数にクラスを与える場合、型はxxxx.class_type.instance_typeと指定する。
  • 参考サイトでは、deferred_typeが必要と書いてあるが無くても問題無さそう。
from numba import jitclass   
from numba import float64

spec = [
     ('polygon', float64[:,:])
]
@jitclass(spec)
class model(object):
    def __init__(self, polygon):
        self.polygon = polygon 
        
    def Area(self):
        p = self.polygon
        s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
        return 0.5*np.abs(np.sum(s))

spec = [
     ('model', model.class_type.instance_type)
]
@jitclass(spec)
class shape(object):
    def __init__(self, polygon):
        self.model = model(polygon)
        
    def Area(self):
        return self.model.Area()
%%timeit
c = shape(polygon)
c.Area()
  • 参考サイト

python - How to nest numba jitclass - Stack Overflow

速度比較:%%timeで計測

numbaにより約10倍高速化する。

  • numbaあり

    13.6 µs ± 169 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

  • numbaなし

    122 µs ± 4.52 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

クラスのリストを含むクラスの高速化

  • クラスshape内のpolygonをクラスpointのリストとして定義する。
  • クラス変数にクラスのリストを与える場合、types.List(xxx.class_type.instance_type, reflected=True)と指定する。
  • reflected=Trueは必須。(詳細は理解していません。詳しい人教えて下さい。)
from numba import jitclass   
from numba import float64
from numba import types

spec = [
     ('coordinate', float64[:])
]
@jitclass(spec)
class point(object):
    def __init__(self, p):
        self.coordinate = p

spec = [
     ('polygon', types.List(point.class_type.instance_type, reflected=True))
]
@jitclass(spec)
class shape(object):
    def __init__(self, polygon):
        self.polygon = [point(polygon[i]) for i in range(len(polygon))]
        
    def Area(self):
        p = self.polygon
        s = np.array( [ p[i-1].coordinate[0]*p[i].coordinate[1] - p[i].coordinate[0]*p[i-1].coordinate[1] for i in range(1, len(p)) ] )
        return 0.5*np.abs(np.sum(s))
%%timeit
c = shape(polygon)
c.Area()
  • 参考サイト

https://groups.google.com/a/continuum.io/forum/#!topic/numba-users/MygFRmKlr24

速度比較:%%timeで計測

numbaにより約6倍高速化する。遅い書き方をしているが、numbaなしの場合あまり遅くならない。numbaありの場合は上の2ケースより遅くなる。

  • numbaあり

    22 µs ± 519 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

  • numbaなし

    115 µs ± 2.58 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

調査中:継承

継承された側のクラスのみにnumbaをつける。

super()メソッドが使えないっぽい。オーバーライドは問題なく使える。

Gist

Jupyter Notebook Viewer


大きなモデルを作る場合はクラスは必須なのでnumbaが使えるとpythonの可能性が広がります。

球磨川流域の雨量について

湯前町の雨量がすごいということだったので、湯前横谷(気象庁の地上雨量観測)のデータをみてみました。 累加雨量が500mmと半端じゃない雨ですが、近年豪雨続きであまり驚かないです。 台風じゃないので予測しづらいのが厄介ですね。

http://www.river.go.jp/kawabou/graphOut/41/1621400100481/20200705/100000_60_S1_D48.png

参考までに東京にこのような雨が降ることを想定してみました。 都市部の排水を考えると短時間降雨が重要となるので元データから10分雨量の時間変化を書いてみました。

f:id:SedimentHydraulics:20200705153654p:plain

現在の東京の排水能力は区部の7割が50mm/hr対応(10分雨量で8.33mm/hr)、一部浸水リスクが高い地域は75mm/hr対応(10分雨量で12.5mm/hr)になっているとのことです。 最終的には区部全域で75mm/hr対応になるとのことですが、これだけの雨が降ると駄目ですね。怖いなあ。

球磨川の浸水推定図が公開されたので洪水ハザードマップと比べてみた

球磨川の浸水推定図が公開されたので洪水ハザードマップと比べてみました。

  • 浸水推定図(2020年7月4日20時作成)

今回の洪水が想定最大規模ほどは大きくないと思うのでこんな感じなのかな。

沖積河川の砂州形成の数値実験

今回の記事はブログレベルの話ではないですが、ずっと悩んでいるのでとりあえずまとめておこうかなというところです。


砂州の定義は以下の記事を参照下さい。

computational-sediment-hyd.hatenablog.jp


支配方程式

流水

  • 平面二次元水深平均流モデル

流砂、河床変動

  • 河床位方程式
  • 掃流砂のみ、単一粒径、平衡流砂量式MPM式
  • 限界掃流力に局所勾配の影響を考慮する。砂州の計算では最も重要な要素でこれがないと砂州が計算できない。 芦田, 江頭, 劉の方法を基本に若干の修正を加えた。
  • 斜面崩壊の影響を考慮する。後述するCASE3のみに影響する。 崩壊の向き、量の考え方は関根の方法を基本とするが、崩壊量は流砂量に付加するのではなく、瞬時に崩壊が生じるものとした。

計算条件

斐伊川下流部をモデルに設定した。初期水深を変えた3ケースの計算を行った。

  • 河床勾配:1/670

  • 河床材料:2mm

  • 川幅:90m

  • 延長:1800m

  • 格子サイズ: 5m

  • 境界条件

  • 初期条件

    • 河床:最大値を粒径とした乱数を与える
    • 初期水深
      • CASE1 : 1.0m
      • CASE2 : 0.6m
      • CASE3 : 0.3m
    • 無次元掃流力:
      • CASE1 : 0.452
      • CASE2 : 0.271
      • CASE3 : 0.136
  • 計算ケースごとの川幅水深比と無次元掃流力の関係を黒木, 岸 の線形安定化解析に基づく砂州形態の分類図にプロットした。

  • CASE1は単列砂州、CASE2、3は複列砂州の領域となる。

計算結果と考察

砂州形成開始時の水深、河床変動量の平面空間分布

  • 計算開始から300時間後の計算結果を下図に示す。
  • 砂州の列数は、CASE1:1列、CASE2:2列、CASE3:3列以上となっており、卓越する成分が明確に把握できる。
  • この状態が上記領域区分の示すところ?(与えた擾乱のうち、どの成分が卓越するか)
CASE1(拡大図はこちら

CASE2(拡大図はこちら

CASE3(拡大図はこちら

概ね平衡状態時の水深、河床変動量の平面空間分布

  • 計算開始から2000時間後の計算結果を下図に示す。
  • CASE1およびCASE2ともに、砂州の列数が1列となっている。一方、CASE3は周期性を持たない複雑な流路を示す。
  • CASE2は初期は列数が2列であったが時間の経過とともに個々の砂州が干渉し、1列の砂州が形成される。
  • CASE3では、初期は複数の流路が出現したが、浮州(ドライベッド)の出現をきっかけに流路の変化が複雑化し、最終的には複数の流路が粗発生する。
CASE1(拡大図はこちら

CASE2(拡大図はこちら

CASE3(拡大図はこちら

砂州形成過程:gif動画

CASE1(拡大図はこちら
  • 500時間後には平衡状態に達してそれ以降は周期的な変化を示す。

CASE2(拡大図はこちら
  • 1000時間後に状態に達する。CASE1と比べて砂州の移動速度が遅い。
  • 2列の砂州が維持できない。なぜ?

CASE3(拡大図はこちら
  • 5700時間後までの計算を実施しているが平衡状態に達していない。周期的な変動ではないため、平衡状態は存在しない可能性がある。
  • 流路変動は側岸浸食によって生じるため、その速度は非常に遅い。
  • もう少し分析が必要。さらに長期の計算を実施予定。

まとめと課題

  • 任意の条件下における砂州の平衡状態を数値計算によって検証した。
  • CASE1は形成条件のとおり、単列砂州が形成される。
  • CASE2は複列砂州条件では安定的な複列砂州は形成されず、長期的にはCASE1と同様に単列砂州となる。では、斐伊川のような多列砂州はなぜ維持されるのか?非定常の影響?他の乱れの影響?もう少し検討が必要。
  • CASE3は、浮州が出現、つまり砂州高>水深となり、複雑な地形を示す。個人的には、これが流路=河川が形成される機構ではないかと考えている。さらにテスト計算を実施予定。

GitHub

Jupyter Notebook Viewer


図の作成に使用したライブラリは以下です。

  • hvplot:計算結果の図化

computational-sediment-hyd.hatenablog.jp

また、はてブロへのhtmlグラフの埋め込みは以下を参照。

computational-sediment-hyd.hatenablog.jp


ご意見、質問等お願いします。

沖積河川の砂州について整理してみる

砂州についてしゃべることは多いのですがまとめたことが無かったので整理しておこうかと。 面白い数値計算結果が出たので、その準備だったりもしもます。


はじめに

  • 沖積河川の砂州(以降、砂州)とは、 川の流れに伴う土砂移動によって形成される川幅スケールの河床地形を示す。
  • 川幅(堤防間の幅)やその曲率等の平面的な河川地形の影響により砂州の形状が異なる。
  • 日本の土砂水理関連の書籍には砂州に関する項目は必ず記述されるが、海外の書籍には記述が少ない。理由は不明。

砂州の分類

砂州の分類

日本では土木学会水理委員会移動床流れの抵抗と河床形状研究小委員会 で定義された次の分類を用いることが多い。

「土砂水理学1 POD版」 p.102より引用

また、河村は次のように整理している。

「土砂水理学1 POD版」 p.101より引用

砂州の特徴

実河川を例にそれぞれの砂州形状の特徴を簡単に示す。

交互砂州(単列砂州)

左右交互に周期的に形成される砂州を示す。下流方向に伝播する。

多列砂州

上述の交互砂州が河道内に2つ以上形成される砂州形状を示す。

網状流路

洪水時に大きく変動しない(砂州上の流れの掃流力が小さい)砂州が部分的に出現することにより、周期性を持たない網状流路が形成される。

画像を右クリックで再生を選択

固定砂州

湾曲部にみられる砂州で経年的に形状が変化しない。砂州を形成する土砂は動いているが砂州形状は変化しない。

  • 例:太田川(左:2008年、右:1961-69年)

砂州の形成機構

  • 平面的な流れ場の影響によって砂州が形成される。
  • 流砂は局所的な地形勾配の影響を強く受ける。
  • 固定砂州を考える場合は、鉛直方向の二次流の流砂への影響を考慮する必要がある。

  • 黒木, 岸 の線形安定化解析より砂州形状は川幅と水深の比によって決定することがわかる。

「21世紀の河川学」p.119より引用

おわりに

  • サステイナブルな川とのつき合い方を考える上で砂州の概念は非常に重要です。
  • 砂州は流水断面を阻害するため、水位を下げる観点からは無い方が良いが、掘削してもわずか数年で元に戻ってしまいます。膨大な費用をかけて掘削したのに次の年には元に戻った事例もあります。
  • また、砂州の環境的な価値(生態系、景観、親水等)も重要視されております。
  • 近年の観測技術の向上により、水面下にも砂州が形成されていることが確認されており(例えば、秋田:水中の河床地形の面的計測とその活用方策について等)、今後より知見が増えていくと思われる。

GitHub

Jupyter Notebook Viewer