2

I have a shapefile with municipality polygons. I would like to count the number of neighbours for each polygon and write it to the attribute table. I have no succeeded yet in doing this, even though I tried to run this script: http://www.qgistutorials.com/de/docs/find_neighbor_polygons.html

2

The script from the linked website writes not the number of neighboring polygons but the summed up value for an attribute from the layer (e.g. the total population for all neighboring countries) to the attribute table.

The easiest way to achieve your goal is probably to change line 74 in the script from neighbors_sum += intersecting_f[_SUM_FIELD] to neighbors_sum += 1

As a result, you will get a column SUM with the number of neighboring polygons.

In case this is not working, check if the polygons in your data meet or overlap, because otherwise the script will not work like expected.

1

Thanks for pointing this out. Another problem was the line

f[_NEW_NEIGHBORS_FIELD] = ','.join(neighbors)

which only works for string fields. I wanted to count on a integer field. The comma-separated string can be built with this line instead of the one above:

f[_NEW_NEIGHBORS_FIELD] = ",".join(map(str, neighbors))

Also, there's another way to count the number of neighbours, like this:

f[_NEW_SUM_FIELD] = len(neighbors)

It just counts the number of items in the neighbours list with len()

The full code is then:

from qgis.utils import iface
from PyQt4.QtCore import QVariant


_NAME_FIELD = 'BFS_NUMMER'


_SUM_FIELD = 'BFS_NUMMER'

# 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]
    # if field is Strings:
    #f[_NEW_NEIGHBORS_FIELD] = ','.join(neighbors)
    # if field is Integer:
    f[_NEW_NEIGHBORS_FIELD] = ",".join(map(str, neighbors))
    # count sum values
    #f[_NEW_SUM_FIELD] = neighbors_sum
    # count list and get number of neighbours
    f[_NEW_SUM_FIELD] = len(neighbors)
    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()
print 'Processing complete.'
  • You should edit the tutorial with the revised script. Or propose it for updating. – Victor Nov 24 '16 at 10:05
1

For QGIS3, there is a simpler way to do this using the aggregate expression. Add a new field with the following expression. Assuming your layer name is 'polygon'

aggregate(
 layer:= 'polygon',
 aggregate:='count',
 expression:=$id,
 filter:=touches($geometry, geometry(@parent))
 )

See more at https://spatialthoughts.com/2019/05/23/neighbor-polygons-aggregate-qgis/

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.