1

I have three simple geometries (download):

1) Province boundaries of a province in Egypt (Al Wadi Al Jadid).
2) A 10 x 10 degree polygon (index = 417).
3) another 10 x 10 degree polygon (index = 416).

enter image description here

In ArcMap, I get valid results for the intersection of the province boundary and gridcell 416. I also get valid results for the intersection of the province boundary and gridcell 417.

In QGIS (And Shapely), The intersection of the province boundary and the gridcell 416 produces a valid result.

In QGIS (And Shapely), The intersection of the province boundary and the gridcell 417 produces unexpected results:

  • Empty result in QGIS
  • GeometryCollection with LineString and Polygon in Shapely.

    GEOMETRYCOLLECTION (LINESTRING (30 22, 30 21.99920082), POLYGON ((30 22, 30 27.62037118927962, 30.53003883 27.6014328, 30.86449623 26.79118919, 31.15589523 26.82732773, 31.29156494 26.83351135, 31.44809151 26.63268852, 31.76193428 26.29763031, 31.95477867 26.13333893, 32.03946686 26.06095505, 32.21749496 25.97387695, 32.24517059 25.98464203, 32.32149506 25.98672104, 32.44264603 26.04069328, 32.67482376 26.10331726, 32.73316956 26.07575607, 32.73933411 26.05672836, 32.72978973 25.98652077, 32.72766876 25.98363495, 32.72722244 25.97396469, 32.68310928 25.91794777, 32.66844177 25.84004784, 32.64354324 25.79011726, 32.53351593 25.69179344, 32.4250946 25.64643478, 32.45020676 25.41184807, 32.45241928 25.39232445, 32.47102356 25.34958458, 32.47856522 25.33220863, 32.51583099 25.24656296, 32.64713287 25.17321396, 32.65511703 25.16874695, 31.75804329 23.47261238, 31.00323486 22.00781059, 30.97751999 22.00852013, 30.96761131 22.00839424, 30.96026993 22.00830269, 30.90060234 22.00754547, 30.84647942 22.00685692, 30.83568764 22.00671959, 30.83117867 22.00736427, 30.82583809 22.00812721, 30.82554626 22.00812531, 30.82500076 22.00811958, 30.68961716 22.00714874, 30.66023064 22.00691986, 30.63969994 22.00675964, 30.57939529 22.00671005, 30.56625938 22.00669861, 30.55422592 22.00668716, 30.53986549 22.00666618, 30.41883087 22.00647736, 30.28344536 22.00617981, 30.21887779 22.00598526, 30.14805603 22.00577164, 30.0126667 22.00526237, 30.0056839 22.00523186, 30.00000191 22.00520325, 30 22)))

what is going on here and how to fix this issue? I can of course remove the LineString in Python but it feels odd that Shapely produces this GeometryCollection in the first place.

Input polygons for shapely:

Gridcell 417

POLYGON ((30 20, 40 20, 40 30, 30 30, 30 20))

Province:

