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

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

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