Bug report #13608

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

Added by Adrian Klink about 5 years ago. Updated almost 5 years ago.

Status:Closed
Priority:Normal
Assignee:-
Category:-
Affected QGIS version:2.10.1 Regression?:No
Operating System: Easy fix?:No
Pull Request or Patch supplied:No Resolution:fixed/implemented
Crashes QGIS or corrupts data:No Copied to github as #:21647

Description

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))

Associated revisions

Revision 34dc3143
Added by Nyall Dawson about 5 years ago

Fix exporting geometry collections to WKT

Child types were incorrectly being dropped when the collection
consisted of mixed geometry types (eg line & polygon) (refs #13608)

History

#1 Updated by Jukka Rahkonen about 5 years ago

Paste also the result geometry from ArcGIS.

#2 Updated by Nyall Dawson about 5 years ago

  • Status changed from Open to Feedback

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.

#3 Updated by Adrian Klink about 5 years ago

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))

#4 Updated by Adrian Klink about 5 years ago

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)

#5 Updated by Giovanni Manghi almost 5 years ago

  • Resolution set to fixed/implemented
  • Status changed from Feedback to Closed

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.

Also available in: Atom PDF