pythonで2次元の領域を扱うのにshapelyが便利だった

背景

装置の座標が地図上の区画に入っているかどうかの判別を行う機会があったのでpythonでそのような処理が出来るライブラリを探したところ、図形を扱うライブラリshapelyで期待通りの処理を組めました。
備忘録として使い方を記事に残します。

使ったもの

python3の実行環境

ライブラリをインストール
pipコマンドでインストールできます。

pip install shapely

ubuntuだとapt版もありますが、それだとバージョンが古くて最新の書き方と異なり参考にした書き方ではエラーになったのでpip版を使うのが良いです。

区画をポリゴンで定義し、点が中に居るか判別

Polygonのcontains関数で判別できます。

from shapely import Polygon, Point

arrPoligon = [
    Polygon([(36.17615365206663, 139.29274439262977),
             (36.176069213304274, 139.29292611226978),
             (36.17615744098333, 139.29299383803968),
             (36.17625487020672, 139.2928342466215)]),
    Polygon([(36.17610006593953, 139.29270885336436),
             (36.17600209524877, 139.29289124355657),
             (36.17592306907745, 139.29284095214328),
             (36.17603565428372, 139.29261698771606)]),
]

for point in arrPoint:
    for polygon in arrPoligon:
        if polygon.contains(point): # <<<< ここです
            print("%s contains %s" % (polygon, point))

判別できました。

POLYGON ((36.17615365206663 139.29274439262977, 36.176069213304274 139.29292611226978, 36.17615744098333 139.29299383803968, 36.17625487020672 139.2928342466215, 36.17615365206663 139.29274439262977)) contains POINT (36.17614534102644 139.29286477985536)
POLYGON ((36.17610006593953 139.29270885336436, 36.17600209524877 139.29289124355657, 36.17592306907745 139.29284095214328, 36.17603565428372 139.29261698771606, 36.17610006593953 139.29270885336436)) contains POINT (36.17601110498835 139.29275883261136)

ちなみに、ポリゴンが意味する数値は、記事を書いている時点で埼玉県深谷市にある株式会社レグミンの駐車場の北と南の区画です。
google mapで右クリックすると緯度経度を取得できるので、拡大率を最大にして4点で1つの区画を作りました。

可視化

shapelyの可視化方法は公式には案内されてないようなので用途に応じた画像描画ライブラリを使うと可視化できます。
自分が探した範囲ではmatplotlibを利用する例が多かったので、それを利用して可視化しました。

from shapely import Polygon, Point
import matplotlib.pyplot as plt

arrPoligon = [
    Polygon([(36.17615365206663, 139.29274439262977),
             (36.176069213304274, 139.29292611226978),
             (36.17615744098333, 139.29299383803968),
             (36.17625487020672, 139.2928342466215)]),
    Polygon([(36.17610006593953, 139.29270885336436),
             (36.17600209524877, 139.29289124355657),
             (36.17592306907745, 139.29284095214328),
             (36.17603565428372, 139.29261698771606)]),
]

arrPoint = [
    Point(36.17614534102644, 139.29286477985536),
    Point(36.17608958979661, 139.2928252172769),
    Point(36.17601110498835, 139.29275883261136),
]

for polygon in arrPoligon:
    plt.fill(*polygon.exterior.xy, color="green")
    # plt.plot(*polygon.exterior.xy, color="green")

for point in arrPoint:
    isPointInSomePolygon = False
    for polygon in arrPoligon:
        if polygon.contains(point):
            print("%s contains %s" % (polygon, point))
            isPointInSomePolygon = True
            break
    plt.scatter(*point.xy, color= "orange" if isPointInSomePolygon else "black")

plt.plot()
plt.show()

ポリゴンの領域を緑色で塗り、ポリゴンの中に居ると判別された点はオレンジ、そうでなければ黒で描画しました。

ポリゴンを塗りつぶしではなく線画で良い場合はfill関数ではなくplot関数を使います。

for polygon in arrPoligon:
    plt.plot(*polygon.exterior.xy, color="green")

描画の設定はpyplotの説明ページが参考になります。
matplotlib.pyplot.plot

おわり

当初は大変になりそうだと思っていた地図上の区画判別処理がshapelyを使うことで楽々判別でき、matplotlibを利用したら描画までできました。
便利でありがたいです。

参考

shapelyの存在を知ったstack overflowの投稿です。
What's the fastest way of checking if a point is inside a polygon in python

shapelyをmatplotlibで描画しているgithubのissueやstack overflowの投稿です。
Spurious splits when intersecting MultiLineString Z with Polygon
How do I plot Shapely polygons and objects using Matplotlib?

描画ライブラリmatportlibの関数の引数説明ページです。
matplotlib.pyplot.plot