在同一图层中寻找相邻的多边形

有些时候,你会想要寻找与图层中的某个多边形相邻的所有多边形。只要弄一点 Python 程序代码,在 QGIS 中就能完成这件事情,甚至做得更多。这里有一个范例程序,让你可以寻找所有与指定的多边形分享边界的多边形,并且把它们的名字加入属性表格中。除此之外,此脚本也把所有邻接多边形的某个指定栏位进行加总运算。

内容说明

为了展示脚本如何运作,这边我们使用全球的国界图层,然后寻找共享国界的国家。我们也会计算与某个国家相邻的全部国家的总人口数。

取得资料

我们要使用 Natural Earth 提供的 Admin 0 - Countries 资料集。

下载 Admin 0 - countries shapefile..

资料来源 [NATURALEARTH]

取得脚本

下载 neighbors.py script 然后储存至硬碟中。

操作流程

  1. 选择:menuselection:Layer –> Add Vector Layer,载入``ne_10m_admin_0_countries`` 图层。
../_images/140.png
  1. 此脚本使用 2 个栏位执行操作,分别是名称与用来加总的栏位。 使用 Identify 工具点选任何图征检察属性,在本例中,会使用**NAME** 栏位当作名称,与 POP_EST 栏位当作人口估计值。
../_images/225.png
  1. 选择 Plugins ‣ Python Console
../_images/316.png
  1. Python Console 视窗中,按下:guilabel:Show Editor 的按钮。
../_images/411.png
  1. Editor 面板中,按下 Open file 钮,找到刚才下载的 neighbors.py 脚本,按下 Open
../_images/511.png
  1. 脚本载入后,或许需要依照图层属性,编辑 _NAME_FIELD 和``_SUM_FIELD`` 的值, 不过如果你也是载入 ne_10m_admin_0_countries 图层的话, 使用预设值即可。如果你有编辑过的话,请先按下:guilabel:Editor 面板中的 Save 钮。 接下来,按下 Run script 钮,程序就会开始执行。
../_images/610.png
  1. 当脚本执行完毕后,右键点选 ne_10m_admin_0_countries 图层然后选择:guilabel:Open Attribute Table
../_images/710.png
  1. 你会发现多了 2 个属性,分别为 NEIGHBORSSUM,这两个属性就是刚才的程序新增的。
../_images/810.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.'