import numpy as np

def get_intersect(a1, a2, b1, b2):
    """ 
    Returns the point of intersection of the lines passing through a2,a1 and b2,b1.
    a1: [x, y] a point on the first line
    a2: [x, y] another point on the first line
    b1: [x, y] a point on the second line
    b2: [x, y] another point on the second line
    """
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

if __name__ == "__main__":
    print get_intersect((0, 1), (0, 2), (1, 10), (1, 9))  # parallel  lines
    print get_intersect((0, 1), (0, 2), (1, 10), (2, 10)) # vertical and horizontal lines
    print get_intersect((0, 1), (1, 2), (0, 10), (1, 9))  # another line for fun

Explanation

Note that the equation of a line is ax+by+c=0. So if a point is on this line, then it is a solution to (a,b,c).(x,y,1)=0 (. is the dot product)

let l1=(a1,b1,c1), l2=(a2,b2,c2) be two lines and p1=(x1,y1,1), p2=(x2,y2,1) be two points.

Finding the line passing through two points:

let t=p1xp2 (the cross product of two points) be a vector representing a line.

We know that p1 is on the line t because t.p1 = (p1xp2).p1=0. We also know that p2 is on t because t.p2 = (p1xp2).p2=0. So t must be the line passing through p1 and p2.

This means that we can get the vector representation of a line by taking the cross product of two points on that line.

Finding the point of intersection:

Now let r=l1xl2 (the cross product of two lines) be a vector representing a point

We know r lies on l1 because r.l1=(l1xl2).l1=0. We also know r lies on l2 because r.l2=(l1xl2).l2=0. So r must be the point of intersection of the lines l1 and l2.

Interestingly, we can find the point of intersection by taking the cross product of two lines.

Can't stand aside,

So we have linear system:

A1 * x + B1 * y = C1
A2 * x + B2 * y = C2

let's do it with Cramer's rule, so solution can be found in determinants:

x = Dx/D
y = Dy/D

where D is main determinant of the system:

A1 B1
A2 B2

and Dx and Dy can be found from matricies:

C1 B1
C2 B2

and

A1 C1
A2 C2

(notice, as C column consequently substitues the coef. columns of x and y)

So now the python, for clarity for us, to not mess things up let's do mapping between math and python. We will use array L for storing our coefs A, B, C of the line equations and intestead of pretty x, y we'll have [0], [1], but anyway. Thus, what I wrote above will have the following form further in the code:

for D

L1[0] L1[1]
L2[0] L2[1]

for Dx

L1[2] L1[1]
L2[2] L2[1]

for Dy

L1[0] L1[2]
L2[0] L2[2]

Now go for coding:

line - produces coefs A, B, C of line equation by two points provided,
intersection - finds intersection point (if any) of two lines provided by coefs.

from __future__ import division 

def line(p1, p2):
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C

def intersection(L1, L2):
    D  = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x,y
    else:
        return False

Usage example:

L1 = line([0,1], [2,3])
L2 = line([2,3], [0,4])

R = intersection(L1, L2)
if R:
    print "Intersection detected:", R
else:
    print "No single intersection point detected"

Shapely has some cool functions. According to this post, this code should work:

from shapely.geometry import LineString
from shapely.geometry import Point

p = Point(0,0)//center
c = p.buffer(0.71).boundary//radius
l = LineString([(0.,0.), (1., 1.)])//line point
i = c.intersection(l)

Apparently here i is the array you are looking for, also, check this post too. Hope this helps.

Here is a function I wrote to find the closest point between two 3d lines

import scipy.optimize
#takes in two lines, the line formed by pt1 and pt2, and the line formed by pt3 and pt4, and finds their intersection or closest point
def fourptsMeetat(pt1,pt2,pt3,pt4):
    #least squares method
    def errFunc(estimates):
        s, t = estimates
        x = pt1 + s * (pt2 - pt1) - (pt3 + t * (pt4 - pt3))
        return x

    estimates = [1, 1]

    sols = scipy.optimize.least_squares(errFunc, estimates)
    s,t = sols.x

    x1 =  pt1[0] + s * (pt2[0] - pt1[0])
    x2 =  pt3[0] + t * (pt4[0] - pt3[0])
    y1 =  pt1[1] + s * (pt2[1] - pt1[1])
    y2 =  pt3[1] + t * (pt4[1] - pt3[1])
    z1 =  pt1[2] + s * (pt2[2] - pt1[2])
    z2 = pt3[2] + t * (pt4[2] - pt3[2])

    x = (x1 + x2) / 2  #halfway point if they don't match
    y = (y1 + y2) / 2  # halfway point if they don't match
    z = (z1 + z2) / 2  # halfway point if they don't match

    return (x,y,z)