【Python】QGISで、細々としたポリゴンを大きいポリゴンに融合してみた。

はじめに

QGISで小さな面積のポリゴンを大きな面積のポリゴンに融合(dissolve)する方法をPythonで実装した。
もっと良い方法や、なにかプラグインがあれば教えていただけると幸いです。

できること

以下の左図のような細々としたポリゴンを、右図のように大きなポリゴンに融合しまとめることができる。
　　　

ソースコード

それぞれのポリゴンの面積を計算し、その面積ごとに順位を決定し、設定した閾値より高い順位のポリゴン(つまり面積の小さいポリゴン)を閾値以下の順位のポリゴンに融合する。融合先ポリゴンはポリゴンの重心間距離が最も小さいポリゴンとしています。
QGISのPythonコンソールで利用することを想定してます。

import qgis
from qgis.core import *
import processing

# レイヤのフィールドの列番号を返す(指定したレイヤが存在しなければ作成して値を返す)
def getFieldsNumber(layer, field_name):
    ret = layer.fields().lookupField(field_name)
    if ret == -1:
        layer_provider=layer.dataProvider()
        layer_provider.addAttributes([QgsField(field_name, QVariant.Double)])
        layer.updateFields()
        ret = layer.fields().lookupField(field_name)
    return ret

# 何種類のポリゴンに融合するか指定
THRESHOLD = 5

aLayer = QgsProject.instance().mapLayersByName('融合前レイヤ名')[0] #レイヤ名を指定して取得
#aLayer = qgis.utils.iface.activeLayer()  # active layerで実行したい人はこちらを使用

fni_area    = getFieldsNumber(aLayer, 'Area')
fni_rank    = getFieldsNumber(aLayer, 'AreaRank')
fni_filter  = getFieldsNumber(aLayer, 'AreaFilter')
fni_x       = getFieldsNumber(aLayer, 'Centroid_x')
fni_y       = getFieldsNumber(aLayer, 'Centroid_y')
fni_nearest = getFieldsNumber(aLayer, 'Nearest')

# 各ポリゴン毎の面積、行番号、重心座標取得
aList= [(feat.geometry().area(), feat.id(), feat.geometry().centroid()) for feat in aLayer.getFeatures()]
aList.sort(reverse = True)  # ソートして面積を昇順に並び替え

# 編集モードを開始
aLayer.startEditing()

# 面積の順位付けと融合するポリゴンの選定
rank = 1
for feat in aList:
    id = feat[1] # 行番号
    aLayer.changeAttributeValue(id, fni_rank, rank)
    if rank <= THRESHOLD:  # 閾値以下であれば1を格納(融合しない)
        aLayer.changeAttributeValue(id, fni_filter, 1)
    else:                  # 融合するポリゴンには0を格納
        aLayer.changeAttributeValue(id, fni_filter, 0)
    aLayer.changeAttributeValue(id, fni_area, feat[0])
    aLayer.changeAttributeValue(id, fni_x, feat[2].asPoint().x())
    aLayer.changeAttributeValue(id, fni_y, feat[2].asPoint().y())
    rank += 1

# 融合先ポリゴンの重心座標を格納(x, y, rank)
isFilter_centroid_lst = []
for i, feat in enumerate(aLayer.getFeatures()):
    isFilter = feat.attributes()[fni_filter]
    if isFilter == 1:
        isFilter_centroid_lst.append((feat.attributes()[fni_x], feat.attributes()[fni_y], feat.attributes()[fni_rank]))
        aLayer.changeAttributeValue(i, fni_nearest, feat.attributes()[fni_rank])

# 融合するポリゴンの最近傍の融合先ポリゴンを計算
for i, feat in enumerate(aLayer.getFeatures()):
    isFilter = feat.attributes()[fni_filter]
    if isFilter == 0:
        x_feat = feat.attributes()[fni_x] #重心座標取得
        y_feat = feat.attributes()[fni_y]
        nearest_id = -1
        nearest_distance = 1000000000000000000
        for lst in isFilter_centroid_lst: #融合先ポリゴンの重心座標で最も近いものを算出
            x, y, id = lst
            distance = (x_feat - x)**2 + (y_feat - y)**2
            if distance < nearest_distance:
                nearest_id = id
                nearest_distance = distance
        aLayer.changeAttributeValue(i, fni_nearest, nearest_id)

#編集モード終了
aLayer.commitChanges()  # save changed and stop editing

#融合処理
outfn = r".\dissolved_{}areas.shp".format(THRESHOLD)
processing.run("gdal:dissolve", {'INPUT':aLayer,'FIELD':"Nearest", 'OUTPUT':outfn}) # Nearestで融合
iface.addVectorLayer(outfn, "", "ogr")

参考文献