実行環境
Windows11
python 3.11.7
geopandas 0.14.3
shapely 2.0.2
はじめに
ラインタイプの GIS データを、ポイントタイプの GIS データで切断する方法について書きます。
GIS ソフトの編集機能でもラインを切れますが、数が多い場合は何らかの方法で一括処理したいところです。
QGIS のドキュメントを漁ってみると、線で分割 というツールがありましたが、点で分割というツールは見当たりませんでした。
ArcGIS Pro には、ポイントでラインを分割 というツールがありましたが、最上位の Advanced ライセンスが必要で試すことができませんでした。
また、Esri Japan の EJPyConv に「ラインをポイントで分断」という名前のツールがあり試してみましたが、私の環境では期待の結果となりませんでした。
QGIS のプラグインとか GDAL のベクタ処理とかも軽く探してみましたが、ズバリなツールやコマンドは見つかりませんでした。
以上のような状況につき（需要ありそうなのにどうしてないのだろうかとか、そもそも別の発想で処理すべきものなのだろうか、などと思いつつも）自作することにしました。
python コード
主に geopandas ライブラリを使っています。
splitLineAtPoint() という名前の関数で、引数は以下の３つです。
| 引数 | 説明 |
|---|---|
| gdf_line | 切断したいラインデータ |
| gdf_pt | 切断する地点を示すポイントデータ |
| tolerance | （ライン-ポイント間の）許容誤差 |
考え方は後述します。コードはこちらです。
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
def splitLineAtPoint(
gdf_line: gpd.GeoDataFrame,
gdf_pt: gpd.GeoDataFrame,
tolerance: float
) -> gpd.GeoDataFrame:
# ポイントから許容誤差でバッファ
gdf_pt['_buffer'] = gdf_pt.buffer(distance=tolerance)
gdf_pt = gdf_pt.explode(index_parts=False)
gdf_pt = gdf_pt.reset_index(drop=False)
# ラインをバッファでくり抜き
gdf_line = gpd.overlay(
df1=gdf_line,
df2=gdf_pt[['_buffer']].set_geometry('_buffer'),
how='symmetric_difference',
keep_geom_type=True
)
gdf_line = gdf_line.explode(index_parts=False)
# apply() 用の関数
def snap_edge(
row: gpd.GeoSeries,
s_pt: gpd.GeoSeries
):
xy_new = list(row['geometry'].coords)
# ラインの終点[-1]と始点[0]を処理
for i in [-1, 0]:
pt_edge = Point(row['geometry'].coords[i])
s_distance = s_pt.distance(pt_edge)
# ポイントから許容誤差範囲内の端点に対して
# 端点の座標をポイントの座標に差し替え
if s_distance.min() < tolerance:
pt_new = s_pt[s_distance.idxmin()]
xy_new[i] = list(pt_new.coords)[0]
return LineString(xy_new)
# 端点をポイントにスナップ
gdf_line['geometry'] = gdf_line.apply(
snap_edge,
axis=1,
s_pt=gdf_pt['geometry']
)
return gdf_line
考え方
splitLineAtPoint() をざっくり説明します。
ラインをポイントのバッファでくり抜いてから、ラインの切り口同士を接合するという処理をしています。
くり抜きの処理は、geopandas の overlay() 関数を how='symmetric_difference' として行っています。
df2 引数でポイントのバッファを指定していますが、このバッファの半径が tolerance です。
ライン上ぴったりにないポイントを切断点として利用するための許容誤差として機能します。
切り口の接合の処理は snap_edge() にもたせていて、切り取ったラインに対して apply() しています。
途中で explode() したり reset_index() したりしているのは、中間データをシンプルなシングルジオメトリに直すためです。
サンプルデータ
サンプルで試してみます。
今回は、国土数値情報の 高速道路時系列データ の2022年（令和4年）のデータを使います。
このデータには、高速道路の路線を示すラインデータ と 接合部（IC や JCT など）を示すポイントデータ が含まれています。
-
圧縮ファイル名
N06-22_GML.zip -
使用するGISデータファイル名
N06-22_HighwaySection.*
N06-22_Joint.*
QGIS で静岡市あたりを表示してみます。
凡例に示している通り、黒い線が路線データ、赤丸が接合部データです。
赤丸が大きく見えますが、ポイントを大きく表示しているだけで、バッファをとっているわけではないです。
路線データには、切れ目がわかりやすいように始終点を白丸にしています。
白丸と赤丸が重なっている地点は、すでにポイントでラインが切断されていることを意味します。
したがって、すべての赤丸のうえに白丸が表示されれば、ポイントでラインを切断できたことになります。
前処理
データは地理座標系なので、あらかじめ投影座標系に変換しておきます。
許容誤差や円を扱ううえで、緯度経度ではなくメートルで処理したいためです。
どのような方法で変換してもよいですが、今回は GDAL で以下の通り投影変換しました。
> ogr2ogr -lco ENCODING=SJIS -s_srs EPSG:6668 -t_srs EPSG:6677 N06-22_HighwaySection_projected.shp N06-22_HighwaySection.shp
> ogr2ogr -lco ENCODING=SJIS -s_srs EPSG:6668 -t_srs EPSG:6677 N06-22_Joint_projected.shp N06-22_Joint.shp
日本測地系2011 を 平面直角座標系第9系 に投影しています。
本処理
処理を実行します。
許容誤差は 10m にしました。
# 入力ファイル
shp_line = r'N06-22_GML\N06-22_HighwaySection_projected.shp'
shp_pt = r'N06-22_GML\N06-22_Joint_projected.shp'
# 許容誤差
tolerance = 10
# 出力ファイル
shp_split_line = 'split_line.shp'
# 読込
gdf_line = gpd.read_file(shp_line)
gdf_pt = gpd.read_file(shp_pt)
# 処理
gdf_split_line = splitLineAtPoint(
gdf_line=gdf_line,
gdf_pt=gdf_pt,
tolerance=tolerance
)
# 出力
gdf_split_line.to_file(shp_split_line)
出力データの確認
出力結果を QGIS で確認します。
すべての赤丸が白丸と重なっているので、ポイントでラインを切断できているようです。
注意点
-
ラインとポイントの座標系を同じ投影座標系で揃えましたが、異なる座標系や地理座標系を扱う場合は注意が必要かと思います
-
切断に使うポイントが密集する場合や、許容誤差が大きすぎる場合は、期待の結果とならないかもしれません
-
許容誤差を 0 にしてしまうとバッファができないため、ラインをまったく切断しなくなります
-
ポイントがライン上ぴったりではない場合、ラインの切り口がライン沿いのポイントにスナップするため、元のラインに対してズレが生じます
これらの課題に対する対処法はいくらでもあるかとは思いますが、主題から逸れるため取りあげていません。
おわりに
ラインタイプの GIS データを、ポイントタイプの GIS データで切断する方法について書きました。
ポリゴンでラインをくり抜いてから接合し直す、という段階を踏む方法を採りましたが、ジオメトリの M値 を使って切断する方法でもいけそうな気がします。
属性情報については触れませんでしたが、今回掲載のコードだと、出力データに入力のラインデータの属性が受け継がれることになります。（文字化けする場合は GeoDataFrame の入出力で文字コードの指定が必要です）