Нахождение соседних полигонов в слое

Примечание

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 это и многое другое. Вот пример сценария, который вы можете использовать, чтобы найти все полигоны, которые имеют общую границу с каждым из полигонов в слое, а также добавить их названия в таблицу атрибутов. В качестве дополнительного бонуса, сценарий также суммирует атрибут на ваш выбор по всем соседним полигонам.

Обзор задачи

Чтобы продемонстрировать, как работает скрипт, мы будем использовать слой полигонов стран и найдем страны, которые имеют общие границы. Мы также хотим вычислить общее население соседних стран для каждой страны.

Получение данных

Мы будем использовать набор данных Admin 0 - Countries из Natural Earth.

Загрузите shape-файл Admin 0 - countries..

Источник данных: [NATURALEARTH]

Получение сценария

Загрузите сценарий neighbors.py и сохраните его на диск.

Методика

  1. Загрузите слой ne_10m_admin_0_countries перейдя к меню Layer ‣ Add Vector Layer.

../_images/1176.png
  1. Скрипт использует 2 поля для выполнения действия. Название поля и поле, которое вы хотите просуммировать. Используйте инструмент Identify, чтобы нажать на любом объекте и просмотреть его атрибуты. В данном случае название поля NAME, и мы хотим просуммировать оценку численности народонаселения из поля POP_EST.

../_images/2144.png
  1. Перейдите к меню Plugins ‣ Python Console.

../_images/380.png
  1. В окне Python Console нажмите на кнопку Show Editor.

../_images/448.png
  1. В панели Editor нажмите на кнопку Open file и перейдите к загруженному сценарию neighbors.py, после чего нажмите Open.

../_images/546.png
  1. После того, как скрипт загружен, вы можете изменить значения _NAME_FIELD и _SUM_FIELD в соответствии с атрибутами вашего собственного слоя. Если вы работаете со слоем ne_10m_admin_0_countries, вы можете оставить всё, как есть. Нажмите: guilabel: save кнопку в: guilabel:` панели Editor`, если вы сделали какие-либо изменения. Теперь нажмите: guilabel: кнопку Run script, для выполнения сценария.

../_images/644.png
  1. После того, как скрипт закончит работу, щелкните правой кнопкой мыши на слое ne_10m_admin_0_countries и выберите Open Атрибут Table.

../_images/744.png
  1. Вы заметите,2 новых атрибута под названием NEIGHBORS и SUM. Они были добавлены сценарием.

../_images/843.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 report any issues with this tutorial, please comment below. (requires GitHub account)