How to count points in polygons in QGIS using Python

In this task we will use python scripting to count the number of points which lie inside a polygon and store the result in a memory layer. The Processing algorithm Count points in polygon does the same thing, but we want to do the task without using it.

We will use two datasets available from Natural Earth:

  • the first one is a point shapefile compilation of world wide airports that is in the public domain (download it here);
  • the second one is a polygon shapefile containing the administrative areas of world wide countries (download it here).

If you download the datasets above, you will see this:

input

The first goal is to count the number of points in each feature of the polygon layer. Then, a new memory layer is generated, containing the same initial attributes plus one additional field that stores the result of the counting operation.

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
##Polygon_Vector=vector

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

layer1 = processing.getObject(Point_Vector)
crs = layer1.crs().toWkt()
layer2 = processing.getObject(Polygon_Vector)

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

points = [feat for feat in layer1.getFeatures()]
polygons = [feat for feat in layer2.getFeatures()]

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

for polygon in polygons:
    inAttr = polygon.attributes() # Input attributes
    poly_geom = polygon.geometry() #Input geometry
    engine = QgsGeometry.createGeometryEngine(poly_geom.geometry())
    engine.prepareGeometry()
    idsList = index.intersects(poly_geom.boundingBox())
    count = 0
    if len(idsList) > 0:
        req = QgsFeatureRequest().setFilterFids(idsList)
        for fit in layer1.getFeatures(req):
            geom = fit.geometry()
            if engine.contains(geom.geometry()):
                count += 1
    outGeom = QgsFeature()
    inAttr.append(count)
    outGeom.setAttributes(inAttr) # Output attributes
    outGeom.setGeometry(poly_geom) # Output geometry
    prov.addFeatures([outGeom]) # Output feature

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

In terms of geometry, the result will be the same but, if you select the Identify tool and click on any of the countries to examine the available attributes, you will see the number of points contained in it. For example, if we click on any point of Australia, we get a total of 17 airports:

identify

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

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