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

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

# Define the coordinates for the LineString
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)

nearest1 = nearest_points(intersection, 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.