在同一圖層中尋找相鄰的多邊形

備註

If you are using QGIS3, there is a simpler and better way to do this analysis using the aggregate function. See my post Find Neighbor Polygons using Summary Aggregate Function in QGIS

有些時候，你會想要尋找與圖層中的某個多邊形相鄰的所有多邊形。只要弄一點 Python 程式碼，在 QGIS 中就能完成這件事情，甚至做得更多。這裡有一個範例程式，讓你可以尋找所有與指定的多邊形分享邊界的多邊形，並且把它們的名字加入屬性表格中。除此之外，此腳本也把所有鄰接多邊形的某個指定欄位進行加總運算。

內容說明

為了展示腳本如何運作，這邊我們使用全球的國界圖層，然後尋找共享國界的國家。我們也會計算與某個國家相鄰的全部國家的總人口數。

取得資料

我們要使用 Natural Earth 提供的 Admin 0 - 國界 資料集。

下載 Admin 0 - 國界（Countries）shapefile

資料來源 [NATURALEARTH]

取得腳本

下載 neighbors.py 腳本 然後儲存至硬碟中。

操作流程

  1. 選擇 圖層 ‣ 加入向量圖層，載入 ne_10m_admin_0_countries 圖層。

../_images/1182.png

  1. 此腳本使用 2 個欄位執行操作，分別是名稱與用來加總的欄位。使用 識別圖徵 工具點選任何圖徵檢察屬性，在本例中，會使用 NAME 欄位當作名稱，與 POP_EST 欄位當作人口估計值。

../_images/2149.png

  1. 選擇 附加元件 ‣ Python 主控台。

../_images/389.png

  1. 在 Python 主控台 視窗中，按下 顯示編輯器 的按鈕。

../_images/456.png

  1. 在 編輯器 面板中，按下 開啟檔案 鈕，找到剛才下載的 neighbors.py 腳本，按下 開啟。

../_images/549.png

  1. 腳本載入後，或許需要依照圖層屬性，編輯 _NAME_FIELD 和 _SUM_FIELD 的值，不過如果你也是載入 ne_10m_admin_0_countries 圖層的話，使用預設值即可。如果你有編輯過的話，請先按下 編輯器 面板中的 儲存 鈕。接下來，按下 執行腳本 鈕，程式就會開始執行。

../_images/647.png

  1. 當腳本執行完畢後，右鍵點選 ne_10m_admin_0_countries 圖層然後選擇 開啟屬性表格。

../_images/747.png

  1. 你會發現多了 2 個屬性，分別為 NEIGHBORS 和 SUM，這兩個屬性就是剛才的程式新增的。

../_images/846.png

以下放上完整的腳本供各位參考，你也可以依照個人需求編輯調整此檔案。

################################################################################
# Copyright 2014 Ujaval Gandhi
#
#This program is free software; you can redistribute it and/or
#modify it under the terms of the GNU General Public License
#as published by the Free Software Foundation; either version 2
#of the License, or (at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#
################################################################################
from qgis.utils import iface
from PyQt4.QtCore import QVariant

# Replace the values below with values from your layer.
# For example, if your identifier field is called 'XYZ', then change the line
# below to _NAME_FIELD = 'XYZ'
_NAME_FIELD = 'NAME'
# Replace the value below with the field name that you want to sum up.
# For example, if the # field that you want to sum up is called 'VALUES', then
# change the line below to _SUM_FIELD = 'VALUES'
_SUM_FIELD = 'POP_EST'

# Names of the new fields to be added to the layer
_NEW_NEIGHBORS_FIELD = 'NEIGHBORS'
_NEW_SUM_FIELD = 'SUM'

layer = iface.activeLayer()

# Create 2 new fields in the layer that will hold the list of neighbors and sum
# of the chosen field.
layer.startEditing()
layer.dataProvider().addAttributes(
        [QgsField(_NEW_NEIGHBORS_FIELD, QVariant.String),
         QgsField(_NEW_SUM_FIELD, QVariant.Int)])
layer.updateFields()
# Create a dictionary of all features
feature_dict = {f.id(): f for f in layer.getFeatures()}

# Build a spatial index
index = QgsSpatialIndex()
for f in feature_dict.values():
    index.insertFeature(f)

# Loop through all features and find features that touch each feature
for f in feature_dict.values():
    print 'Working on %s' % f[_NAME_FIELD]
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.
    intersecting_ids = index.intersects(geom.boundingBox())
    # Initalize neighbors list and sum
    neighbors = []
    neighbors_sum = 0
    for intersecting_id in intersecting_ids:
        # Look up the feature from the dictionary
        intersecting_f = feature_dict[intersecting_id]

        # For our purpose we consider a feature as 'neighbor' if it touches or
        # intersects a feature. We use the 'disjoint' predicate to satisfy
        # these conditions. So if a feature is not disjoint, it is a neighbor.
        if (f != intersecting_f and
            not intersecting_f.geometry().disjoint(geom)):
            neighbors.append(intersecting_f[_NAME_FIELD])
            neighbors_sum += intersecting_f[_SUM_FIELD]
    f[_NEW_NEIGHBORS_FIELD] = ','.join(neighbors)
    f[_NEW_SUM_FIELD] = neighbors_sum
    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()
print 'Processing complete.'

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)