shapely: Intersection gives wrong result

Expected behavior and actual behavior.

I ran into a case where an intersection between two (valid) (multi)polygons results in a wrong result. One example results in an “invalid” wrong result, another one results in a “valid” wrong result.

Probably the problem is in geos… but because my test case is using shapely I thought it was more logical to post it here.

Polygon 1 in the example below is the culprit:

Steps to reproduce the problem.

import matplotlib.pyplot as plt
import geopandas as gpd
import shapely
import shapely.wkt

# Calculate intersection
print(shapely.geos.geos_version_string)
wkt1_str = "MultiPolygon (((66697.40120137332996819 185279.95469107336248271, 66698.375 185273.625, 66697.375 185280.125, 66697.40120137332996819 185279.95469107336248271)),((66705.875 185229.625, 66704.875 185234.625, 66703.875 185240.625, 66702.875 185245.375, 66706.90145934304746334 185229.45392344283754937, 66710.375 185228.875, 66705.875 185229.625)))"
wkt2_str = "MultiPolygon (((66720 185280, 66720 185200, 66680 185200, 66680 185280, 66720 185280)))"
wkt3_str = "MultiPolygon (((66720 185320, 66720 185280, 66680 185280, 66680 185320, 66720 185320)))"
poly1 = shapely.wkt.loads(wkt1_str)
print(f"poly1.is_valid: {poly1.is_valid}")
poly2 = shapely.wkt.loads(wkt2_str)
print(f"poly2.is_valid: {poly2.is_valid}")
poly3 = shapely.wkt.loads(wkt3_str)
print(f"poly3.is_valid: {poly3.is_valid}")
poly1_intersect_poly2 = poly1.intersection(poly2)
print(f"poly1_intersect_poly2.is_valid: {poly1_intersect_poly2.is_valid}")
poly1_intersect_poly3 = poly1.intersection(poly3)
print(f"poly1_intersect_poly3.is_valid: {poly1_intersect_poly3.is_valid}")

# Plot
print(poly1_intersect_poly2.wkt)
fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True)
gpd.GeoSeries([poly1]).boundary.plot(ax=ax[0][0], edgecolor='blue')
ax[0][0].set_title('polygon 1')
ax[0][0].ticklabel_format(useOffset=False)
ax[0][0].tick_params(axis='x', rotation=45)
gpd.GeoSeries([poly2]).boundary.plot(ax=ax[0][1], edgecolor='green')
ax[0][1].set_title('polygon 2')
ax[0][1].tick_params(axis='x', rotation=45)
gpd.GeoSeries([poly1_intersect_poly2]).boundary.plot(ax=ax[0][2], edgecolor='red')
ax[0][2].set_title('intersection poly1 + poly2')
ax[0][2].tick_params(axis='x', rotation=45)

gpd.GeoSeries([poly1]).boundary.plot(ax=ax[1][0], edgecolor='blue')
ax[1][0].set_title('polygon 1')
ax[1][0].ticklabel_format(useOffset=False)
ax[1][0].tick_params(axis='x', rotation=45)
gpd.GeoSeries([poly3]).boundary.plot(ax=ax[1][1], edgecolor='green')
ax[1][1].set_title('polygon 3')
ax[1][1].tick_params(axis='x', rotation=45)
gpd.GeoSeries([poly1_intersect_poly3]).boundary.plot(ax=ax[1][2], edgecolor='red')
ax[1][2].set_title('intersection poly1 + poly3')
ax[1][2].tick_params(axis='x', rotation=45)
plt.show()

Operating system

Windows 10

Shapely version and provenance

shapely 1.8.0 from conda-forge geos 3.10.0 from conda-forge