MULTIPOLYGON (((31.75804329 23.47261238, 31.00323486 22.00781059, 30.97751999 22.00852013, 30.96761131 22.00839424, 30.96026993 22.00830269, 30.90060234 22.00754547, 30.84647942 22.00685692, 30.83568764 22.00671959, 30.83117867 22.00736427, 30.82583809 22.00812721, 30.82554626 22.00812531, 30.82500076 22.00811958, 30.68961716 22.00714874, 30.66023064 22.00691986, 30.63969994 22.00675964, 30.57939529 22.00671005, 30.56625938 22.00669861, 30.55422592 22.00668716, 30.53986549 22.00666618, 30.41883087 22.00647736, 30.28344536 22.00617981, 30.21887779 22.00598526, 30.14805603 22.00577164, 30.0126667 22.00526237, 30.0056839 22.00523186, 30.00000191 22.00520325, 30 22, 30 21.99920082, 29.98200226 22, 29.87727547 22.00465012, 29.8767395 22.0046463, 29.86609077 22.00458908, 29.74188995 22.00526619, 29.69901085 22.00547409, 29.66639709 22.0056324, 29.60650444 22.00546074, 29.56506729 22.00531578, 29.47110939 22.00498962, 29.46170425 22.00494957, 29.46073151 22.00494576, 29.33572388 22.00441933, 29.32574654 22.00437546, 29.31726456 22.00433731, 29.31509018 22.0050621, 29.31305313 22.00573921, 29.31183815 22.00614357, 29.28542709 22.00580978, 29.24757767 22.00532913, 29.24201584 22.00536346, 29.20068169 22.00560951, 29.12313843 22.0060463, 29.06529045 22.00637054, 28.94155884 22.00696564, 28.92990112 22.00702095, 28.79450989 22.00757027, 28.77450752 22.00764275, 28.76717186 22.00766945, 28.75491142 22.00753403, 28.67557907 22.00665092, 28.6591301 22.00670242, 28.56130219 22.00696564, 28.52373505 22.00706673, 28.38834953 22.00733185, 28.35404587 22.00737381, 28.25295067 22.00749969, 28.12548065 22.00756073, 28.12142754 22.00720596, 28.11759949 22.00687027, 28.10774612 22.00601768, 28.10703659 22.00595665, 28.06148911 22.00563049, 28.05128288 22.00606346, 28.01960182 22.00740433, 27.99394989 22.0073719, 27.98229599 22.00735664, 27.91341209 22.00723648, 27.84690094 22.00712013, 27.80693817 22.00702858, 27.79823494 22.00700951, 27.76398659 22.00814247, 27.73838615 22.00764847, 27.71154594 22.00712967, 27.68975639 22.00670433, 27.67487144 22.00641251, 27.63397217 22.00765228, 27.57618523 22.00754547, 27.56340981 22.00751495, 27.44078827 22.00723076, 27.30685806 22.00680351, 27.30540276 22.0067997, 27.3038044 22.00679398, 27.27376366 22.00667953, 27.17000961 22.0062809, 27.03462982 22.00564957, 26.94745064 22.0051918, 26.9069519 22.00542068, 26.89923859 22.00546455, 26.76384926 22.00617027, 26.63140106 22.00676346, 26.62846184 22.00677681, 26.53366661 22.00713921, 26.49551773 22.00611115, 26.49487686 22.0061512, 26.49308968 22.00625992, 26.49061966 22.00641251, 26.47603035 22.00731087, 26.39429855 22.00611305, 26.3770237 22.00586128, 26.37558937 22.0058403, 26.37342262 22.0058403, 26.35775185 22.00584984, 26.22236061 22.0058403, 26.20599556 22.00582695, 26.08696938 22.00573158, 26.05418015 22.00568962, 25.96693993 22.00576591, 25.95157814 22.00577927, 25.88501167 22.00579453, 25.81619263 22.00580978, 25.68079948 22.00573921, 25.54541016 22.00556946, 25.51301765 22.0055027, 25.41002083 22.00528908, 25.27463531 22.00491714, 25.27050781 22.00490379, 25.1392498 22.00444031, 25.02123833 22, 24.99999809 21.99920082, 24.99999809 22.0658474, 24.99999809 22.08250046, 24.99999809 22.13006401, 24.99999809 22.16579819, 24.99999809 22.24909973, 24.99999809 22.33239937, 24.99999809 22.38243866, 24.99999809 22.41570282, 24.99999809 22.49899864, 24.99999809 22.58230019, 24.99999809 22.634758, 24.99999809 22.66559982, 24.99999809 22.74889946, 24.99999809 22.8321991, 24.99999809 22.8870182, 24.99999809 22.91550064, 24.99999809 22.99880028, 24.99999809 23.01312447, 24.99999809 23.08209991, 24.99999809 23.13921356, 24.99999809 23.16540146, 24.99999809 23.24869728, 24.99999809 23.26527977, 24.99999809 23.33200073, 24.99999809 23.39133453, 24.99999809 23.41530037, 24.99999809 23.49860191, 24.99999809 23.51736641, 24.99999809 23.5357933, 24.99999809 23.58189964, 24.99999809 23.64338493, 24.99999809 23.66519928, 24.99999809 23.74850082, 24.99999809 23.76938248, 24.99999809 23.83180237, 24.99999809 23.89536095, 24.99999809 23.91509819, 24.99999809 23.99839973, 24.99999809 24.08169937, 24.99999809 24.11018181, 24.99999809 24.16500092, 24.99999809 24.24830055, 24.99999809 24.36219406, 24.99999809 24.47298431, 24.99999809 24.49819946, 24.99999809 24.58149719, 24.99999809 24.61426544, 24.99999809 24.66480064, 24.99999809 24.69965363, 24.99999809 24.74810028, 24.99999809 24.83140182, 24.99999809 24.91469955, 24.99999809 25, 24.99999809 25.08329964, 24.99999809 25.11854935, 24.99999809 25.16659927, 24.99999809 25.24990082, 24.99999809 25.33320045, 24.99999809 25.37077141, 24.99999809 25.41650009, 24.99999809 25.49980164, 24.99999809 25.58309937, 24.99999809 25.66640091, 24.99999809 25.74970055, 24.99999809 25.83300018, 24.99999809 25.91629982, 24.99999809 25.99959946, 24.99999809 26.08289909, 24.99999809 26.16620255, 24.99999809 26.24950027, 24.99999809 26.33279991, 24.99999809 26.41609955, 24.99999809 26.49939919, 24.99999809 26.58270073, 24.99999809 26.66392326, 24.99999809 26.66600037, 24.99999809 26.7493, 24.99999809 26.83259964, 24.99999809 26.91589737, 24.99999809 26.97445869, 24.99999809 26.99920082, 24.99999809 27.00765419, 25.00005531 27.00765419, 25.00005531 27.04971123, 25.00005531 27.08250046, 25.00005531 27.16580009, 25.00005531 27.24909973, 25.00005531 27.33239937, 25.00005531 27.38955688, 25.00005531 27.41570091, 25.00005531 27.45453262, 25.00005531 27.49900055, 25.00005531 27.51569748, 25.00005531 27.58229828, 25.00005531 27.64191437, 25.00005531 27.65014839, 27.36088943 27.67578125, 28.25244141 27.67483521, 28.46664238 27.67460823, 28.53467369 27.67272758, 30.53003883 27.6014328, 30.86449623 26.79118919, 31.15589523 26.82732773, 31.29156494 26.83351135, 31.44809151 26.63268852, 31.76193428 26.29763031, 31.95477867 26.13333893, 32.03946686 26.06095505, 32.21749496 25.97387695, 32.24517059 25.98464203, 32.32149506 25.98672104, 32.44264603 26.04069328, 32.67482376 26.10331726, 32.73316956 26.07575607, 32.73933411 26.05672836, 32.72978973 25.98652077, 32.72766876 25.98363495, 32.72722244 25.97396469, 32.68310928 25.91794777, 32.66844177 25.84004784, 32.64354324 25.79011726, 32.53351593 25.69179344, 32.4250946 25.64643478, 32.45020676 25.41184807, 32.45241928 25.39232445, 32.47102356 25.34958458, 32.47856522 25.33220863, 32.51583099 25.24656296, 32.64713287 25.17321396, 32.65511703 25.16874695, 31.75804329 23.47261238)))

