趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

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

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

国土地理院の標高タイルから標高値を取得する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

またまた地図ネタです.

先日ネタにでましたミャンマーのある河川の河口部について,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

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

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

ずっと悩んでいるのですが,改めて.

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

nbviewer

問題は下流端流量です.

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

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

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

上流端の河積(水深)も同様の問題が生じます.

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


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

1DUpWindFlux_schemecheck

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

今まで,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

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

github.com

またまた更新です.

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

計算結果

case1

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

大きいものはこちら

case2

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

大きいものはこちら

case3

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

大きいものはこちら

case4

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

大きいものはこちら

感想

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

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

先日お話させて頂きました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))も駄目でした.

lambda式

駄目です.ちょっと不便

pandas

これも駄目です.

dictionary型

駄目です.

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

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

感想

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

rigid-lid近似

連投になります.

直交格子でLES + 河床変動計算をやりたいと考えており, iricにあるNabiさんのLESモデルについて調べていると,河床変動は出来ないのですが,関連する論文を調べているとすでに河床変動計算も完成しているようです.

この辺りです.

https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2012WR011911 https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/wrcr.20303 https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/wrcr.20457

水面の処理の話が殆ど書いていないので,よくよく読んでいるとrigid-lid近似としていると書いてありました.

これって,海の計算でよく使う「海面の鉛直変位は認めない」という条件です.要するに,水位は一定ということです.

河床変動って,河床面の摩擦力なので,垂直抗力(重力)が重要なので,このような近似を設けても良いのかなと疑問に思います.