趣味で計算流砂水理

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

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

MENU

オープンデータによる河川流解析：平面二次元流解析

R

ようやく今年1つ目の記事です．

今年はオープンデータを活用したフィールドの計算をやろうとかなと思ってます． 第一弾は，平面二次元流の計算結果です．まだ，ソースは上げてないですが週末には上げます．

斐伊川ですがなかなかいい感じです．実はドライベッドで質量が保存しないところに結構苦労しました．

大きいものはこちら

2019/2/4 uplaod

github.com

徐々に解説を書きます．

www.xsim.info

openfoamの条件ファイルを作るwebアプリですがご存知でしたか？

こんなのが作りたいです．

地理院標高タイルから標高値を取得するpythonモジュール : libgsidem2el

R

またまた地図ネタです．今回はちょっと力作です．コードは短いですが，結構便利です．

国土地理院の標高タイルから標高値を取得するpythonモジュール libgsidem2el を作ってみました．

標高タイルとは， https://maps.gsi.go.jp/development/ichiran.html#dem を見るとわかるように，タイル単位で配信される数値標高データです．

これがあると，基盤地図情報からダウンロードしなくても標高値を取得できます．

コードは以下です．

github.com

Usage

import module

import libgsidem2el as gsi

constructor

dem = gsi.libgsidem2el('DEM5A')

標高値を取得

  • 引数は経度，緯度，ズームレベルです.
  • ズームレベルは，DEM5A,B:15, DEM10B:0-14, DEMGM:0-8です．
el = dem.getEL(lon, lat, zoom = 15)
  • getELメソッドで取得，DEM5AB，10Bは[m],DEMGMは[cm]の標高値を返す．
  • それぞれの地図の対象範囲外は'outside'を返す
  • 水面等により欠測の場合は，'e'を返す（地理院タイルの使用）.
  • 一度読み込んだタイルはクラスの変数demdfに保存するため，近接するデータの取得が高速化される反面，大領域を一度に処理する場合はメモリを消費する．

これを使って計算の準備をしています．お正月休みには準備できると思います．

ミャンマー某河川の河口部：ETOPO1 + geoviews

R

またまた地図ネタです．

先日ネタにでましたミャンマーのある河川の河口部について，ETOPO1とgeoviewsで絵を書いてみました．

大きいものはこちら

  • スケールがわかりにくいですが，河口から300km程度は，10m程度の浅い水深の範囲が続いています．
  • 河口から50kmくらいがぐちゃぐちゃな感じです．この辺りが活発に河床変動が生じているような気がします．
  • スケールが大きくて面白いです．
  • landsatのデータもかなり面白いので，後日絵を書いてみます．

ETOPO1 + geoviewsの作り方

驚くほど簡単でした．手順をまとめておきます．

モジュールのロード

以下のとおり，geoviewsとxarrayを読みます．

xarrayは多次元データの解析用のモジュールで非常に高速です．後日記事を書きます．

import xarray as xr
import geoviews as gv
import geoviews.tile_sources as gts

gv.extension('bokeh')

ETOPO1をロード

以下から，ETOPO1のnetCDFデータをダウンロードします．これが一番時間がかかります．

https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/

リードは一行です．gzのままいけちゃいます．

xr_ensemble = xr.open_dataset('ETOPO1_Ice_g_gmt4.grd.gz')

データセットの編集

  • 系列名を変更:geoviewsで複数の地図をリンクする際に軸名が同じである必要があるため．
  • 描画範囲をスライス：全部は重すぎます．
xr_ensemble = xr_ensemble.rename({'x':'Longitude', 'y':'Latitude'})
xr_ensemble = xr_ensemble.sel(Longitude=slice(95,99), Latitude=slice(14,18))

xarrayをgeoviews.Datasetに変更

この形式以外でも描けますが，geoviews.Datasetは高速なので，データが多い場合は有効です．

