趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

衛星画像処理用スクリプト:備忘録

最近衛星画像の処理をする機会が増えましたが,データ容量が大きいため,GISソフトウェアなんか使っていたらいつまでも終わらないので,GDALやpythonで処理をします。

結構煩雑なので備忘録的にまとめておきます。

ちなみにフリーで手に入るものは,

ALOS AVNIR-2 https://www.eorc.jaxa.jp/ALOS/alos-ori/index.html

LANDSAT8(他にもいろんなダウンロードサイトがあります) https://landbrowser.airc.aist.go.jp/landbrowser/index.html

等があります。

回転したgeotifのRGB合成

衛星画像は緯度に対して回転した(傾いた)geotiffがほとんどですが,gdalが対応してません。なので,worldfileを使って以下のようなに処理が必要。

rem worldfile(tfw file)を作成
gdal_translate -co "TFW=YES" band.tif out.tif

rem worldfileを設定(作業終了後は必ずunsetする)
set GDAL_GEOREF_SOURCES=out.tfw
rem vrtファイルを作成
gdalbuildvrt -separate out.vrt bandred.tif bandgreen.tif bandblue.tif
rem RGB合成したtifを作成
gdal_translate out.vrt out.tif

参考サイト

https://gis.stackexchange.com/questions/271995/how-to-get-gdaltranslate-to-create-world-file

https://stackoverflow.com/questions/52334933/merging-multiple-bands-together-through-gdal-correctly

https://gis.stackexchange.com/questions/336367/error-gdalbuildvrt-does-not-support-rotated-geo-transforms/336674

pansharpen

LANDSAT等は,パンシャープンを行う。

gdal_pansharpen panchro.tif bandred.tif bandgreen.tif bandblue.tif out.tif

参考サイト

https://gdal.org/programs/gdal_pansharpen.html

大規模tifのmerge

基本的には使わない。gdalbuildvrtでvrtを作る。

圧縮とかbigtiffとかを指定。

gdal_merge -o out.tif -ot Uint16 -co COMPRESS=LZW -co BIGTIFF=YES file1.tif file2.tif ......

参考サイト

https://qiita.com/naru-T/items/b87c83fad05634919532

ラスタ演算

pythonがかなり楽です。

import xarray as xr

# この一行でgeotiffをxarrayで読み込む
d = xr.open_rasterio('*.tif')
# バンドを取得,シングルでもマルチでも同じ
band = d.values

# こんな感じで計算
ind = (band[2] - band[4]) / (band[2] + band[4])

ちなみに結果をgeotiffで出力する場合,rioxarrayを使う。

import rioxarray

xrds = xr.Dataset({'value':(['y','x'], ind)} 
                    , coords={'x':d['x'], 'y':d['y']}
                    , attrs= d.attrs ) 

# うまくepsgが付かないことがある。
xrds.rio.set_crs("epsg:****")
xrds['value'].rio.to_raster('out.tif')

