Faster PointsInPolygon(+Unique/Weighted) algorithms

nyalldawson · nyalldawson · commit 596b56c44130 · 2015-11-13T13:24:08.000+11:00
- Avoid use of seperate feature requests for every point
- Use GEOS prepared geometries when testing for point containment

Quick stats:

~1500 point layer

BEFORE: 17 seconds
AFTER: 3 seconds

~900k point layer

BEFORE: 30 mins = canceled at 20%
AFTER: 2.5 mins = 100% complete
diff --git a/python/plugins/processing/algs/qgis/PointsInPolygon.py b/python/plugins/processing/algs/qgis/PointsInPolygon.py
@@ -76,26 +76,25 @@ def processAlgorithm(self, progress):
         geom = QgsGeometry()
 
         current = 0
-        hasIntersections = False
 
         features = vector.features(polyLayer)
         total = 100.0 / float(len(features))
         for ftPoly in features:
             geom = ftPoly.geometry()
+            engine = QgsGeometry.createGeometryEngine(geom.geometry())
+            engine.prepareGeometry()
+
             attrs = ftPoly.attributes()
 
             count = 0
-            hasIntersections = False
             points = spatialIndex.intersects(geom.boundingBox())
             if len(points) > 0:
-                hasIntersections = True
-
-            if hasIntersections:
-                for i in points:
-                    request = QgsFeatureRequest().setFilterFid(i)
-                    ftPoint = pointLayer.getFeatures(request).next()
-                    tmpGeom = QgsGeometry(ftPoint.geometry())
-                    if geom.contains(tmpGeom):
+                request = QgsFeatureRequest().setFilterFids(points)
+                fit = pointLayer.getFeatures(request)
+                ftPoint = QgsFeature()
+                while fit.nextFeature(ftPoint):
+                    tmpGeom = ftPoint.geometry()
+                    if engine.contains(tmpGeom.geometry()):
                         count += 1
 
             outFeat.setGeometry(geom)
diff --git a/python/plugins/processing/algs/qgis/PointsInPolygonUnique.py b/python/plugins/processing/algs/qgis/PointsInPolygonUnique.py
@@ -80,26 +80,25 @@ def processAlgorithm(self, progress):
         geom = QgsGeometry()
 
         current = 0
-        hasIntersections = False
 
         features = vector.features(polyLayer)
         total = 100.0 / float(len(features))
         for ftPoly in features:
             geom = ftPoly.geometry()
+            engine = QgsGeometry.createGeometryEngine(geom.geometry())
+            engine.prepareGeometry()
+
             attrs = ftPoly.attributes()
 
             classes = []
-            hasIntersections = False
             points = spatialIndex.intersects(geom.boundingBox())
             if len(points) > 0:
-                hasIntersections = True
-
-            if hasIntersections:
-                for i in points:
-                    request = QgsFeatureRequest().setFilterFid(i)
-                    ftPoint = pointLayer.getFeatures(request).next()
+                request = QgsFeatureRequest().setFilterFids(points)
+                fit = pointLayer.getFeatures(request)
+                ftPoint = QgsFeature()
+                while fit.nextFeature(ftPoint):
                     tmpGeom = QgsGeometry(ftPoint.geometry())
-                    if geom.contains(tmpGeom):
+                    if engine.contains(tmpGeom.geometry()):
                         clazz = ftPoint.attributes()[classFieldIndex]
                         if clazz not in classes:
                             classes.append(clazz)
diff --git a/python/plugins/processing/algs/qgis/PointsInPolygonWeighted.py b/python/plugins/processing/algs/qgis/PointsInPolygonWeighted.py
@@ -86,27 +86,25 @@ def processAlgorithm(self, progress):
         geom = QgsGeometry()
 
         current = 0
-        hasIntersections = False
-
         features = vector.features(polyLayer)
         total = 100.0 / float(len(features))
         for ftPoly in features:
             geom = ftPoly.geometry()
+            engine = QgsGeometry.createGeometryEngine(geom.geometry())
+            engine.prepareGeometry()
+
             attrs = ftPoly.attributes()
 
             count = 0
-            hasIntersections = False
             points = spatialIndex.intersects(geom.boundingBox())
             if len(points) > 0:
-                hasIntersections = True
-
-            if hasIntersections:
                 progress.setText(unicode(len(points)))
-                for i in points:
-                    request = QgsFeatureRequest().setFilterFid(i)
-                    ftPoint = pointLayer.getFeatures(request).next()
+                request = QgsFeatureRequest().setFilterFids(points)
+                fit = pointLayer.getFeatures(request)
+                ftPoint = QgsFeature()
+                while fit.nextFeature(ftPoint):
                     tmpGeom = QgsGeometry(ftPoint.geometry())
-                    if geom.contains(tmpGeom):
+                    if engine.contains(tmpGeom.geometry()):
                         weight = unicode(ftPoint.attributes()[fieldIdx])
                         try:
                             count += float(weight)
