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:


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):


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

layer1 = processing.getObject(Point_Vector)
crs =
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

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:

for polygon in polygons:
    inAttr = polygon.attributes() # Input attributes
    poly_geom = polygon.geometry() #Input geometry
    engine = QgsGeometry.createGeometryEngine(poly_geom.geometry())
    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()
    outGeom.setAttributes(inAttr) # Output attributes
    outGeom.setGeometry(poly_geom) # Output geometry
    prov.addFeatures([outGeom]) # Output feature

# Add the layer to the Layers panel

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:


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


Leave a Reply

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

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

Google+ photo

You are commenting using your Google+ 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 )


Connecting to %s