Skip to content

Commit

Permalink
Better clip for very separate features
Browse files Browse the repository at this point in the history
  • Loading branch information
nyalldawson committed Aug 3, 2016
1 parent 71ebdb8 commit d1aa03a
Showing 1 changed file with 65 additions and 47 deletions.
112 changes: 65 additions & 47 deletions python/plugins/processing/algs/qgis/Clip.py
Expand Up @@ -60,70 +60,88 @@ def defineCharacteristics(self):
self.addOutput(OutputVector(Clip.OUTPUT, self.tr('Clipped'))) self.addOutput(OutputVector(Clip.OUTPUT, self.tr('Clipped')))


def processAlgorithm(self, progress): def processAlgorithm(self, progress):
sourceLayer = dataobjects.getObjectFromUri( source_layer = dataobjects.getObjectFromUri(
self.getParameterValue(Clip.INPUT)) self.getParameterValue(Clip.INPUT))
maskLayer = dataobjects.getObjectFromUri( mask_layer = dataobjects.getObjectFromUri(
self.getParameterValue(Clip.OVERLAY)) self.getParameterValue(Clip.OVERLAY))


writer = self.getOutputFromName(self.OUTPUT).getVectorWriter( writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(
sourceLayer.pendingFields(), source_layer.fields(),
sourceLayer.dataProvider().geometryType(), source_layer.dataProvider().geometryType(),
sourceLayer.dataProvider().crs()) source_layer.crs())

inFeatA = QgsFeature()
inFeatB = QgsFeature()
outFeat = QgsFeature()


# first build up a list of clip geometries # first build up a list of clip geometries
clip_geoms = [] clip_geoms = []
for maskFeat in vector.features(maskLayer, QgsFeatureRequest().setSubsetOfAttributes([])): for maskFeat in vector.features(mask_layer, QgsFeatureRequest().setSubsetOfAttributes([])):
clip_geoms.append(maskFeat.geometry()) 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: if len(clip_geoms) > 1:
combined_clip_geom = QgsGeometry.unaryUnion(clip_geoms) combined_clip_geom = QgsGeometry.unaryUnion(clip_geoms)
single_clip_feature = False
else: else:
combined_clip_geom = clip_geoms[0] combined_clip_geom = clip_geoms[0]
single_clip_feature = True


# use prepared geometries for faster insection tests # use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(combined_clip_geom.geometry()) engine = QgsGeometry.createGeometryEngine(combined_clip_geom.geometry())
engine.prepareGeometry() engine.prepareGeometry()


input_features = [f for f in vector.features(sourceLayer, QgsFeatureRequest().setFilterRect(combined_clip_geom.boundingBox()))] tested_feature_ids = set()
total = 100.0 / len(input_features)
for current, in_feat in enumerate(input_features): for i, clip_geom in enumerate(clip_geoms):
if not in_feat.geometry(): input_features = [f for f in vector.features(source_layer, QgsFeatureRequest().setFilterRect(clip_geom.boundingBox()))]
continue

if single_clip_feature:
if not engine.intersects(in_feat.geometry().geometry()): total = 100.0 / len(input_features)
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 = 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.'))
else: else:
# clip geometry totally contains feature geometry, so no need to perform intersection total = 0
new_geom = in_feat.geometry()

for current, in_feat in enumerate(input_features):
try: if not in_feat.geometry():
outFeat = QgsFeature() continue
outFeat.setGeometry(new_geom)
outFeat.setAttributes(in_feat.attributes()) if in_feat.id() in tested_feature_ids:
writer.addFeature(outFeat) # don't retest a feature we have already checked
except: continue
ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
self.tr('Feature geometry error: One or more ' tested_feature_ids.add(in_feat.id())
'output features ignored due to '
'invalid geometry.')) if not engine.intersects(in_feat.geometry().geometry()):
continue continue


progress.setPercentage(int(current * total)) 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 = 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.'))
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 del writer

0 comments on commit d1aa03a

Please sign in to comment.