Merge pull request #3304 from nyalldawson/processing_clip

[FEATURE] Optimise processing clip algorithm
qgis · Aug 3, 2016 · c695571 · c695571
2 parents 420917b + d1aa03a
commit c695571
Show file tree

Hide file tree

Showing 4 changed files with 79 additions and 68 deletions.
diff --git a/python/plugins/processing/algs/qgis/Clip.py b/python/plugins/processing/algs/qgis/Clip.py
@@ -60,77 +60,88 @@ def defineCharacteristics(self):
         self.addOutput(OutputVector(Clip.OUTPUT, self.tr('Clipped')))
 
     def processAlgorithm(self, progress):
-        layerA = dataobjects.getObjectFromUri(
+        source_layer = dataobjects.getObjectFromUri(
             self.getParameterValue(Clip.INPUT))
-        layerB = dataobjects.getObjectFromUri(
+        mask_layer = dataobjects.getObjectFromUri(
             self.getParameterValue(Clip.OVERLAY))
 
         writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
-            layerA.pendingFields(),
-            layerA.dataProvider().geometryType(),
-            layerA.dataProvider().crs())
-
-        inFeatA = QgsFeature()
-        inFeatB = QgsFeature()
-        outFeat = QgsFeature()
-
-        index = vector.spatialindex(layerB)
-
-        selectionA = vector.features(layerA)
-
-        total = 100.0 / len(selectionA)
-
-        for current, inFeatA in enumerate(selectionA):
-            geom = inFeatA.geometry()
-            attrs = inFeatA.attributes()
-            intersects = index.intersects(geom.boundingBox())
-            first = True
-            found = False
-            if len(intersects) > 0:
-                for i in intersects:
-                    layerB.getFeatures(
-                        QgsFeatureRequest().setFilterFid(i)).nextFeature(
-                            inFeatB)
-                    tmpGeom = inFeatB.geometry()
-                    if tmpGeom.intersects(geom):
-                        found = True
-                        if first:
-                            outFeat.setGeometry(QgsGeometry(tmpGeom))
-                            first = False
-                        else:
-                            cur_geom = outFeat.geometry()
-                            new_geom = QgsGeometry(cur_geom.combine(tmpGeom))
-                            if new_geom.isEmpty() or 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.'))
-                                break
-
-                            outFeat.setGeometry(QgsGeometry(new_geom))
-                if found:
-                    cur_geom = outFeat.geometry()
-                    new_geom = QgsGeometry(geom.intersection(cur_geom))
+            source_layer.fields(),
+            source_layer.dataProvider().geometryType(),
+            source_layer.crs())
+
+        # first build up a list of clip geometries
+        clip_geoms = []
+        for maskFeat in vector.features(mask_layer, QgsFeatureRequest().setSubsetOfAttributes([])):
+            clip_geoms.append(maskFeat.geometry())
+
+        # are we clipping against a single feature? if so, we can show finer progress reports
+        if len(clip_geoms) > 1:
+            combined_clip_geom = QgsGeometry.unaryUnion(clip_geoms)
+            single_clip_feature = False
+        else:
+            combined_clip_geom = clip_geoms[0]
+            single_clip_feature = True
+
+        # use prepared geometries for faster intersection tests
+        engine = QgsGeometry.createGeometryEngine(combined_clip_geom.geometry())
+        engine.prepareGeometry()
+
+        tested_feature_ids = set()
+
+        for i, clip_geom in enumerate(clip_geoms):
+            input_features = [f for f in vector.features(source_layer, QgsFeatureRequest().setFilterRect(clip_geom.boundingBox()))]
+
+            if single_clip_feature:
+                total = 100.0 / len(input_features)
+            else:
+                total = 0
+
+            for current, in_feat in enumerate(input_features):
+                if not in_feat.geometry():
+                    continue
+
+                if in_feat.id() in tested_feature_ids:
+                    # don't retest a feature we have already checked
+                    continue
+
+                tested_feature_ids.add(in_feat.id())
+
+                if not engine.intersects(in_feat.geometry().geometry()):
+                    continue
+
+                if not engine.contains(in_feat.geometry().geometry()):
+                    cur_geom = in_feat.geometry()
+                    new_geom = combined_clip_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))
+                        int_com = in_feat.geometry().combine(new_geom)
+                        int_sym = in_feat.geometry().symDifference(new_geom)
+                        new_geom = 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 '
-                                                       'output features ignored due to '
-                                                       'invalid geometry.'))
-                        continue
-
-            progress.setPercentage(int(current * total))
+                else:
+                    # clip geometry totally contains feature geometry, so no need to perform intersection
+                    new_geom = in_feat.geometry()
+
+                try:
+                    out_feat = QgsFeature()
+                    out_feat.setGeometry(new_geom)
+                    out_feat.setAttributes(in_feat.attributes())
+                    writer.addFeature(out_feat)
+                except:
+                    ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
+                                           self.tr('Feature geometry error: One or more '
+                                                   'output features ignored due to '
+                                                   'invalid geometry.'))
+                    continue
+
+                if single_clip_feature:
+                    progress.setPercentage(int(current * total))
+
+            if not single_clip_feature:
+                # coarse progress report for multiple clip geometries
+                progress.setPercentage(100.0 * i / len(clip_geoms))
 
         del writer
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>