実行環境
Windows10
Python 3.9.13
pyproj 3.3.1
rdp 0.8
pyshp 2.3.1
はじめに
rdp というライブラリで、シェープファイルのポリゴンの頂点を間引きます。
Ramer-Douglas-Peucker というアルゴリズムが使われていて、ライブラリ名はそのイニシャルからきています。
シェープファイルを扱うのであれば shapely ライブラリにも simplify というメソッドがあり、rdp と同様のアルゴリズムで間引き処理ができます。
ただ今回は、GIS データ以外にも応用が利くように rdp を取り上げたいと思います。
Ramer-Douglas-Peucker アルゴリズム
大層な名前でとっつきにくいですが、ライブラリのドキュメント に掲載のアニメーションで大体分かります。
原始的で直観的で処理時間のかかるやり方です。
無理やり要約するなら、「間引きをして生じる線分」と「もともとの線の最遠の頂点」の間の距離が許容値 ε 以下なら間引く、これの繰り返しです。
間引き後の座標は間引き前の座標のサブセットになります。
ウィキペディアだと英語版になりますが、詳しい説明がありました。
rdp 関数の構文
座標を受け取って、間引いた座標を返します。
rdp.rdp(
M,
epsilon=0,
dist=pldist,
algo='iter',
return_mask=False
)
rdp 関数の引数
M
頂点の座標を python リストまたは np.array で指定します。
epsilon
アルゴリズムで規定するところのイプシロン ε を指定します。
この値は距離なので、平面直角座標では必然的に m 単位です。
値が大きいほど間引き処理が発生しやすくなります。
デフォルトが 0 で、一直線のラインの中に複数の頂点がある場合だけ間引くことになります。
dist
距離を計測する方法を指定したい場合は、関数で指定します。
デフォルトの pldist は、ライブラリの中で定義されている距離を計算する関数です。
座標値通り幾何学的に計算するだけなら、指定の必要はありません。
algo
'iter' または 'recursive' を指定します。
関数の中身で while を使う方法か再帰関数を使う方法かを選べるということのようです。
今回の例では、出力データに違いはありませんでした。
（処理速度が違うとかあるのかもしれないです。）
return_mask
この引数を True に指定すると、座標ではなく bool が返るようになります。
返ってくる bool は、残る座標が True、間引く座標が False です。
テストケース
国土数値情報の行政区域データ から新宿区の行政区域ポリゴンを抜き出して、間引き処理します。
感度分析のために ε を 5 通り試してみます。
国土数値情報は地理座標系なので、座標を投影座標系に変換して距離を m で扱うようにしています。
コード
import pyproj
from pyproj.enums import WktVersion
import rdp
import shapefile
fp_zip = r'shp\N03-20220101_13_GML.zip' # 東京都の行政区域ポリゴン
fp_shp = r'shp\polygon_rdp.shp' # 出力ファイル
fp_prj = r'shp\polygon_rdp.prj' # 出力ファイル
sf_r = shapefile.Reader(fp_zip)
sf_w = shapefile.Writer(fp_shp, shapeType=5) # 5: ポリゴン
poly = sf_r.shape(12) # レコード12が新宿区
coords = [coord for coord in poly.points] # 頂点座標
print(f'vertex={len(coords): >4}') # 頂点数
trans = pyproj.Transformer.from_crs(
6668, # 世界測地系2011（緯度経度）
6677, # 世界測地系2011の平面直角座標系第9系（m）
always_xy=True
)
coords_trans = trans.transform(*list(zip(*coords))) # 座標変換
list_epsilon = [0, 1, 5, 20, 100] # ε リスト
sf_w.field('epsilon', 'N') # フィールド追加
for epsilon in list_epsilon:
# 間引き
coords_rdp = rdp.rdp(
list(zip(*coords_trans)),
epsilon=epsilon
)
sf_w.poly([coords_rdp]) # ジオメトリ書き込み
sf_w.record(epsilon) # 属性書き込み
print(f'epsilon={epsilon: >3}, vertex={len(coords_rdp): >4}')
sf_w.close()
# 座標系ファイル作成
with open(fp_prj, 'w') as f:
wkt = pyproj.CRS.from_epsg(6677).to_wkt(WktVersion.WKT1_GDAL)
f.write(wkt)
結果
標準出力で、もともとの頂点数と ε ごとの間引き後頂点数を表示しています。
vertex=1277
epsilon= 0, vertex=1277
epsilon= 1, vertex= 501
epsilon= 5, vertex= 301
epsilon= 20, vertex= 170
epsilon=100, vertex= 51
あくまでも新宿区のポリゴンの話ですが、1m 許容すれば、頂点数が 39% (501/1277) になります。
100m まで妥協すれば、わずか 4% (51/1277) です。
次にポリゴンの形状を比較します。
背景に地理院タイルを使って、QGISで表示しています。
黒線がもともとの行政区域の境界です。
区全体のスケールで見た場合、ε=100（赤）では間引きが目立ちますが、ε=20（オレンジ）では処理前とほぼ一致しているように見えます。
新宿駅周辺を見るスケールだと、ε=20（オレンジ）の間引きが確認できます。
最後に
python ライブラリ rdp で、シェープファイルのポリゴンの頂点を間引きました。
いまどき、ベクタデータの頂点程度で間引きが必要になることはそれほどないのかもしれないですが。
（点群データならともかく）
なお今回は Ramer-Douglas-Peucker を使いましたが、別の間引きのアルゴリズムとかカーブフィッティングとかでデータ量を削減する方法もあるので、データや目的に応じて適切な方法を使いましょう。