Intersection point doesn't intersect either of the lines

Question

I have a road network and I need to split two roads (lines) at their intersection. However, when finding the intersection, this point doesn't intersect with either of the lines. The distance between the intersection and the lines is not 0. Here's a reproductible example:

### Example of shapely and floating points, and how it doesn't work: 

from shapely.geometry import Point, LineString
from shapely.ops import nearest_points
from shapely.ops import snap

# Define the coordinates for the LineString

coords1 = [
    (73.0867321, 33.6644942),
    (73.08678639972806, 33.66455575008311),
    (73.08684069953381, 33.66461730014245),
    (73.08689499941724, 33.66467885017808),
    (73.08694929937838, 33.66474040018994),
    (73.08700359941724, 33.66480195017806),
    (73.08705789953379, 33.664863500142445),
    (73.08711219972804, 33.6649250500831),
    (73.0871665, 33.6649866)
]

coords2 = [
    (73.08690239999999, 33.6650841), (73.08694961275046, 33.665018912562815), 
    (73.08699682542931, 33.664953725107694), (73.08704403803665, 33.664888537634624), 
    (73.08709125057243, 33.66482335014359), (73.08713846303664, 33.66475816263463), 
    (73.08718567542934, 33.6646929751077), (73.08723288775042, 33.66462778756282), 
    (73.08728009999999, 33.6645626), (73.0873316, 33.664484), (73.08737451271202, 33.66442327505189),
    (73.0899173, 33.6613607), (73.08996137509759, 33.661297200023455), 
    (73.09000545013011, 33.66123370003129), (73.09004952509757, 33.66117020002346), 
    (73.0900936, 33.6611067), (73.09013767509688, 33.661043400023274), 
    (73.09018175012917, 33.660980100031026), (73.0902258250969, 33.66091680002327), 
    (73.09026989999998, 33.6608535), (73.09031404732161, 33.66081331574375), 
    (73.09035819461442, 33.66077313147986), (73.0904023418784, 33.66073294720829),
    (73.09157271456951, 33.65874581072471), (73.09160897388341, 33.65869747008384),
    (73.09164523315653, 33.658649129432376), (73.09168149238894, 33.658600788770315),
    (73.09171775158057, 33.65855244809768), (73.0917540107315, 33.658504107414444),
    (73.09179026984168, 33.658455766720635), (73.09182652891113, 33.658407426016254),
    (73.09186278793982, 33.65835908530127), (73.09189904692778, 33.658310744575715),
    (73.09193530587498, 33.65826240383956), (73.09197156478149, 33.65821406309283),
    (73.09200782364721, 33.6581657223355), (73.09204408247223, 33.65811738156759),
    (73.09208034125646, 33.65806904078909), (73.0921166, 33.6580207)
]

# Create the LineString object
line2 = LineString(coords2)

line1 = LineString(coords1)

print(line1.intersects(line2))
intersection  = line1.intersection(line2)

print(intersect.intersects(line1))

print(line2.project(intersection))

p = line2.interpolate(line2.project(intersection))
p.intersects(line2)

snap_precision =  0.0001
snap_line1 = snap(intersection, line1, snap_precision)
snap_line2 = snap(intersection, line2, snap_precision)

 
np0, np1 = nearest_points(intersection, line1)
print(np1.intersects(line1))

I tried finding the nearest point to the intersection within the line, but doesn't seem to work . I found this relevant post where line1.interpolate(line1.project(pt)) is suggested, however this doesn't work either.

Answer 1

That's a floating's accuracy issue.

To save yourself the trouble, I would suggest doing a unary_union and getting its geoms :

from shapely import unary_union

uu = unary_union(list(map(LineString, [coords1, coords2])))

lines = [l for l in uu.geoms] # << `uu` is a MultiLineString

Output (lines) :

[
    <LINESTRING (73.087 33.664, 73.087 33.665, 73.087 33.665, 73.087 33.665, 73....>,
    <LINESTRING (73.087 33.665, 73.087 33.665, 73.087 33.665)>,
    <LINESTRING (73.087 33.665, 73.087 33.665, 73.087 33.665, 73.087 33.665, 73....>,
    <LINESTRING (73.087 33.665, 73.087 33.665, 73.087 33.665, 73.087 33.665, 73....>
]