2

I am trying to script the clipping of overlapping areas from a set of polygons. After having fetched features (outerFeature) a loop (innerFeature) is to check if other features intersect with area > 0. I then try to clip the intersecting area. Apparently difference does not give the right result.

#Intersecting features
for outerFeature in pr.getFeatures():
    #Looping through all features and check if intersect
    for innerFeature in pr.getFeatures():
        geom = None
        if outerFeature.id() != innerFeature.id():
            geom = outerFeature.geometry().intersection( innerFeature.geometry() )
            #If intersection area > 0, then construct difference and update geometry
            if geom.area() > 0:
                print("intersecting and clipping")
                outerFeature.setGeometry( geom.difference( outerFeature.geometry() ) )
                vl.updateExtents()

Any advise is highly welcome.

Best regards

Sebastian

  • Is there a reason you're trying to do this programmatically? There are QGIS tools to address this issue: gis.stackexchange.com/questions/201668/… – Alex Aug 23 at 16:10
  • The QGIS tool in mention does not allow for automated resolution to my understanding. I first want to clean programmatically and then proceed with topological editing. – Sebastian Aug 27 at 8:37
  • The geometry checker tool should have automated resolution. See the bottom of the help file here: docs.qgis.org/3.4/en/docs/user_manual/plugins/core_plugins/… – Alex Aug 27 at 13:39
  • I saw that and tried it, but it did not resolve the conflicts. Suggested is to remove "Overlap with sewershed_raw at feature 31" -> Error is obsolete. This occurs for most of the several errors. When checking again they reoccur also in the same spots. – Sebastian Aug 27 at 18:41
0

It is not clear for me what do you want to do with above code but, when you have overlapped features in the same layer you can use itertools python module for avoiding this kind of lines:

if outerFeature.id() != innerFeature.id()

for repeated features. So, your code can be written as follow for producing a memory layer of all intersection of primary features.

import itertools

layer = iface.activeLayer()

features = [ feature for feature in layer.getFeatures() ]

n = layer.featureCount()

list = range(n)

intersections = [features[i].geometry().intersection(features[j].geometry()).asWkt()
                 for i,j in itertools.combinations(list, 2)
                 if features[i].geometry().intersects(features[j].geometry()) ]

epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(intersections)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(intersections[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

By using shapefile with overlapped features of following image:

enter image description here

after running above code, I got following result:

enter image description here

However, these interceptions can also intercepts each other or with another primary features. In this case, for a complete clipping, you can use 'saga:polygonselfintersection' method from Processing in a very concise script as follow:

import processing

polys = QgsVectorLayer("/home/zeito/pyqgis_data/overlapped_features.shp","self_intersections","ogr")
output = "/home/zeito/pyqgis_data/intersection.shp"
processing.runAndLoadResults("saga:polygonselfintersection", 
                {'POLYGONS':polys, 
                 'INTERSECT':output})

After running above code, result obtained was:

enter image description here

In attributes table was created a field (named ID) where it was placed all possible intersections for primary features.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.