[FEATURE] Optimise processing clip algorithm

Before the algorithm was written to optimise clipping a few features against thousands of mask features. The revised algorithm is optimised for clipping thousands of input features against a few mask features. Given that this second operation is much more likely, it makes sense to optimise for this use case. I've also applied some other optimisations like taking advantage of spatial indexes on the providers, using prepared geometries and also only applying an intersection operation if the geometry isn't wholly contained by the mask geometry. Benchmarks: clipping roads layer with 1 million lines against 2 polygons before: 5 mins 30 seconds after: 10 seconds clipping address layer with 5 million points against 2 polygons before: 50 minutes after: 30 seconds
qgis · Aug 3, 2016 · 71ebdb8 · 71ebdb8
1 parent 8182f29
commit 71ebdb8
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 65 deletions.
diff --git a/python/plugins/processing/algs/qgis/Clip.py b/python/plugins/processing/algs/qgis/Clip.py
@@ -60,76 +60,69 @@ def defineCharacteristics(self):
         self.addOutput(OutputVector(Clip.OUTPUT, self.tr('Clipped')))
 
     def processAlgorithm(self, progress):
-        layerA = dataobjects.getObjectFromUri(
+        sourceLayer = dataobjects.getObjectFromUri(
             self.getParameterValue(Clip.INPUT))
-        layerB = dataobjects.getObjectFromUri(
+        maskLayer = dataobjects.getObjectFromUri(
             self.getParameterValue(Clip.OVERLAY))
 
         writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
-            layerA.pendingFields(),
+            sourceLayer.pendingFields(),
-            layerA.dataProvider().geometryType(),
+            sourceLayer.dataProvider().geometryType(),
-            layerA.dataProvider().crs())
+            sourceLayer.dataProvider().crs())
 
         inFeatA = QgsFeature()
         inFeatB = QgsFeature()
         outFeat = QgsFeature()
 
-        index = vector.spatialindex(layerB)
+        # first build up a list of clip geometries
-
+        clip_geoms = []
-        selectionA = vector.features(layerA)
+        for maskFeat in vector.features(maskLayer, QgsFeatureRequest().setSubsetOfAttributes([])):
-
+            clip_geoms.append(maskFeat.geometry())
-        total = 100.0 / len(selectionA)
+
-
+        if len(clip_geoms) > 1:
-        for current, inFeatA in enumerate(selectionA):
+            combined_clip_geom = QgsGeometry.unaryUnion(clip_geoms)
-            geom = inFeatA.geometry()
+        else:
-            attrs = inFeatA.attributes()
+            combined_clip_geom = clip_geoms[0]
-            intersects = index.intersects(geom.boundingBox())
+
-            first = True
+        # use prepared geometries for faster insection tests
-            found = False
+        engine = QgsGeometry.createGeometryEngine(combined_clip_geom.geometry())
-            if len(intersects) > 0:
+        engine.prepareGeometry()
-                for i in intersects:
+
-                    layerB.getFeatures(
+        input_features = [f for f in vector.features(sourceLayer, QgsFeatureRequest().setFilterRect(combined_clip_geom.boundingBox()))]
-                        QgsFeatureRequest().setFilterFid(i)).nextFeature(
+        total = 100.0 / len(input_features)
-                            inFeatB)
+        for current, in_feat in enumerate(input_features):
-                    tmpGeom = inFeatB.geometry()
+            if not in_feat.geometry():
-                    if tmpGeom.intersects(geom):
+                continue
-                        found = True
+
-                        if first:
+            if not engine.intersects(in_feat.geometry().geometry()):
-                            outFeat.setGeometry(QgsGeometry(tmpGeom))
+                continue
-                            first = False
+
-                        else:
+            if not engine.contains(in_feat.geometry().geometry()):
-                            cur_geom = outFeat.geometry()
+                cur_geom = in_feat.geometry()
-                            new_geom = QgsGeometry(cur_geom.combine(tmpGeom))
+                new_geom = combined_clip_geom.intersection(cur_geom)
-                            if new_geom.isEmpty() or new_geom.isGeosEmpty() or not new_geom.isGeosValid():
+                if new_geom.wkbType() == Qgis.WKBUnknown or QgsWKBTypes.flatType(new_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
-                                ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
+                    int_com = in_feat.geometry().combine(new_geom)
-                                                       self.tr('GEOS geoprocessing error: One or '
+                    int_sym = in_feat.geometry().symDifference(new_geom)
-                                                               'more input features have invalid '
+                    new_geom = int_com.difference(int_sym)
-                                                               'geometry.'))
+                    if new_geom.isGeosEmpty() or not new_geom.isGeosValid():
-                                break
-
-                            outFeat.setGeometry(QgsGeometry(new_geom))
-                if found:
-                    cur_geom = outFeat.geometry()
-                    new_geom = QgsGeometry(geom.intersection(cur_geom))
-                    if new_geom.wkbType() == Qgis.WKBUnknown or QgsWKBTypes.flatType(new_geom.geometry().wkbType()) == QgsWKBTypes.GeometryCollection:
-                        int_com = QgsGeometry(geom.combine(cur_geom))
-                        int_sym = QgsGeometry(geom.symDifference(cur_geom))
-                        new_geom = QgsGeometry(int_com.difference(int_sym))
-                        if new_geom.isGeosEmpty() or not new_geom.isGeosValid():
-                            ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
-                                                   self.tr('GEOS geoprocessing error: One or more '
-                                                           'input features have invalid geometry.'))
-                            continue
-                    try:
-                        outFeat.setGeometry(new_geom)
-                        outFeat.setAttributes(attrs)
-                        writer.addFeature(outFeat)
-                    except:
                         ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
