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