在同一圖層中尋找相鄰的多邊形¶
備註
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 腳本
然後儲存至硬碟中。
操作流程¶
選擇
,載入ne_10m_admin_0_countries
圖層。

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

選擇
。

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

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

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

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

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

以下放上完整的腳本供各位參考,你也可以依照個人需求編輯調整此檔案。
################################################################################
# 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.'