はじめに
はじめまして！データサイエンティストの渡邊です。株式会社GEOTRA(@GEOTRA )にてGIS（Geographic Information System：地理情報システム）データを用いた分析を行っています。
本記事ではgeopandasでポリゴンを任意の形状に分割する方法を紹介します。
背景
- すでに存在するポリゴンで表現されるエリアを、新たに複数個のエリアに分割して各エリア間の移動人数を調べたい
- そのために隙間なく綺麗に分割されたポリゴンを作成する必要がある
- 一方で、分割後のポリゴンの頂点が多かったり、そもそもの分割数が多くなったりすると、手作業でポリゴンを作成するのは大変
考え方
- クッキー生地を型でくり抜いていくイメージ
- 元のポリゴンがクッキー生地、分割するためのポリゴンが型
- 元のポリゴン$P$を分割するためのポリゴン $P_{1}, ..., P_{n}$ を用意する
- $P_{1}, ..., P_{n}$ は大雑把に作って良い
- 重複はいくらあってもいいが隙間は作らないように注意
- 順番にも気をつける必要がある
- 順番次第では予期せぬ分割のされ方をする
- なんらかの手段でナンバリングしてあげると良い
- $P$と $P_1$の重複部分を抽出し、それを$P$から切り抜く(これを$P'_1$とする)
- 切り抜かれて残ったポリゴン$P_{tmp}$に対して、同様に$P_2$との重複部分を抽出しそれを$P_{tmp}$から切り抜く。これを繰り返す。
- こうして隙間・重複のない分割ポリゴン${P^{\prime}}_{i} (i=1,...,n)$が得られる
実装
ポリゴン作成方法とPythonでの読み込み方法は以下参照
- このポリゴンを3分割することを考える
- まずGoogle Earth上でポリゴン $P_{1}, ..., P_{n}$ を定義しKMLファイルをエクスポート
- ①②③の順で切り抜くために左メニューのポリゴンの順番を整える
- Google Earthだとこのポリゴンの編集が用意
- PythonでKMLファイルを読み込み分割
-
geopandas.GeoSeries.geometry.intersection()- geometry同士の交差部分を計算する
-
geopandas.GeoSeries.geometry.difference()- geometry同士の差分を計算する
-
geopandas.GeoSeries.geometry.intersection()の逆で、交差しない部分を抽出する
-
from pykml import parser # 0.2.0
import geopandas as gpd # 0.10.2
from shapely.geometry import Polygon # 1.8.0
P: gpd.GeoSeries # 元のポリゴン
# 分割用ポリゴンの読み込み
path = 'hogehoge/filename.kml'
with open(path, 'r', encoding='utf-8') as f:
doc = parser.parse(f).getroot().Document
polygons = [] # 分割用ポリゴンのリスト
for p in doc.Placemark:
loc = p.Polygon.outerBoundaryIs.LinearRing.coordinates
loc = str(loc).replace('\n', '').replace('\t', '').split(' ')[:-1]
loc = list(map(lambda x: tuple(map(float, x.split(',')[:2])), loc))
polygons.append(Polygon(loc))
P_tmp = P.copy() # 分割されるポリゴン
divided_polygons = [] # 分割後のポリゴンを格納するリスト
for polygon in polygons:
# P_tmpと分割用ポリゴンの交差部分を抽出
P_intersect = gpd.GeoDataFrame(P_tmp.geometry.intersection(polygon), geometry=0)
divided_polygons.append(P_intersect[0][0]) # リストに追加
# 交差部分を元のポリゴンから取り除く
P_tmp = P_tmp.geometry.difference(P_intersect[0][0])
- kepler.glで確認すると、重複・隙間なく分割できているのが確認できる
from keplergl import KeplerGl # 0.3.2
Map = KeplerGl(height=800)
Map.add_data(data=gpd.GeoDataFrame({'geometry':divided_polygons}), name='分割後のPolygon')
Map
おわりに
最後までお読みくださりありがとうございました。
GEOTRAのQiitaでは、引き続き地理空間情報に関する情報を中心に技術系のテクニックを発信していきます。今後も皆さんのお役に立てるコンテンツを配信できればと思っておりますので、皆様のいいね・フォローをお待ちしています！
また、こちらの記事を読んでGEOTRAに興味を持ってくださった方は、是非Wantedlyから会社紹介記事もご覧ください！GEOTRAは、学生インターンを含め、現在一緒に働く仲間を募集しております。興味をお持ちの方は、Wantedly上の募集ページやストーリー、各種サイトをご確認ください！
企業HP： https://www.geotra.jp/
note： https://note.com/2022geotra/
Wantedly： https://www.wantedly.com/stories/s/geotra_introduction
LinkedIn： https://www.linkedin.com/company/%E6%A0%AA%E5%BC%8F%E4%BC%9A%E7%A4%BEgeotra/