How to count points within the current Print composer atlas feature in QGIS 2.8?

Question

In the print composer when generating an atlas, does anyone know how to show he feature count per each atlas coverage feature?

For example, I have generated an atlas where each map is a different province (in Cambodia) containing a village point layer. When I click the sigma sign at the bottom of the legend for villages, the feature count will calculate the entire country's number of villages, not for the individual province's number of villages.

Does anyone know how to get a feature count in for the province showing in the atlas?

Answer 1

This works for me: use the aggregate() expression in a text label, with within($geometry,@atlas_geometry) as a filter for points within your atlas feature.

The full syntax in your text label would be

[% aggregate('LAYER_NAME','count',"UNIQUE_NON_NULL_ATTRIBUTE_FOR_COUNTING",within($geometry,@atlas_geometry)) %]

That will count all the points in your specified layer - by an unique non-null attribute, so each point is distinct - but only if the point geometry is within the atlas geometry.

The [% %] is required to make the expression work. And note the quotation marks - you must use double quotes for the attribute name.

As a bonus you can add further filters, so to address one of the comments above, you can modify the filter to say within($geometry,@atlas_geometry) AND "village" = 'A', or whatever combination.

I haven't extensively tested this out, especially on large/complex datasets, but it seems to do the job.

Answer 2

You need to create a script counting point within the current atlas geometry.

In the expression function editor add this code:

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args="auto", group='Custom')
def countPointsInCurrentAltlasFeature(pointLayerName, geomAtlas, feature, parent):
    # If point geom is empty, return 0
    if (geomAtlas is None):
        return 0

    # Get point layer reference from layername
    pointLayer = QgsMapLayerRegistry.instance().mapLayersByName(pointLayerName)[0]

    # Raise if layer not found
    if pointLayer is None:
        raise Exception("Layer not found: " + pointLayerName)

    # Count point within current Atlas feature
    countPoint = 0
    for pointFeature in pointLayer.getFeatures():
        pointGeom = pointFeature.geometry()
        if (pointGeom is None):
            continue
        if pointGeom.within(geomAtlas):
            countPoint += 1

    return countPoint

Use the new function as an expression for a Print Composer label like:

Villages in province: [% countPointsInCurrentAltlasFeature( 'Villages', $atlasgeometry)%]