趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

河川断面データを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