Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #3304 from nyalldawson/processing_clip 
[FEATURE] Optimise processing clip algorithm
  • Loading branch information
@nyalldawson
nyalldawson committed Aug 3, 2016
2 parents 420917b + d1aa03a commit c695571
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 68 deletions.
137 changes: 74 additions & 63 deletions python/plugins/processing/algs/qgis/Clip.py
View file
Expand Up @@ -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
2 changes: 1 addition & 1 deletion python/plugins/processing/tests/testdata/expected/clip_lines_by_multipolygon.gml
View file
Expand Up @@ -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>
Expand Down
4 changes: 2 additions & 2 deletions python/plugins/processing/tests/testdata/expected/clip_multipolygons_by_polygons.gml
View file
Expand Up @@ -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>
Expand Down
4 changes: 2 additions & 2 deletions python/plugins/processing/tests/testdata/expected/clip_polys_by_multipolygon.gml
View file
Expand Up @@ -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>
Expand Down

0 comments on commit c695571

Please sign in to comment.