Encontrar Polígonos Vecinos en una capa

Nota

Si estás usando QGIS3, hay una manera más simple y mejor de hacer este análisis usando la función agregar. Vea my publicación Find Neighbor Polygons using Summary Aggregate Function in QGIS

Hay algunos casos de uso donde quieres encontrar todos los polígonos que son vecinos de cada uno de los polígonos en una capa. Con un pequeño script python, puedes alcanzar esto y mucho más en QGIS. Aquí tiene un script ejemplo que puede usar para encontrar todos los polígonos que comparten borde con cada uno de los polígonos en una capa y también agregar sus nombres a la tabla de atributos. Como un bono adicional, el script también acumula la suma de un atributo de su elección de todos los polígonos vecinos.

Vista general de la tarea

Para demostrar como funciona el script, usaremos una capa de polígonos de países y encontraremos países que comparte el límite. También queremos calcular la población total de los vecinos del país.

Obtener los datos

Usaremos el conjunto de datos Admin 0 - Countries de Natural Earth.

Descargar el archivo shape Admin 0 - countries..

Fuente de Datos [NATURALEARTH]

Obtener el script

Descargue el neighbors.py script y guárdelo en su disco.

Procedmiento

  1. Cargue la capa ne_10m_admin_0_countries yendo a Capa ‣ Agregar Capa Vectorial.

../_images/1174.png
  1. El script usa 2 campos para realizar la acción. Un campo nombre y un campo que quieres que se sume. Usa la herramienta Identificar para hacer clic en cualquier elemento y examinar los atributos. En este caso el campo nombre es NAME y queremos sumar los estimados de población a partir del campo POP_EST.

../_images/2142.png
  1. Vaya a Complementos ‣ Consola de Python.

../_images/379.png
  1. En la ventan Python Console, clic al botón Mostrar Editor

../_images/447.png
  1. En el panel Editor, clic al botón Abrir archivo y navegue al script descargado neighbors.py y clic Abrir.

../_images/545.png
  1. Una vez que el script está cargado, puedes querer cambiar los valores _NAME_FIELD y _SUM_FIELD para que coincidan con los atributos de tu propia capa. Si estás trabajando con la capa ne_10m_admin_0_countries, puedes dejarlos tal como están. Clic al botón Guardar en el panel Editor si has hecho algún cambio. Ahora clic al botón Ejecutar script para ejecutar el script.

../_images/643.png
  1. Una vez que el script finalizar, clic-derecho a la capa ne_10m_admin_0_countries y selecciona Abrir Tabla de Atributos.

../_images/743.png
  1. Notarás 2 nuevos atributos llamados NEIGHBORS y SUM. Estos fueron agregados por el script.

../_images/842.png

Abajo está el script completo de referencia. Puedes modificarlo para que se ajuste a tus necesidades.

################################################################################
# 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.'
comments powered by Disqus