Update:

It seems that the south-western part of gridcell 417 causes the issue. The gridcell is the same as the province boundary. Vertices 25, 26, 27 all have x coordinates of 30.0000

enter link description here

Update 2:

I fixed the issue with applying a buffer(0) in Shapely. No solution for QGIS yet.

Note that this function can produce empty geometry. You can remove the empty geometries using:

gdf_nonempty = gdf.loc[gdf.geometry.is_empty == False]
  • You should make a bug report for the issue in QGIS: www.issues.qgis.org – csk Dec 5 '18 at 19:59
2

The geometry of province has 6 millimeters long section that (almost) follows the boundary of polygon #417.

enter image description here

You can find that place from this corner

enter image description here

But the linear intersection is just a part of the common boundary. The rest of it forms a spike into the polygon part. Spike is not a topology error because boundary does not self-intersect but it is perhaps not what you want either.

It is sometimes hard to deside how to fix this kind of topology issues. Simplifying the source data by tolerance of 1 cm might resolve this case.

enter image description here

  • Instead of spitting out an empty shapefile, QGIS could at least capture the error or provide the option to keep the polygon part only. – RutgerH Dec 3 '18 at 12:47
  • 1
    As far as I know QGIS does not support GeometryCollections and it will have difficult times when processing algorithms create those. In your case the result from Shapely is correct and you will get the same result with software like PostGIS, SpatiaLite, or OpenJUMP. Perhaps QGIS should have a special tool "Intersection area" that would create only polygons. While such tool does not exist you should perhaps try OpenJUMP which deals well with GeometryCollections. – user30184 Dec 3 '18 at 13:13

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.