python の shapely を使って複数のポリゴンを可能な限り結合する時のメモ

準備

インストールする。

pip install shapely

インポートする。

from shapely.geometry import Polygon
from shapely.ops import unary_union

前提

(django の話になってしまうが)モデル Area に multipolygon フィールドがあり、型が models.MultiPolygonField() である事。

データは下記のような値であること。

MULTIPOLYGON(((139.770167304785 35.7297682113412,139.794543220313 35.6851633103714,139.827502204688 35.7141593351991,139.807932807715 35.7425875003111,139.770167304785 35.7297682113412)),((139.799349738867 35.7253088444921,139.748537971289 35.7241939637722,139.818575813086 35.6807014476555,139.850504829199 35.7216854251056,139.799349738867 35.7253088444921)))

コード

メイン関数

def get_union_polygons(request):
    """
    可能な限りエリアのマルチポリゴンを結合して返す。
    """
    areas = Area.objects.all()
    # 結合できる shaply_polygon に変換する。
    shaply_polygons = []
    for area in areas:
        multipolygon = area.multipolygon
        for polygon in multipolygon:
            tuple_polygon = polygon.tuple[0]
            list_polygon = [coordinate for coordinate in tuple_polygon]
            shaply_polygons.append(Polygon(list_polygon))
    # 可能な限りポリゴンを結合する。
    merged_polygons = []
    skip_indexes = []
    for index, shaply_polygon in enumerate(shaply_polygons):
        if index not in skip_indexes:
            # ポリゴンを結合する。
            merged_polygon, merged_indexes = merge_polygons(shaply_polygon, shaply_polygons, index)
            merged_polygons.append(merged_polygon)
            # 一度でも結合したポリゴンは、再結合は不要なので飛ばす。
            skip_indexes += merged_indexes
    # shaply から tuple に変換してエンコードする。
    encode_polygons = []
    for merged_polygon in merged_polygons:
        tuple_polygon = shaply_to_tuple(merged_polygon)
        encode_polygon = polyline.encode(tuple_polygon, geojson=True)
        encode_polygons.append({
            "encodePolygon": [encode_polygon],
        })
    return HttpResponse(json.dumps(encode_polygons))

結合する関数

def merge_polygons(target_polygon, polygons, index):
    """
    第一引数のポリゴンに可能な限り、第二引数のポリゴンを結合して返す。
    """
    # 一応。
    merged_indexes = [index]
    # 自分自身との結合は不要なので次のポリゴンから比較するため +1 する。
    start_index = index + 1
    # 指定できる最後のインデックスは、配列の要素数 -1 番目まで。
    last_index = len(polygons) - 1
    i = start_index
    while i <= last_index:
        # 一度結合した場合再結合は不要。
        if i in merged_indexes:
            i += 1
            continue
        # 結合した場合、結合できるポリゴンが出てくる可能性があるのでもう一度最初から結合できるかを検証する。
        if target_polygon.intersects(polygons[i]):
            merged_indexes.append(i)
            target_polygon = unary_union([target_polygon, polygons[i]])
            i = start_index
            continue
        i += 1
    return target_polygon, merged_indexes

変換する関数

def shaply_to_tuple(shaply_polygon):
    """
    shaply を tuple に変換する。
        【変換前】
            POLYGON ((139.8122815591052 35.69123599843862, 139.7913388711169 35.65776849358565, 139.8356275063708 35.65581578929326, 139.8122815591052 35.69123599843862))
        【変換後】
            ((139.8122815591052, 35.69123599843862), (139.7913388711169, 35.65776849358565), (139.8356275063708, 35.65581578929326), (139.8122815591052, 35.69123599843862))
    """
    shaply_polygon.wkt
    shaply_polygon = shaply_polygon.wkt
    shaply_polygon = shaply_polygon.replace("((", "")
    shaply_polygon = shaply_polygon.replace("))", "")
    shaply_polygon = shaply_polygon.replace("POLYGON ", "")
    shaply_polygon = shaply_polygon.split(", ")
    coordinates = []
    for x in shaply_polygon:
        latlng = x.split(" ")
        coordinate = (float(latlng[0]), float(latlng[1]))
        coordinates.append(coordinate)
    return tuple(coordinates)