How to count duplicate points in QGIS using Python

In this task we will use python scripting to count the number of duplicate geometries from a point vector layer and store the result in a memory layer.

Starting, from something like this (overlapping points are not visible!):

input

we will see how it’s easy (and fast) to do the task using a QgsSpatialIndex.

The first step is to count the number of duplicate geometries from the input point vector. Then, a new memory layer is generated, containing the same initial attributes plus two additional fields:

  • the “COUNT” field, which stores the result of the counting operation (i.e. the total number of point features which lie on the same position);
  • the “DUPLICATES” field, which store the number of duplicate features (i.e. the result of the counting operation decreased by one because the using of the Spatial Index considers the current point itself).

We can use the following code (if you don’t know how to create a new script from Processing, please see here):

##Point_Vector=vector
##COUNT_field_name=string COUNT
##DUPLICATES_field_name=string DUPLICATES

from qgis.core import *
from qgis.PyQt.QtCore import QVariant

layer = processing.getObject(Point_Vector)
crs = layer.crs().toWkt()

# Create the output layer
outLayer = QgsVectorLayer('Point?crs='+ crs, 'outLayer' , 'memory')
prov = outLayer.dataProvider()
fields = layer.pendingFields() # Fields from the input layer
fields.append(QgsField(COUNT_field_name, QVariant.Int, '', 10, 0))
fields.append(QgsField(DUPLICATES_field_name, QVariant.Int, '', 10, 0))
prov.addAttributes(fields) # Add input layer fields to the outLayer
outLayer.updateFields()

points = [feat for feat in layer.getFeatures()]

index = QgsSpatialIndex() # this spatial index contains all the features of the point layer
for point in points:
    index.insertFeature(point)

for point in points:
    inAttr = point.attributes() # Input attributes
    tmp_geom = point.geometry() #Input geometry
    idsList = index.intersects(tmp_geom.boundingBox())
    count = len(idsList)
    outGeom = QgsFeature()
    inAttr.append(count)
    duplicates = count - 1
    inAttr.append(duplicates)
    outGeom.setAttributes(inAttr) # Output attributes
    outGeom.setGeometry(tmp_geom) # Output geometry
    prov.addFeatures([outGeom]) # Output feature
# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

You can save the new vector layer by right clicking on its name in the Layers Panel and then by clicking “Save as…”.

In terms of geometry, the result will be the same but, if you open the attribute table, you will see the result of the task:

attributes

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s