kdims = ['Longitude', 'Latitude']
vdims = ['z']
xr_dataset = gv.Dataset(xr_ensemble, kdims=kdims, vdims=vdims)

レンダリング

今回はESRI mapとリンクさせてます．

%%output size=100 filename='out' fig='html'
%opts Image {+framewise} 
xr_dataset.to(gv.Image, ['Longitude', 'Latitude']).redim.range(z=(-50, 10)).options(cmap='jet',colorbar=True,width=450, height=400) \
+ gts.ESRI.options(width=400, height=400) * gts.StamenLabels.options(level='annotation')

ソースコードは以下です．

github.com

ETOPOメモ

備忘録です．

ETOPO 5/2/1について

データセットの解像度を示す.ETOPO1が最新.

  • ETOPO5 : 5分 緯度5分=~9.2km
  • ETOPO2 : 2分 緯度2分=~3.7km
  • ETOPO1 : 1分 緯度1分=~1.8km

http://ofgs.aori.u-tokyo.ac.jp/~okino/gmtscripts/datainfo.html

Ice Surface or Bedrock

https://www.ngdc.noaa.gov/mgg/global/etopo1sources.html

南極とグリーンランドの標高値を氷上か基岩かどちらを使うかの差だけです．基本は，絵的には前者を使ったほうが良いですね．

コロケート格子の境界条件の設定の悩み

R

ずっと悩んでいるのですが，改めて．

非常にシンプルな一次元不定流計算の定常解です． バックウォーターの計算です．

nbviewer

問題は下流端流量です．

境界なので片側差分で計算にすると，少し大きな値になってしまいます．

圧力項dHdxの差分がまずいです．

単純に一つ上の断面と同値にすると問題ないのですが，非定常の場合はあまりよろしくないですよね．

上流端の河積（水深）も同様の問題が生じます．

一般的にこのような場合はどうやって解決するのでしょうか．．．．高次精度とか？

なんかgistがうまく表示されない．

1DUpWindFlux_schemecheck

はてブロにインタラクティブなグラフを埋め込む：外部htmlのとりこみ　CDNとかiframeとか

R

今まで，javascriptを直で書いたりとかしてましたが，よく考えたら，htmlをそのままiframeで埋め込むのが楽じゃないかという話です．

結果簡単にできました．

基本的な流れは，

  1. Holoviewsでグラフを作成し，htmlにエクスポート
  2. githubに公開
  3. github上のhtmlをCDN化
  4. ブログ上にiframeで取り込み

1.2.は常識なので，3.4.をメモしておきます．

その前に

CDNサーバrawgitが終了

愛用していたRawGitが終了するとのことです.新たなアップロードは不可でこれまで上げたものは，2019/2までは参照可能とのことです．これまでブログに上げたものにも影響がでますね．．．

rawgit.com

github上のhtmlの公開方法：github pagesの活用

以下を参照．簡単です．

[setting] - [GitHub Pages] - [source] : master branch に変更でsave

www.tam-tam.co.jp

変更が10分くらい反映されません．

いろんなサイトにindex.htmlを作らないと駄目だよと書いてますが，無くても全然大丈夫でした．

  • index.htmlありの場合 : https://[username].github.io/[repositoryname]
  • index.htmlなしの場合 : https://[username].github.io/[repositoryname]/[htmlfilename]

前回の記事の内容をリンクすると，以下の感じです．

https://computational-sediment-hyd.github.io/1DBedVariationMixed/case01.html

iframeとして埋め込む方法

以下を参照．

mmorley.hatenablog.com

上記サイトのとおりですが，普通に埋め込むと，

<iframe src="https://computational-sediment-hyd.github.io/1DBedVariationMixed/case01.html" width="600" height="300" ></iframe>

なので，下のようにスケールを調整．大きめのサイズで取り込んで0.6倍してます．

<div style="width:600px;height:200px;"><div style="width:1000px;"><iframe src="https://computational-sediment-hyd.github.io/1DBedVariationMixed/case01.html" frameborder="0" width="1000" height="333" style="transform:scale(0.6); transform-origin: left top;"></iframe></div></div>

