Possible Bug in GEOS when calculating intersection of two polygons with inner rings?

[qgib]

Author Name: Adrian Klink (@aklink)
Original Redmine Issue: 13608
Affected QGIS version: 2.10.1

---

As far as I understand, GEOS is responsible for calculationg intersections between polygons in QGIS. When I try to intersect two polygons with inner rings (holes), intersection fails. Is this a possible bug in GEOS module in QGIS? Using ArcGIS creates a valid geometry. I have tried using QGIS 2.10 Pisa 64bit Windows 7 and QGIS 2.10 Pisa 32bit Windows 8.1, using Vector -> Intersection via menu and/or Python scripting. Either I get no geometry or invalid GeometryCollection.

Can someone double check, please?

Here are the polygons (as WKT String, can be imported into QGIS via QuickWKT Plugin):

Polygon ((32363196.675240271 5467660.3907237295, 32363272.961243533 5467641.608177687, 32363283.202260282 5467556.0808803849, 32363090.88836604 5467558.0382346781, 32363089.014254488 5467625.9334499743, 32363196.675240271 5467660.3907237295),(32363160.68400000035762787 5467605.09290000051259995, 32363150.1810000017285347 5467599.87690000049769878, 32363135.58999999985098839 5467593.77390000037848949, 32363124.54699999839067459 5467587.28889999911189079, 32363160.21000000089406967 5467594.69590000063180923, 32363157.51599999889731407 5467596.38389999978244305, 32363156.73499999940395355 5467598.88590000011026859, 32363158.03199999779462814 5467602.30690000019967556, 32363160.68400000035762787 5467605.09290000051259995))

intersection with:

Polygon ((32363196.675240271 5467660.3907237295, 32363272.961243533 5467641.608177687, 32363283.202260282 5467556.0808803849, 32363090.88836604 5467558.0382346781, 32363089.014254488 5467625.9334499743, 32363196.675240271 5467660.3907237295),(32363262.36500000208616257 5467580.39589999988675117,32363173.82899999991059303 5467597.48589999973773956, 32363160.21000000089406967 5467594.69590000063180923, 32363124.54699999839067459 5467587.28889999911189079, 32363108.6939999982714653 5467581.39589999988675117, 32363262.36500000208616257 5467580.39589999988675117))

[qgib]

Author Name: Jukka Rahkonen (Jukka Rahkonen)

---

Paste also the result geometry from ArcGIS.

[qgib]

Author Name: Nyall Dawson (@nyalldawson)

---

Thanks - that's helped me track down a possibly related bug in geometry collections. Looks like that intersection operation results in a collection of a string and polygon.

Can you please confirm which tool/processing algorithm/menu item you are using to perform the intersection? That will probably need to be updated to handle this case.

---

• status_id was changed from Open to Feedback

[qgib]

Author Name: Adrian Klink (@aklink)

---

Nyall Dawson wrote:

Thanks - that's helped me track down a possibly related bug in geometry collections. Looks like that intersection operation results in a collection of a string and polygon.

Can you please confirm which tool/processing algorithm/menu item you are using to perform the intersection? That will probably need to be updated to handle this case.

Originally this WKT was taken from a Shapefile. This problem (empty geometry) occured, when using:

QGIS Menu -> Vector -> Geoprocessing Tools -> Intersect

The second approach (resulting in invalid GeometryCollection) was using Python Scripting in QGIS:

    layer1 = iface.addVectorLayer(first_file, "first_layer", "ogr")
    layer2 = iface.addVectorLayer(second_file, "second_layer", "ogr")
    iter1 = layer1.getFeatures()
        for feature1 in iter1:
            geom1 = feature1.geometry()
            for feature2 in iter2:
                geom2 = feature2.geometry()
                if geom1.intersects(geom2):
                    geomintersect = geom1.intersection(geom2)
                    print "%s" % (geomintersect.exportToWkt())

Resulting Geometry in ArcGIS is (Polygon only, since input have been 2 polygons and intersecting line of inner rings has no area, thus is no polygon):

Polygon ((32363272.96124353259801865 5467641.60817768704146147, 32363283.20226028189063072 5467556.08088038489222527, 32363090.88836603984236717 5467558.03823467809706926, 32363089.01425448805093765 5467625.93344997335225344, 32363196.67524027079343796 5467660.39072373043745756, 32363272.96124353259801865 5467641.60817768704146147),(32363108.6939999982714653 5467581.39590000081807375, 32363262.36500000208616257 5467580.39589999988675117, 32363173.82899999991059303 5467597.48589999973773956, 32363160.21000000089406967 5467594.69590000063180923, 32363157.51599999889731407 5467596.38389999978244305, 32363156.73499999940395355 5467598.88590000104159117, 32363158.03199999779462814 5467602.30690000019967556, 32363160.68400000035762787 5467605.09289999958127737, 32363150.1810000017285347 5467599.87690000049769878, 32363135.58999999985098839 5467593.77389999944716692, 32363124.54699999839067459 5467587.28889999911189079, 32363108.6939999982714653 5467581.39590000081807375))

[qgib]

Author Name: Adrian Klink (@aklink)

---

Adrian Klink wrote:

Polygon only, since input have been 2 polygons and intersecting line of inner rings has no area, thus is no polygon

Since I mentioned that ArcGIS and QGIS behave different on intersections of polygons (result is polygon only in ArcGIS) there is also a second case (I currently have no example, but I may provide one next week):

Two polygons touching each other at a line string behave different in QGIS and ArcGIS:

• ArcGIS: no intersection (touching line string has no area - thus no polygon)
• QGIS: intersection resulting in a line string (same behaviour as with above case which was resulting in a GeometryCollection of Polygon and LineString, but this time the outer ring and not the inner rings are affected)

[qgib]

Author Name: Giovanni Manghi (@gioman)

---

this was caused by a change in the qgis code that left several tools in the vector menu "broken" when the result contained a collection, this was fixed recently in qgis master.

---

• resolution was changed from to fixed/implemented
• status_id was changed from Feedback to Closed