Counting points in polygon by time in QGIS

Question

I want to count the points in the polygon by time category.

However, the point have time category.

When I want to use Vector > Analysis Tools > Count points in polygon option, it does not consider the time category. It counts all the points and does not give separate by time.

Shortly, how could I count the points in a polygon by date/time?

Answer 1

Because I think this is a very interesting question and a good use case, I wrote this little processing tool. Copy paste the code to a .py file and place it in C:\Users\Username\AppData\Roaming\QGIS\QGIS3\profiles\default\processing\scripts (respectively Linux/Mac paths). You can then find in your processing toolbox within scripts -> from gisse -> Count Points in Polygon by Datetime.

A modified version of this script is now also part of the ProcessX-PlugIn you can get in the official QGIS Plug-In repository. You can find it in your processing toolbox in ProcessX -> Vector - Creation -> Create Timepolygon with Pointcount. Alternatively you may also want to try Count Features in Features with Condition or Count Points in Polygons With Condition algorithms.

It requires:

• A field of type DateTime in your Pointslayer
• A given start- and end datetime you can enter as string
• A given interval in seconds

It does:

• Count the Points in Polygons by the given datetime of the points
• Creates a new polygon-copy for each datetime-range and adds the pointcount, fromdatetime and todateime as attributes
• If Point- and Polygonlayer are not in the same CRS, the Pointlayer is temporarily reprojected to the polygonlayers crs (no output is generated for that one)

The script:

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsField, QgsFeature, QgsProcessing, QgsExpression, QgsGeometry, QgsPoint, QgsFields,
                       QgsWkbTypes,
                       QgsFeatureSink, QgsFeatureRequest, QgsProcessingAlgorithm, QgsSpatialIndex,
                       QgsProcessingParameterFeatureSink, QgsProcessingParameterDateTime, QgsProcessingParameterField,
                       QgsProcessingParameterFeatureSource, QgsProcessingParameterEnum, QgsProcessingParameterString,
                       QgsProcessingParameterNumber)
import processing
from datetime import *
from math import ceil


