4

I am trying to detect if a given point(x,y) is in a polygon of n*2 array. But it seems that some points on the borders of the polygon return that it's not include.

def point_inside_polygon(x,y,poly):

    n = len(poly)
    inside =False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xinters = (y-p1y)*(p2x-p1x)/float((p2y-p1y))+p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside
Georgy
  • 12,464
  • 7
  • 65
  • 73
Rabih Assaf
  • 95
  • 1
  • 1
  • 8
  • 1
    are those coordinates integers or floating? and python 2 or python 3? – Jean-François Fabre Sep 23 '16 at 12:36
  • float and python 2. i also changed the division to float and it gives false on some points on the boundary – Rabih Assaf Sep 23 '16 at 14:24
  • if it's float, then it's bad to compare `p1x == p2x`: it could be equal or not, there's a precision loss problem. – Jean-François Fabre Sep 23 '16 at 14:27
  • Does this answer your question? [Checking if a point is inside a polygon](https://stackoverflow.com/questions/16625507/checking-if-a-point-is-inside-a-polygon) – Georgy Mar 04 '20 at 09:14
  • Does this answer your question? [How can I determine whether a 2D Point is within a Polygon?](https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon) – Fabrizio Mar 04 '20 at 13:01

3 Answers3

7

Here is a simple way with multiple options

from shapely.geometry import Point, Polygon

# Point objects(Geo-coordinates)
p1 = Point(24.952242, 60.1696017)
p2 = Point(24.976567, 60.1612500)

# Polygon
coords = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)]
poly = Polygon(coords)

Method 1:

Using within Function:

Syntax: point.within(polygon)

# Check if p1 is within the polygon using the within function
print(p1.within(poly))
# True

# Check if p2 is within the polygon
print(p2.within(poly))
# False

Method 2:

Using contains Function:

Syntax: polygon.contains(point)

# Check if polygon contains p1 
print(poly.contains(p1))
# True

# Check if polygon contains p2
print(poly.contains(p2))
# False

To check if point is lying on the boundary of Polygon:

Using touches Function:

Syntax: polygon.touches(point)

poly1 = Polygon([(0, 0), (1, 0), (1, 1)])
point1 = Point(0, 0)

poly1.touches(point1)
# True

If you want to speed up the process then please use

import shapely.speedups
shapely.speedups.enable()

and

Make use of Geopandas


Reference:

  1. https://automating-gis-processes.github.io/CSC18/lessons/L4/point-in-polygon.html#point-in-polygon-using-geopandas

  2. https://shapely.readthedocs.io/en/latest/manual.html

  3. https://streamhacker.com/2010/03/23/python-point-in-polygon-shapely/

Suhas_Pote
  • 3,620
  • 1
  • 23
  • 38
2

You may use contains_point function from matplotlib.path with small negative and positive radius (small trick). Something like this:

import matplotlib.path as mplPath
import numpy as np

crd = np.array([[0,0], [0,1], [1,1], [1,0]])# poly
bbPath = mplPath.Path(crd)
pnts = [[0.0, 0.0],[1,1],[0.0,0.5],[0.5,0.0]] # points on edges
r = 0.001 # accuracy
isIn = [ bbPath.contains_point(pnt,radius=r) or bbPath.contains_point(pnt,radius=-r) for pnt in pnts]

The result is

[True, True, True, True]

By default (or r=0) all the points on borders are not included, and the result is

[False, False, False, False]
Sergey
  • 487
  • 3
  • 7
2

Here is the correct code that includes edges:

def point_inside_polygon(x, y, poly, include_edges=True):
    '''
    Test if point (x,y) is inside polygon poly.

    poly is N-vertices polygon defined as 
    [(x1,y1),...,(xN,yN)] or [(x1,y1),...,(xN,yN),(x1,y1)]
    (function works fine in both cases)

    Geometrical idea: point is inside polygon if horisontal beam
    to the right from point crosses polygon even number of times. 
    Works fine for non-convex polygons.
    '''
    n = len(poly)
    inside = False

    p1x, p1y = poly[0]
    for i in range(1, n + 1):
        p2x, p2y = poly[i % n]
        if p1y == p2y:
            if y == p1y:
                if min(p1x, p2x) <= x <= max(p1x, p2x):
                    # point is on horisontal edge
                    inside = include_edges
                    break
                elif x < min(p1x, p2x):  # point is to the left from current edge
                    inside = not inside
        else:  # p1y!= p2y
            if min(p1y, p2y) <= y <= max(p1y, p2y):
                xinters = (y - p1y) * (p2x - p1x) / float(p2y - p1y) + p1x

                if x == xinters:  # point is right on the edge
                    inside = include_edges
                    break

                if x < xinters:  # point is to the left from current edge
                    inside = not inside

        p1x, p1y = p2x, p2y

    return inside

Updated: fixed a bug