Intersection of polygons gives incorrect results, depending on order of arguments

[user2699]

Picture shows it, code to reproduce below. Both geometries are valid. The adjustment is setting precision to a very small value on one geometry. The only difference between sub-plots is the order used in calculating the intersection in the second plot, and a tiny rounding difference in the third.

Code to reproduce:

import shapely
from matplotlib import pyplot as plt
from shapely.plotting import plot_polygon, plot_line

# Define shapes
shape1 = shapely.from_wkt(
    'POLYGON ((0.3081742230052605 0.8574917215741684, 0.3235558676951882 0.959764955575765, 0.3268384579963268 0.9761955555650583, 0.3309851357533714 0.9899402381062132, 0.3364277995673436 1.0032251410108646, 0.3431165154642001 1.015928381315176, 0.3509899175067092 1.0279334125460768, 0.3599757707991199 1.039130093982343, 0.3699916342095888 1.0494157011457903, 0.3809456167302097 1.058695868251819, 0.392737220535371 1.0668854539727242, 0.4052582630037556 1.073909322570786, 0.4183938692448725 1.079703033234622, 0.4320235260241182 1.0842134312944256, 0.4460221874170914 1.0873991358919737, 0.4602614220492758 1.0892309196312477, 0.4746105913957159 1.0896919767265496, 0.4913319068024928 1.0886253811328888, 1.6700718959842467 0.9565165637600046, 1.57979867649902 0.6026736686526396, 1.091045353275972 0.0772874773100065, 1.0880332290753003 0.0751954761072001, 1.0807249769395977 0.0690039561787506, 1.0740426168917647 0.0621416291321494, 1.0680474565569975 0.0546714537082506, 1.063366321369268 0.0475338562857705, 1.0334449647384247 0.0153698445990219, 0.7231700686154534 0.6266171384924346, 0.7149755148483247 0.6410255508351703, 0.7068170547053731 0.6526513145597219, 0.697594844851672 0.6634526894489446, 0.6873916890776646 0.6733326928059559, 0.6762991988285427 0.6822026146859029, 0.6644169706492398 0.6899828143991633, 0.6497523731524246 0.6977095626115124, 0.3081742230052605 0.8574917215741684))',
)
shape2 = shapely.from_wkt(
'POLYGON ((1.108561208420377086625308038492221385240554809570312500000000 0.096116174702266674034767390821798471733927726745605468750000, 1.579798676499019949304170040704775601625442504882812500000000 0.602673668652639604736975798005005344748497009277343750000000, 1.445650745099990475139861700881738215684890747070312500000000 0.076855572327130111665027811795880552381277084350585937500000, 1.444900805558562373320796723419334739446640014648437500000000 0.076572789294606147247890248763724230229854583740234375000000, 1.429693401346880499147573573281988501548767089843750000000000 0.072597961830626389634346651291707530617713928222656250000000, 1.414153318559506544715986819937825202941894531250000000000000 0.070237829803077911350328577100299298763275146484375000000000, 1.398451511191880936024745096801780164241790771484375000000000 0.069518356651499413345618449966423213481903076171875000000000, 1.380137763162443853559580020373687148094177246093750000000000 0.070602770087456412362314495112514123320579528808593750000000, 1.160083364028043684257340828480664640665054321289062500000000 0.095265653076203160587454021879239007830619812011718750000000, 1.150524405216126089257500098028685897588729858398437500000000 0.095875386407627882778825778586906380951404571533203125000000, 1.140950960368624311769281121087260544300079345703125000000000 0.095567779425867860343402071521268226206302642822265625000000, 1.131450861503619176318125028046779334545135498046875000000000 0.094345654285582614129168632643995806574821472167968750000000, 1.122111267723040706201231841987464576959609985351562500000000 0.092220223431104253464241082838270813226699829101562500000000, 1.113017865567405006999024408287368714809417724609375000000000 0.089210986727341734514595827931771054863929748535156250000000, 1.112913182538546097433140857901889830827713012695312500000000 0.089164814272908357595248673987953225150704383850097656250000, 1.108561208420377086625308038492221385240554809570312500000000 0.096116174702266674034767390821798471733927726745605468750000))'
)

# Define plotting function
def plot_intersection(ax, A, B):
    plot_polygon(A, ax=ax, facecolor='none', edgecolor='r', hatch="/")
    plot_polygon(B, ax=ax, facecolor='none', edgecolor='b', hatch="\\")
    plot_polygon(shapely.intersection(A, B), ax=ax, linewidth=5, linestyle="-.", edgecolor='g', alpha=0.5, hatch="||||")
    return ax

# Make sure there's nothing wrong with the geometry
assert shape1.is_valid
assert shape2.is_valid

# And do some clean up just in case...
shape1 = shape1.buffer(0)
shape2 = shape2.buffer(0)
shape1 = shape1.normalize()
shape2 = shape2.normalize()

# Plot outline of shapes and shade their intersection
fig, [ax1, ax2, ax3] = plt.subplots(1, 3, sharex=True, sharey=True)

ax1.set_title("Wrong behavior")
plot_intersection(ax1, shape1, shape2)
ax2.set_title("Expected behavior")
plot_intersection(ax2, shape2, shape1)
ax3.set_title("Behavior, after \"adjusting\" geometry")
plot_intersection(ax3, shape1, shapely.set_precision(shape2, 1e-300))
plt.show()

Version info:

Python 3.13.7 (main, Mar  3 2026, 12:19:54) [GCC 15.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import shapely
>>> shapely.__version__
'2.1.2'
>>> shapely.geos_version_string
'3.13.1'

On Ubuntu 25.10