Sunday, 28 July 2013

Point in Polygon by Ray Casting

An old classic.   Note that points that fall on the boundary are not considered inside. Tests for boundary conditions would have to be added.  The script runs in standalone mode.

'''
PointInPolygonRayCasting.py
 Determine if a point is inside a given polygon or not
 Polygon is a list of (x,y) pairs. This fuction
 returns True or False.  The algorithm is called
 "Ray Casting Method".

'''
def point_in_poly(pnt, poly):
  '''Point in Polygon test using the Ray Casting Method.
     pnt = [x, y]
     poly = a list of pnts
     Sources: various but
http://paulbourke.net/geometry/insidepoly/  or
             
http://geospatialpython.com/2011_01_01_archive.html are good'''
  n = len(poly)
  inside = False
  x, y = pnt
  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)/(p2y - p1y) + p1x
          if p1x == p2x or x <= xinters:
            inside = not inside
    p1x, p1y = p2x, p2y
  return inside

if __name__ == "__main__":
  ''' Test
      Note point on boundary fails ie (0,10), floats don't matter
  '''
  polygon = [(0.0,0.0),(0.0,10.0),(10.0,10.0),(10.0,0.0),(0.0,0.0)]
  print "\nPoint in Polygon testing:  Note, points on boundary are not considered inside"
  print "\nPolygon boundary: ", polygon
  #point_x = 0.01 #5
  #point_y = 0.01 #5
  pnts = [[0,0],[0.001,0], [0.001,0.001],[1,5], [-10,-10],[9.9999999,9.9999999]]
  #pnt = [point_x, point_y]
  for pnt in pnts:## Call the fuction with the points and the polygon
    print "point: ", pnt, " in polygon? ", point_in_poly(pnt, polygon)

No comments:

Post a Comment