Using matplotlib to solve Point in Polygone

Question

I am looking for an algorithm to check if a point is within a polygon or not.

I am currently using mplPath and contains_point() but it doesn't seem to work in some cases.

EDIT 16 Sept 2016:

Okay so I imporved my code by simply checking if the point where also on the edges. I still have some issues with the rectangle and bow-tie example though:

NEW CODE:

#for PIP problem
import matplotlib.path as mplPath
import numpy as np
#for plot
import matplotlib.pyplot as plt


def plot(poly,points):
    bbPath = mplPath.Path(poly)
    #plot polygon
    plt.plot(*zip(*poly))

    #plot points
    xs,ys,cs = [],[],[]
    for point in points:
        xs.append(point[0])
        ys.append(point[1])
        color = inPoly(poly,point)
        cs.append(color)
        print point,":", color
    plt.scatter(xs,ys, c = cs , s = 20*4*2)

    #setting limits
    axes = plt.gca()
    axes.set_xlim([min(xs)-5,max(xs)+50])
    axes.set_ylim([min(ys)-5,max(ys)+10])

    plt.show()

def isBetween(a, b, c): #is c between a and b ?
    crossproduct = (c[1] - a[1]) * (b[0] - a[0]) - (c[0] - a[0]) * (b[1] - a[1])
    if abs(crossproduct) > 0.01 : return False   # (or != 0 if using integers)

    dotproduct = (c[0] - a[0]) * (b[0] - a[0]) + (c[1] - a[1])*(b[1] - a[1])
    if dotproduct < 0 : return False

    squaredlengthba = (b[0] - a[0])*(b[0] - a[0]) + (b[1] - a[1])*(b[1] - a[1])
    if dotproduct > squaredlengthba: return False

    return True

def get_edges(poly):
    # get edges
    edges = []
    for i in range(len(poly)-1):
        t = [poly[i],poly[i+1]]
        edges.append(t)
    return edges

def inPoly(poly,point):
    if bbPath.contains_point(point) == True:
        return 1
    else:
        for e in get_edges(poly):
            if isBetween(e[0],e[1],point):
                return 1
    return 0
# TESTS ========================================================================
#set up poly
polys = {
1 : [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10],[10,10]], # test rectangulary shape
2 : [[20,10],[10,20],[30,20],[20,10]], # test triangle
3 : [[0,0],[0,10],[20,0],[20,10],[0,0]], # test bow-tie
4 : [[0,0],[0,10],[20,10],[20,0],[0,0]] # test rect
}

#points to check
points = {
1 : [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)], # rectangulary shape test pts
2 : [[20,10],[10,20],[30,20],[-5,0],[20,15]] , # triangle  test pts
3 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]],  # bow-tie shape test pts
4 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]]  # rect shape test pts
}

#print bbPath.contains_points(points) #0 if outside, 1 if inside
for data in zip(polys.itervalues(),points.itervalues()):
    plot(data[0],data[1])

Outputs from new code:

OLD CODE:

#for PIP problem
import matplotlib.path as mplPath
import numpy as np
#for plot
import matplotlib.pyplot as plt

#set up poly
array = np.array([[10,10],[10,50],[50,50],[50,80],[100,80],[100,10]])
bbPath = mplPath.Path(array)

#points to check
points = [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)]


print bbPath.contains_points(points) #0 if outside, 1 if inside

#plot polygon
plt.plot(*zip(*array))

#plot points
xs,ys,cs = [],[],[]
for point in points:
    xs.append(point[0])
    ys.append(point[1])
    cs.append(bbPath.contains_point(point))
plt.scatter(xs,ys, c = cs)

#setting limits
axes = plt.gca()
axes.set_xlim([0,120])
axes.set_ylim([0,100])

plt.show()

I come up with the following . As you can see the three points surrounded in red are indicated as being outside the polygon (in blue) when I would expect them to be inside.

I also tried changing the radius value of the path bbPath.contains_points(points, radius = 1.) but that didn't make any difference.

Any help would be welcome.

EDIT :

screenshots from the algorithm proposed in the answer to this question seem to show that it fails for other cases.

Answer 1

Okay so I finally managed to get it done using shapely instead.

#for PIP problem
import matplotlib.path as mplPath
import numpy as np
#for plot
import matplotlib.pyplot as plt
import shapely.geometry as shapely

class MyPoly(shapely.Polygon):
    def __init__(self,points):
        super(MyPoly,self).__init__(points)
        self.points = points
        self.points_shapely = [shapely.Point(p[0],p[1]) for p in points]

def convert_to_shapely_points_and_poly(poly,points):
    poly_shapely = MyPoly(poly)
    points_shapely = [shapely.Point(p[0],p[1]) for p in points]
    return poly_shapely,points_shapely

def plot(poly_init,points_init):
    #convert to shapely poly and points
    poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)

    #plot polygon
    plt.plot(*zip(*poly.points))

    #plot points
    xs,ys,cs = [],[],[]
    for point in points:
        xs.append(point.x)
        ys.append(point.y)
        color = inPoly(poly,point)
        cs.append(color)
        print point,":", color
    plt.scatter(xs,ys, c = cs , s = 20*4*2)

    #setting limits
    axes = plt.gca()
    axes.set_xlim([min(xs)-5,max(xs)+50])
    axes.set_ylim([min(ys)-5,max(ys)+10])

    plt.show()


def isBetween(a, b, c): #is c between a and b ?
    crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
    if abs(crossproduct) > 0.01 : return False   # (or != 0 if using integers)

    dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y)
    if dotproduct < 0 : return False

    squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
    if dotproduct > squaredlengthba: return False

    return True

def get_edges(poly):
    # get edges
    edges = []
    for i in range(len(poly.points)-1):
        t = [poly.points_shapely[i],poly.points_shapely[i+1]]
        edges.append(t)
    return edges


def inPoly(poly,point):
    if poly.contains(point) == True:
        return 1
    else:
        for e in get_edges(poly):
            if isBetween(e[0],e[1],point):
                return 1
    return 0


# TESTS ========================================================================
#set up poly
polys = {
1 : [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10],[10,10]], # test rectangulary shape
2 : [[20,10],[10,20],[30,20],[20,10]], # test triangle
3 : [[0,0],[0,10],[20,0],[20,10],[0,0]], # test bow-tie
4 : [[0,0],[0,10],[20,10],[20,0],[0,0]], # test rect clockwise
5 : [[0,0],[20,0],[20,10],[0,10],[0,0]] # test rect counter-clockwise
}

#points to check
points = {
1 : [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)], # rectangulary shape test pts
2 : [[20,10],[10,20],[30,20],[-5,0],[20,15]] , # triangle  test pts
3 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]],  # bow-tie shape test pts
4 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]],  # rect shape test pts
5 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]]  # rect shape test pts
}

#print bbPath.contains_points(points) #0 if outside, 1 if inside
for data in zip(polys.itervalues(),points.itervalues()):
    plot(data[0],data[1])