もっと複雑なグラフを作ろうと思うとサーバが必要ですね．

githubの新機能GitHub Actionsがベータ版ですが公開されたようで，自動でデプロイまでやってくれるらしいので，サーバを立てるのもいいですね．

www.atmarkit.co.jp

一次元混合粒径掃流砂河床変動計算の更新

R

github.com

またまた更新です．

ほぼ完成かなと思います．

計算結果

case1

勾配変曲点設けた計算です．

大きいものはこちら

case2

勾配変曲点設けた計算で,上流側河道をガチガチのアーマーコートにしてます．

大きいものはこちら

case3

河床勾配は一定で，狭窄部を中央に設けた計算です．川幅100mから50mです．

大きいものはこちら

case4

下流端水位をE.L.5mで固定してます．

大きいものはこちら

感想

  • 混合粒径の計算は，大粒径により変動が制限されますが，粒度分布は顕著に変わるところが面白いです．
  • 前回記事のとおり，速度が早くなったのでいくらでも回せます．
  • Holoviewの癖が強い．ブログ表示用にHTMLにする際は，データを間引かないと駄目です．いろいろ読んでるとBokehサーバを立てるともっと高度なものが作れるみたいです．今度やってみます．
  • ブログへの埋め込み等は後日記載します．
  • 次は浮遊砂ですね．

numba.jitの衝撃：5分 ⇒ 3秒

R

先日お話させて頂きましたnumba.jitの件です．私の手続きが間違ってました． 前回の河床変動の計算時間ですが，タイトルのとおり，超高速化します．100倍です． github更新しておきました．

github.com

numba.jitの使い方

でも，実問題はそんなわけにはいきません．

以下のサイトを参考にしました．

http://yutori-datascience.hatenablog.com/entry/2014/12/10/123157

http://www7a.biglobe.ne.jp/~thor/novel/column/14.html

ポイントは，公式 http://numba.pydata.org/ にもあるように，

@numba.jit(nopython=True)
or
@numba.njit

と書くことです．nopythonオプションは，最適化できないものはエラーとして返すというものです．これがないと最適化できるものはするけど，できないものはしないのでほんとに最適化できているかわからないということです．

特に，後述の高速化できないものがループに含まれるとそのループ自体が最適化されないので，ほとんど早くなりません．

一応並列化もやってみましたが，一次元なので，逆に遅くなりましたが，大規模だと早くなるようです．

@numba.jit(nopython=True, parallel=True)
# 並列化loopをrange(imax)から，
numba.prange(imax)
# に変更

まだ，やってないですがclassの最適化もあるようです．ちょっと面倒です．

https://stackoverflow.com/questions/38682260/how-to-nest-numba-jitclass

numba.jitで使えないもの

文字列型

一切駄目です．例えば，print(str(i))も駄目でした． 引数は駄目です。printは使えます。

lambda式

駄目です．ちょっと不便 使えます。

pandas

これも駄目です．

dictionary型

駄目です．

通常のリスト型およびその処理

すべて，numpyにしましょう．例えば，通常sumも使えません．np.sumならO.K.です

感想

  • もとのpythonのコードをほとんど手を加えずに，超高速化出来るので便利です．野望が広がります．
  • pythonらしいコードは書けないので，必要なdefだけ最適化すれば良いと思います．それでも，リスト内包表記や配列の演算などは普通に使えると十分です．
  • なので，静的な言語の知識は若干必要かも．我々のように高級言語から入った人間にはそんなに苦労する話じゃないですが．
  • numbaのエラー文がめちゃくちゃわかりにくい．エラー箇所が特定できないので，ちょっとコーディングが大変.

2020/3/7追記

より詳しい解説を下記事にしました．

computational-sediment-hyd.hatenablog.jp

computational-sediment-hyd.hatenablog.jp

2020/06追記

computational-sediment-hyd.hatenablog.jp