3

Like in a title. Let's say I have normal polygon that overlaps other. And in my case, when the other polygon has holes or is splitted (like on attached image), normal methods "intersection" or "intersects" don't work.

One of my try was something like this, but it failed.

"some_feature" is that normal polygon.

# if polygon has holes, create new one without them
if len(feature.geometry().asPolygon()) > 1:
    feat = QgsFeature()
    geom = feature.geometry().asPolygon()
    feat.setGeometry(QgsGeometry.fromPolygon([geom[0]]))
    print(some_feature.geometry().intersection(feat.geometry()).area())

Below are some images, that show cases I'm talking about.

On the left is case with holes and on the right is one object as two polygons.

Does anyone have idea, how to achieve such thing?

Joseph
  • 75,746
  • 7
  • 171
  • 282
DIV
  • 886
  • 1
  • 6
  • 18

1 Answers1

3

I don't know why normal methods "intersection" or "intersects" don't work in your script for polygons with holes or multi parts but, in my script, they did. I used polygon layers of next image (where attributes table belongs to multipart layer):

enter image description here

My complete script is:

registry = QgsMapLayerRegistry.instance()

polygon3_test_holes = registry.mapLayersByName('polygon3_test_holes')
polygon7_test = registry.mapLayersByName('polygon7_test')
polygon_multipart = registry.mapLayersByName('polygon_multipart')

featp3_hole = polygon3_test_holes[0].getFeatures().next()
featp7_test = polygon7_test[0].getFeatures().next()
feat_multip = polygon_multipart[0].getFeatures().next()

epsg = polygon3_test_holes[0].crs().postgisSrid()

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

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

prov = mem_layer.dataProvider()

if featp3_hole.geometry().intersects(featp7_test.geometry()):
    geom = featp3_hole.geometry().intersection(featp7_test.geometry())
    feat = QgsFeature()
    feat.setAttributes([1])
    feat.setGeometry(geom)
    prov.addFeatures([feat])

if featp3_hole.geometry().intersects(feat_multip.geometry()):

    geom = featp3_hole.geometry().intersection(feat_multip.geometry())
    feat = QgsFeature()
    feat.setAttributes([2])
    feat.setGeometry(geom)
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running the script at the Python Console of QGIS, I got memory layer (violet) of next image (where it was selected second feature):

enter image description here

Memory layer is a valid layer (it was corroborated directly at command line).

>>>layer = iface.activeLayer()  #memory layer
>>>new_feat = layer.getFeatures().next()
>>>new_feat.geometry().isGeosValid()
True

Editing Note:

To validate feature geometry:

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

for feat in feats:

    if feat.geometry().isGeosValid():
        pass
    else:
        print "geometry is not valid for ID: ", feat.attribute("ID")
xunilk
  • 29,891
  • 4
  • 41
  • 80
  • First of all, thanks for your time, xunilk. It helped me in some way, but not solved a problem. Anyway, now I know at least, that a polygons cause problems. I've used your script on new created polygons and result is like yours, but when I do the same with my "original" layers, I got nothing. So, now I have to figure out whats wrong with my polygons :/ – DIV Dec 20 '16 at 22:07
  • See my Editing Note. I put a script in it to validate geometry for each feature. Change feat.attribute("ID") for adapting to attribute name field in your case. – xunilk Dec 21 '16 at 11:47
  • It works. Some of my polygons are not valid. Do you know how I can fix this? Found something like v.clean and it's good for some of them, but maybe you know better way? – DIV Dec 21 '16 at 12:50
  • You can also use 'Topology Checker' plugin of QGIS. – xunilk Dec 21 '16 at 14:26