-                                               self.tr('Feature geometry error: One or more '
+                                               self.tr('GEOS geoprocessing error: One or more '
-                                                       'output features ignored due to '
+                                                       'input features have invalid geometry.'))
-                                                       'invalid geometry.'))
+            else:
-                        continue
+                # clip geometry totally contains feature geometry, so no need to perform intersection
+                new_geom = in_feat.geometry()
+
+            try:
+                outFeat = QgsFeature()
+                outFeat.setGeometry(new_geom)
+                outFeat.setAttributes(in_feat.attributes())
+                writer.addFeature(outFeat)
+            except:
+                ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
+                                       self.tr('Feature geometry error: One or more '
+                                               'output features ignored due to '
+                                               'invalid geometry.'))
+                continue
 
             progress.setPercentage(int(current * total))
 

diff --git a/python/plugins/processing/tests/testdata/expected/clip_lines_by_multipolygon.gml b/python/plugins/processing/tests/testdata/expected/clip_lines_by_multipolygon.gml
@@ -23,7 +23,7 @@
   </gml:featureMember>
   <gml:featureMember>
     <ogr:clip_lines_by_multipolygon fid="lines.3">
-      <ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>3,1 4,1</gml:coordinates></gml:LineString></ogr:geometryProperty>
+      <ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>4,1 3,1</gml:coordinates></gml:LineString></ogr:geometryProperty>
     </ogr:clip_lines_by_multipolygon>
   </gml:featureMember>
   <gml:featureMember>

diff --git a/python/plugins/processing/tests/testdata/expected/clip_multipolygons_by_polygons.gml b/python/plugins/processing/tests/testdata/expected/clip_multipolygons_by_polygons.gml
@@ -13,15 +13,15 @@
 
   <gml:featureMember>
     <ogr:clip_multipolygons_by_polygons fid="multipolys.0">
-      <ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:4326"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>4.0,1.666666666666667 4,1 2,1 2,2 3,2 4.0,1.666666666666667</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
+      <ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:4326"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>2 1, 2 2, 3 2, 4 1.66666666666667, 4 1, 2 1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
       <ogr:Bname>Test</ogr:Bname>
       <ogr:Bintval>1</ogr:Bintval>
       <ogr:Bfloatval>0.123</ogr:Bfloatval>
     </ogr:clip_multipolygons_by_polygons>
   </gml:featureMember>
   <gml:featureMember>
     <ogr:clip_multipolygons_by_polygons fid="multipolys.1">
-      <ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:4326"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>8,1 8,0 7,0 7,1 8,1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
+      <ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:4326"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>7 0, 7 1, 8 1, 8 0, 7 0</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
     </ogr:clip_multipolygons_by_polygons>
   </gml:featureMember>
   <gml:featureMember>

diff --git a/python/plugins/processing/tests/testdata/expected/clip_polys_by_multipolygon.gml b/python/plugins/processing/tests/testdata/expected/clip_polys_by_multipolygon.gml
@@ -21,14 +21,14 @@
   </gml:featureMember>
   <gml:featureMember>
     <ogr:clip_polys_by_multipolygon fid="polys.3">
-      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>8,0 7,0 7,1 8,1 8,0</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
+      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>7 1, 8 1, 8 0, 7 0, 7 1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
       <ogr:name>ASDF</ogr:name>
       <ogr:intval>0</ogr:intval>
     </ogr:clip_polys_by_multipolygon>
   </gml:featureMember>
   <gml:featureMember>
     <ogr:clip_polys_by_multipolygon fid="polys.5">
-      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>3,2 4.0,1.666666666666667 4,1 2,1 2,2 3,2</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
+      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>2 1, 2 2, 3 2, 4 1.66666666666667, 4 1, 2 1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
       <ogr:name>elim</ogr:name>
       <ogr:intval>2</ogr:intval>
       <ogr:floatval>3.33</ogr:floatval>