tellus(https://www.tellusxdp.com/ja/)がローカルのjupyterでも使えるようになったらいいのな。

デカルト座標系による実河川の平面二次元計算

実河川の平面二次元は一般座標系又は直交曲線座標系で解くことが一般的ですが,格子生成が面倒です。

そこで,デカルト座標系で計算してみました。

計算条件

  • 境界に合うようにメッシュを回転させました。
  • メッシュサイズは2.5x2.5m。メッシュ数は10万程度
  • 標高は地理院の5mメッシュでいつものlibで作りました。

computational-sediment-hyd.hatenablog.jp

  • 初期はドライベッドで上流端から10m3/sを流して定常解を得るまでの2時間計算しました。

計算結果

  • かなりいい感じです。計算時間も21分で実用的かと思います。まあ,メッシュの9割は無駄ですが。。。
  • ちょっと雰囲気を見る程度の計算には便利そうです。

ソースコードとか

github.com

nbviewerは以下です。

流れの計算

https://nbviewer.jupyter.org/github/computational-sediment-hyd/D2Dmodel-of-NaturalRiver-in-CartesianGrid/blob/master/D2Dmodel.ipynb

格子生成

https://nbviewer.jupyter.org/github/computational-sediment-hyd/D2Dmodel-of-NaturalRiver-in-CartesianGrid/blob/master/make-elevation.ipynb

図化

https://nbviewer.jupyter.org/github/computational-sediment-hyd/D2Dmodel-of-NaturalRiver-in-CartesianGrid/blob/master/make-figure.ipynb

河川断面データをgeojsonで取り扱う

河川断面データは河川のシミュレーションを行う上で最も重要データですが,geometry情報として扱われておらず,計算の際に非常に手間です.

そこで,geojsonで取り扱う方法を考えてみました.

対象として,荒川下流部の以下の測線を選びました.

例のツールで地理院5mメッシュより断面を生成.

computational-sediment-hyd.hatenablog.jp

※このツールは結構良い出来だと思ってましたが,こちら国土数値情報の鉄道データを3Dで表示してみる - Qiitaのほうが遥かに上手でした.勉強します.

手順はギハブに上げますが,測線上に左岸から5m間隔で配置した点群に標高値を付与したxyzの点群データを基にgeojsonを作成します.

点群のndarrayは,

point3d
# ---------------------------------------------------------
array([[-3.36773568e+03, -2.63076896e+04, -1.20000000e-01],
.........................................
       [-3.29085113e+03, -2.69732637e+04,  6.00000000e-02]])

geojsonは以下のような感じでLineStringで扱うのがいい気がします. multipointでもいいですがこちらのほうが扱いやすいです. sections.geojsonでエクスポートしておきます.

from shapely.geometry import mapping, LineString

j = {
"type": "FeatureCollection",
"name": "area",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::6677" } },
"features": [
{ "type": "Feature", "properties": { "id": 1 }, "geometry":mapping( LineString(point3d) ) }
]
}

with open("sections.geojson", "w") as f:
    json.dump(j, f)
# ---------------------------------------------------------
{'crs': {'properties': {'name': 'urn:ogc:def:crs:EPSG::6677'}, 'type': 'name'},
 'features': [{'geometry': {'coordinates': ((-3367.73567716,
                                             -26307.68964812,
                                             -0.12),
                                            (-3367.161911906202,
                                             -26312.65661836689,
                                             -0.09),
....................................
                                            (-3291.4248984048904,
                                             -26968.296690956955,
                                             0.21),
                                            (-3290.8511331510927,
                                             -26973.263661203848,
                                             0.06)),
                            'type': 'LineString'},
               'properties': {'id': 1},
               'type': 'Feature'}],
 'name': 'area',
 'type': 'FeatureCollection'}

wktで見ると,LINESTRING Zになります.

LineString(point3d).wkt
# ---------------------------------------------------------
'LINESTRING Z (-3367.73567716 -26307.68964812 -0.12, -3367.161911906202 -26312.65661836689 -0.09,........................................)

例えば,断面図を図化するときは以下のような感じです. geopandasを使うのは趣味ですがワンライナーでgeojsonが読み込めます.

import geopandas as gpd
import holoviews as hv
import numpy as np

hv.extension('bokeh', logo=False)

gdf = gpd.read_file("sections.geojson")
point = np.array( gdf.geometry[0].coords[:] )
# drop z
pointxy = point[:,:2]
L = np.linalg.norm( pointxy - pointxy[0],ord=2,axis=1 )
z = point[:,-1]
f = hv.Curve((L,z)).options(width=400, height=300)

愛用のholoviewsで簡単に書けます.

次は数値計算編を書きます.

githubは以下です.

github.com

nbviewerは以下です.

Jupyter Notebook Viewer

平面二次元計算による急拡水路の河床変動

久しぶりに完全に趣味で計算しました.

急拡水路の河床変動計算を平面二次元でやってみました.

計算条件

  • 掃流砂のみ,M.P.M.式
  • 単一粒径 5mm 無次元掃流力は1くらい
  • 初期水深は上流端で1mくらい.上流端流量は一定,平衡給砂
  • 下流端は等流
  • dt:0.02sec,dx,dy:5m
  • 計算時間は500時間

計算結果

初期 f:id:SedimentHydraulics:20190803164708p:plain 500時間後 f:id:SedimentHydraulics:20190803164731p:plain

詳細な計算結果はこちら

  • まだまだ平衡状態には至ってない感じですね.
  • 河道中央に想像以上に堆積します.
  • 以下の記事による簡易式だと拡幅部の堆積量は35cmになります.もうちょっと流し続ければいいのかな.
  • ソースとかはgithubに上げました.

github.com

その他

  • 計算時間は私のPC(後述)で23時間かかりました.
  • CFLはかなり余裕を持ってます.流れだけだとdtを5倍くらいにできます.でも,河床変動の場合,大きくするとたまにドライベッドがでるのでそのときに跳ねてしまいます.上手くdtを可変にすればかなり速度はあがると思います.
  • 部分的にはかなりの急勾配がでてます.崩落モデルをいれると変わるかも.
  • 乱流モデルによっては剥離域の形状は変わるかも.

新しいPC(ubuntu18.04 インテル® Core™ i9-9900K)が引くほど早い.ノートPCのi7-8550Uと比較すると,シングルコアの計算で40%程度速くなります.並列化が楽しみです.

周期境界の三重対角:Sherman-Morrison formula

2019/7/31にあげた記事に間違いがあり修正しました.(2019/8/1) アジャイル型なのでお許しください.


computational-sediment-hyd.hatenablog.jp

以前悩んでいた話ですが,勉強したのでまとめておきます. 1次元版ですが,ソースコードをつけてあります.

Jupyter Notebook Viewer

github.com

上手に式変形することにより,TDMAを2回繰り返す形になります.

これで周期境界のsimple法にアプローチできます.

備忘録:proxy環境下でのubuntu18.04のセットアップ

ubuntu用のPCを購入したのでセットアップ時の備忘録です. 毎日やっていればすぐにできるかもしれないですが,久しぶりにやると忘れているのでなんだかんだで半日くらいはかかりますね. GUI付きのubuntu18.04を基本にしてます.

IPの固定

proxy環境だと必須になると思うので. IPアドレスデフォルトゲートウェイDNSネームサーバ等の設定は, GUIでもCLIでもどっちでも同じくらいの手間です. 以下を参考にしました

jyn.jp

プロキシ設定

以前にも書きましたが,以下あたりを参考にしました。

blog.syo-ko.com

qiita.com

windowsからの接続の設定

xrdpによるリモートデスクトップの設定

xrdpは最初から入っているので問題ないですが,細かい設定がいろいろと必要です.

  • 共有の設定:GUIのほうが楽そうです. jikuu.site
  • new_cursorsの無効化,~/.xsessionrcの作成 www.hiroom2.com
  • Authentication Requiredダイアログの回避(上のサイトのやり方だと上手くいきませんでした.)

https://www.serotoninpower.club/archives/207#i-4

これで,リモートデスクトップからxrdpのXorgで接続可能です.細かい設定はまだまだ必要だと思います.

sshの設定

WSLからsshで接続するために設定しています. ubuntusshのサーバーソフトがインストールされていないのでインストールが必要です.

aquarius-train.hatenablog.com

その他

時刻合わせ

proxy環境だと同期がかからない.いろいろな方法があるけどかなり面倒. 同期ではないけど,

ntpdate ntp.nict.jp

でとりあえずはO.K.

kawatama.net

mDNSを活用

windows10 build 1803からmDNSが使えるようになったので,IPアドレスではなく,OSが異なってもPC名.localで接続できます.(なぜか,WSLからは使えない・・・)

kokufu.blogspot.com

IPアドレスからPCのホスト名(NetBIOS名)を見つけたい場合は以下を参照.

qiita.com


インテル Distribution for Pythonってオープンソースなんですね. 普通にanacondaから使えます.

https://jp.xlsoft.com/documents/intel/python/2018/greeneltch_fastpython.pdf#page=14