はじめに
 ポリゴンの作成
 Pythonでのポリゴン読み込み
 おわりに

はじめに

はじめまして！データサイエンティストの渡邊です。株式会社GEOTRA(@GEOTRA )にてGIS（Geographic Information System：地理情報システム）データを用いた分析を行っています。

GISデータを取り扱う際に必ずと言っていいほど出てくるのが、ポリゴンという概念です。ポリゴンは、2次元平面上の閉じた図形のことを指し、複数の頂点を結んで形成されます。例えば、特定のエリアにどれくらい人がいるかを集計する際には、エリアを表現するポリゴンを作成し、そのポリゴン内の座標を抽出することになります。

本記事ではGoogle Earthを用いてポリゴンを作成し、Pythonで読み込む方法を紹介します。

ポリゴンの作成

座標を取得するためにGoogle Earthを用います。利点として、描画したポリゴンを保存できることが挙げられます。利点として以下が挙げられます。

複数のポリゴンを1つの「プロジェクト」としてまとめられる
ポリゴンに名前をつけて保存できる
後からポリゴンの形状を変更したくなったときに、過去に作成したポリゴンの頂点の位置を調整するだけで済む

手順

Pythonでのポリゴン読み込み

上記でポリゴンの情報を含むKMLファイルを得ることができました。KML（Keyhole Markup Language）ファイルとは地理空間情報を記述するためのファイル形式の一種で、XMLベースの形式で地理的な領域、地図上のポイントやライン、領域を定義することができます。

<Placemark id="0EAD842F1130682B2C06">
		<name>東京駅</name>
		<LookAt>
			<longitude>139.7501841111797</longitude>
			<latitude>35.68636344078155</latitude>
			<altitude>37.25960696330495</altitude>
			<heading>0</heading>
			<tilt>0</tilt>
			<gx:fovy>30</gx:fovy>
			<range>5924.133578336332</range>
			<altitudeMode>absolute</altitudeMode>
		</LookAt>
		<styleUrl>#__managed_style_0B7D4AFCD7306829CA29</styleUrl>
		<Polygon>
			<outerBoundaryIs>
				<LinearRing>
					<coordinates>
						139.7677959733631,35.67777874606371,0 139.7691946410166,35.68285571318568,0 139.7668705305633,35.68367674107219,0 139.7652487360985,35.67875838784805,0 139.7677959733631,35.67777874606371,0 
					</coordinates>
				</LinearRing>
			</outerBoundaryIs>
		</Polygon>
	</Placemark>

手順

pykmlを使うのが便利
kmlファイルの処理
- Document.Placemark に各ポリゴンの情報が保存されている
- Document.Placemark.Polygon.outerBoundaryIs.LinearRing.coordinatesに各ポリゴンの頂点の座標が保存されている
- ポリゴンの頂点の座標のx座標, y座標, z座標が空白区切りで保存されている
  - x1,y1,z1 x2,y2,z2 ... xn,yb,zn
  - z座標は基本的に使わないので除く必要がある

from pykml import parser
from shapely.geometry import Polygon

# kmlファイルの読み込み
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: #すべてのPlacemark(=ポリゴン)について処理を行う
    #ポリゴンの頂点の座標を取得
    loc = p.Polygon.outerBoundaryIs.LinearRing.coordinates
    #分割してx,y座標のみ取得, 各座標のtupleのリストを得る
    loc = str(loc).replace('\n', '').replace('\t', '').split(' ')[:-1] # 最後の頂点の座標は最初の頂点と同じなので除く
    loc = list(map(lambda x: tuple(map(float, x.split(',')[:2])), loc)) #floatに変換&z座標を除く
    polygons.append(Polygon(loc)) # リストに追加

可視化

確認のため得られたポリゴンをGISツールを用いて確認していきます。ここではGISツールとしてkepler.glを用います。

まず、得られたポリゴンをgeopandasを通してcsv形式でアウトプットします。

import geopandas as gpd
gpd.GeoDataFrame({'geometry':polygons}).to_csv('polygons.csv')

あとはkepler.glでcsvファイルを読み込むだけです。Google Earthで作成した2つのポリゴンが抽出できているのが確認できます。

kepler.glのライブラリを使えばjupyter notebook上でも確認できます

from keplergl import KeplerGl

Map = KeplerGl()
Map.add_data(data=gpd.GeoDataFrame({'geometry':polygons}))
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/