class CountPointsInPolygonByTime(QgsProcessingAlgorithm):
    POLYGON_LYR = 'POLYGON_LYR'
    POINT_LYR = 'POINT_LYR'
    DATETIME_FIELD = 'DATETIME_FIELD'
    START_DATETIME = 'START_DATETIME'
    END_DATETIME = 'END_DATETIME'
    INTERVALSEC = 'INTERVALSEC'
    OUTPUT = 'OUTPUT'

    def initAlgorithm(self, config=None):

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.POLYGON_LYR, self.tr('Polygon'), [QgsProcessing.TypeVectorPolygon], 'polygons'))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.POINT_LYR, self.tr('Point'), [QgsProcessing.TypeVectorPoint], 'points'))
        self.addParameter(
            QgsProcessingParameterField(
                self.DATETIME_FIELD, self.tr('Datetime Field'), 'datetime', 'POINT_LYR'))
        self.addParameter(
            QgsProcessingParameterString(
                self.START_DATETIME, self.tr('Start Datetime in YYYY-MM-DD HH:MM:SS format'), '2020-01-01 00:00:00'))
        self.addParameter(
            QgsProcessingParameterString(
                self.END_DATETIME, self.tr('End Datetime in YYYY-MM-DD HH:MM:SS format'), '2020-01-10 23:59:59'))
        self.addParameter(
            QgsProcessingParameterNumber(
                self.INTERVALSEC, self.tr('Interval in Seconds'), 0,
                86400))  # Indicator as number. 0=Int, 1 would be double; 1=default number
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT, self.tr('TimePolygons with Pointcount')))  # Output

    def processAlgorithm(self, parameters, context, feedback):
        lyr_polygons = self.parameterAsLayer(parameters, self.POLYGON_LYR, context)
        lyr_points = self.parameterAsLayer(parameters, self.POINT_LYR, context)
        fld_time = self.parameterAsString(parameters, self.DATETIME_FIELD, context)
        start_date_string = self.parameterAsString(parameters, self.START_DATETIME, context)
        end_date_string = self.parameterAsString(parameters, self.END_DATETIME, context)
        intervalsec = self.parameterAsInt(parameters, self.INTERVALSEC, context)

        if lyr_polygons.sourceCrs() != lyr_points.sourceCrs():
            reproj = processing.run("native:reprojectlayer",
                                        {'INPUT': lyr_points,
                                         'TARGET_CRS': lyr_polygons.sourceCrs(),
                                         'OUTPUT': 'memory:Reprojected'})
            lyr_points = reproj['OUTPUT']

        fields = lyr_polygons.fields()
        fields.append(QgsField('from_datetime', QVariant.DateTime))
        fields.append(QgsField('to_datetime', QVariant.DateTime))
        fields.append(QgsField('pointcount', QVariant.Int, len=0))

        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                               fields, lyr_polygons.wkbType(),
                                               lyr_polygons.sourceCrs())

        start_date = datetime.strptime(start_date_string, '%Y-%m-%d %H:%M:%S')
        end_date = datetime.strptime(end_date_string, '%Y-%m-%d %H:%M:%S')
        total_seconds = int((end_date - start_date).total_seconds())
        idx_points = QgsSpatialIndex(lyr_points.getFeatures())

        required_iterations = ceil(total_seconds / intervalsec)
        total = 100.0 / (lyr_polygons.featureCount() * required_iterations) if lyr_polygons.featureCount() else 0
        current = 0

        for current_interval in range(0, total_seconds, intervalsec):
            current_start_datetime = start_date + timedelta(seconds=current_interval)
            current_end_datetime = (
                        start_date + timedelta(seconds=current_interval + intervalsec) - timedelta(seconds=1))
            for polygon in lyr_polygons.getFeatures():
                current += 1
                new_feat = QgsFeature(fields)
                new_feat.setGeometry(polygon.geometry())
                idx = 0
                for attr in polygon.attributes():
                    new_feat[idx] = attr
                    idx += 1
                new_feat['from_datetime'] = current_start_datetime.strftime('%Y-%m-%d %H:%M:%S')
                new_feat['to_datetime'] = current_end_datetime.strftime('%Y-%m-%d %H:%M:%S')
                new_feat['pointcount'] = 0
                for pointid in idx_points.intersects(polygon.geometry().boundingBox()):
                    point = lyr_points.getFeature(pointid)
                    if feedback.isCanceled():
                        break
                    if point[fld_time] >= current_start_datetime and point[fld_time] <= current_end_datetime:
                        if point.geometry().intersects(polygon.geometry()):
                            new_feat['pointcount'] += 1
                            idx_points.deleteFeature(
                                point)  # dont count a point twice, removing it from the index speeds up the code around 25%

                if feedback.isCanceled():
                    break

                sink.addFeature(new_feat, QgsFeatureSink.FastInsert)
                feedback.setProgress(int(current * total))

        return {self.OUTPUT: dest_id}  # Return result of algorithm

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return CountPointsInPolygonByTime()

    def name(self):
        return 'CountPointsInPolygonByTime'

    def displayName(self):
        return self.tr('Count Points in Polygon by Datetime')

    def group(self):
        return self.tr('FROM GISSE')

    def groupId(self):
        return 'from_gisse'

    def shortHelpString(self):
        return self.tr('This Algorithm counts points in polygons by a given datetime condition')

Answer 2

Use Join attributes by location with your Points as Base layer and Polygons as Join Layer

Then use Statistics By Category with Joined Layer as the Input Vector Layer

id as Field to Calculate Statistics On

For the Fields with categories parameter use id_2 (which is the id field from the polygons) and date

The output will be a table that has a count for each polyon id and each date (plus several other statistics fields which you can ignore or delete). This table can be joined back to your polygon layer using id_2 and a selected date. Otherwise, you could pivot this table and join the whole pivot table back.

Answer 3

Using QGIS expressions

create a new field in the polygon layer, and populate it with this expression:

array_to_string(
    with_variable('arraydate',
        aggregate(
        layer:='Point', 
        aggregate:='array_agg', 
        expression:="date",
        order_by:="date",
        filter:= intersects( $geometry, geometry(@parent))),
array_foreach(
      array_distinct(@arraydate),
             @element||': '||array_count(@arraydate,@element))
))

Answer 4

Although it creates new data files, which may not be desirable, I would:

1. Use the Split Vector Layer tool to split the point file into multiple layers each with points with a unique date
2. Run Count Points in Polygon for each unique date layer

You'll end up with a polygon file with multiple fields, each representing the point count for that specific date.