How to find point of intersection of two line segments in Python?
Question: I have data with one independent variable x and two dependent variables y1 and
y2 as shown below:

    
    
     x      y1      y2
    -1.5    16.25   1.02
    -1.25   17      1.03
    -1      15      1.03
    -0.75   9       1.09
    -0.5    5.9     1.15
    -0.25   5.2     1.17
     0      4.77    1.19
    +0.25   3.14    1.35
    +0.5    2.5     1.54
    +0.75   2.21    1.69
    +1      1.91    1.96
    +1.25   1.64    2.27
    +1.5    1.52    2.56
    +1.75   1.37    3.06
    +2      1.24    4.12
    +2.25   1.2     4.44
    +2.5    1.18    4.95
    +2.75   1.12    6.49
    +3      1.07    10
    

So, here the value of `x where y1 = y2` is somewhere around `+1`. How do I
read the data and calculate this in python?


Answer: The naive solution goes like this:

    
    
    txt = """-1.5    16.25   1.02
    -1.25   17      1.03
    -1      15      1.03
    -0.75   9       1.09
    -0.5    5.9     1.15
    -0.25   5.2     1.17
     0      4.77    1.19
    +0.25   3.14    1.35
    +0.5    2.5     1.54
    +0.75   2.21    1.69
    +1      1.91    1.96
    +1.25   1.64    2.27
    +1.5    1.52    2.56
    +1.75   1.37    3.06
    +2      1.24    4.12
    +2.25   1.2     4.44
    +2.5    1.18    4.95
    +2.75   1.12    6.49
    +3      1.07    10"""
    
    import numpy as np
    # StringIO behaves like a file object, use it to simulate reading from a file
    from StringIO import StringIO   
    
    x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
    p1 = np.poly1d(np.polyfit(x, y1, 1))
    p2 = np.poly1d(np.polyfit(x, y2, 1))
    print 'equations: ',p1,p2
    
    #y1 and y2 have to be equal for some x, that you solve for :
    #  a x+ b = c x + d  --> (a-c) x= d- b
    
    a,b=list(p1)
    c,d=list(p2)
    x=(d-b)/(a-c)
    print 'solution x= ',x
    

output:

    
    
    equations:   
    -3.222 x + 7.323  
    1.409 x + 1.686
    solution x=  1.21717324767
    

But then you plot the 'lines':

    
    
    import matplotlib.pyplot as p
    %matplotlib inline
    p.plot(x,y1,'.-')
    p.plot(x,y2,'.-')
    

[![enter image description
here](http://i.stack.imgur.com/77clE.png)](http://i.stack.imgur.com/77clE.png)

And you realize you can't use a linear assumption but for a few segments.

    
    
    x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
    x,y1,y2=x[8:13],y1[8:13],y2[8:13]
    p1 = np.poly1d(np.polyfit(x, y1, 1))
    p2 = np.poly1d(np.polyfit(x, y2, 1))
    print 'equations: ',p1,p2
    a,b=list(p1)
    c,d=list(p2)
    x0=(d-b)/(a-c)
    print 'solution x= ',x0
    p.plot(x,y1,'.-')
    p.plot(x,y2,'.-')
    

Output:

    
    
    equations:   
    -1.012 x + 2.968  
    1.048 x + 0.956
    solution x=  0.976699029126
    

[![enter image description
here](http://i.stack.imgur.com/TBmK1.png)](http://i.stack.imgur.com/TBmK1.png)

Even now one could improve by leaving two more points out (looking very
linear, but that can be coincidental for a few points).

    
    
    x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
    x1,x2=x[8:12],x[9:13]
    y1,y2=y1[8:12],y2[9:13]
    p1 = np.poly1d(np.polyfit(x1, y1, 1))
    p2 = np.poly1d(np.polyfit(x2, y2, 1))
    print 'equations: ',p1,p2
    a,b=list(p1)
    c,d=list(p2)
    x0=(d-b)/(a-c)
    print 'solution x= ',x0
    
    import matplotlib.pyplot as p
    %matplotlib inline
    p.plot(x1,y1,'.-')
    p.plot(x2,y2,'.-')
    

Output:

    
    
    equations:   
    -1.152 x + 3.073  
    1.168 x + 0.806
    solution x=  0.977155172414
    

[![enter image description
here](http://i.stack.imgur.com/9ml2F.png)](http://i.stack.imgur.com/9ml2F.png)

Possibly better would be to use more points and apply a 2nd order
interpolation `np.poly1d(np.polyfit(x,y1,2))` and then solve the equality for
two 2nd order polynomials, which I leave as an exercise (quadratic equation)
for the reader.

