Pages

Tuesday, August 23, 2011

Point in Polygon 2: Walking the line

Credit: SpikedMath.com
This post is a follow-up to my original article on testing if a point is inside a polygon.  Reader Sebastian V. pointed out the ray-casting alogrithm I used does not test to see if the point is on the edge of the polygon or one of the verticies.  He was even nice enough to send a PHP script he found which uses an indentical ray-casting method and includes a vertex and edge test as well.

These two checks are relatively simple however whether they are necessary is up to you and how you apply this test.  There are some cases where a boundary point would not be considered for inclusion.  Either way now you have an option.  This function could even be modified to optionally check for boundary points.

# Improved point in polygon test which includes edge
# and vertex points

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"
      
   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:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

# Test a vertex for inclusion
poligono = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016

print point_in_poly(lat, lon, poligono)

# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 3
y = 1
print point_in_poly(x, y, poly2)
You can download this script here.

18 comments:

  1. Hi, i'm new in python so now i have this problem.

    Whe i put this
    poly = [(-33.400138,-70.589486),(-33.399673,-70.587834),(-33.402180,-70.587255),(-33.401858,-70.589293),(-33.400138,-70.589486)]
    it works fine, but when i put this
    poly="-33.385909,-70.591407|-33.387306,-70.593467|-33.389349,-70.590076|-33.387342,-70.588703|-33.385909,-70.591407"
    and this to transform
    poligono = []
    pol_split = poly.split("|")
    contar = len(pol_split)
    for i in range(contar):
    coord_split = pol_split[i].split(",")
    lat_geo = float(coord_split[0])
    lon_geo = float(coord_split[1])
    poligono.append( (lat_geo, lon_geo) )

    print poligono

    Is not the same and the result it's bad.
    So, can you help me please, to transform this
    poly="-33.385909,-70.591407|-33.387306,-70.593467|-33.389349,-70.590076|-33.387342,-70.588703|-33.385909,-70.591407"
    to this
    poly = [(-33.400138,-70.589486),(-33.399673,-70.587834),(-33.402180,-70.587255),(-33.401858,-70.589293),(-33.400138,-70.589486)]

    thanks

    ReplyDelete
  2. sorry, i want to transform the format, not the numbers

    ReplyDelete
  3. Do this:

    poly="-33.385909,-70.591407|-33.387306,-70.593467|-33.389349,-70.590076|-33.387342,-70.588703|-33.385909,-70.591407"

    poligono = [tuple(c.split(",")) for c in p.split("|")]

    ReplyDelete
  4. Whoops - small typo. Do this instead:

    poly="-33.385909,-70.591407|-33.387306,-70.593467|-33.389349,-70.590076|-33.387342,-70.588703|-33.385909,-70.591407"

    poligono = [tuple(c.split(",")) for c in poly.split("|")]

    The construct [x for x in aList] is called a python list comprehension.

    ReplyDelete
  5. thanks for the answer, but i have the same problem

    poly_sv = "-33.399274,-70.589905|-33.398629,-70.587759|-33.402068,-70.586987|-33.401925,-70.589046|-33.399274,-70.589905"

    With the formula i got

    [('-33.399274', '-70.589905'), ('-33.398629', '-70.587759'), ('-33.402068', '-70.586987'), ('-33.401925', '-70.589046'), ('-33.399274', '-70.589905')]

    And i need
    [(-33.399274, -70.589905)...
    without ' '

    it's only, because the function doesn`t read well the polygon

    i don`t know how to resolve ir

    ReplyDelete
  6. pairs = poly.split("|")

    coords = [c.split(",") for c in pairs]

    poligono = [(float(i[0]), float(i[1])) for i in coords]

    ReplyDelete
  7. thanks, you really save my day.
    it works perfect.

    thanks

    ReplyDelete
  8. I love the internet. Thanks!

    ReplyDelete
  9. I am having a devil of a time getting this point in poly to work with new Kansas Legislative District Shapefiles (can download at http://redistricting.ks.gov/_Plans/district_court.html) I'm experimenting with Kansas Senate shapefile.

    I can load file, and len(sr) = 40 (as it should - 40 senate districts)

    but when I try to run this point_in_poly function I get this error:

    File "/Users/Jake/scripts/rrgh.py", line 36, in point_in_poly
    if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
    AttributeError: _ShapeRecord instance has no attribute '__getitem__'

    Any insight?

    ReplyDelete
    Replies
    1. Jake,

      I put a very basic working example using your shapefile here that tests a point in District 33:
      http://pastebin.com/YsJ7ZGz8

      You'll probably want to do something more sophisticated but I hope this sample is simple enough to get you started. Let me know if you're still having troubles and thanks for reading the blog.

      Delete
  10. Hi, love this blog!

    I am currently attempting to use this script, but I am getting hung up where I think I need to do some looping. I looked at the working example you pasted above and have downloaded shapefile.py

    I am trying to figure out how to loop through the individual polygons I have in a shapefile. In your working example, you choose just one polygon by setting it to r.shape(10).points. I need to go through each polygon but am not sure how to set that up.

    Essentially I will run through each polygon and determine which points I have sit on vertices/boundaries, if any, but I am currently not sure how to loop though the "r.shape().points" so that each time, a new polygon is checked.

    Any help would be great! I'm pretty novice in Python, but love the shapefile.py and this script!

    ReplyDelete
  11. Alison,

    You use the "shapes()" method which returns a list of shapes. Let's say you have a list of points called "myPoints" and a function that checks if a point is on an edge called "edge()" that returns True or False...

    edgePoints = []
    for p in myPoints:
    for polygon in r.shapes():
    if edge(p, polygon):
    edgePoints.append(p)

    ReplyDelete
  12. Here's a formatted version:

    http://pastebin.com/tZE57xfT

    ReplyDelete
    Replies
    1. Thanks for this code and the help! I will see what I can come up with, though it is all a little daunting to me.

      Delete
  13. Hi Joel,
    Love the blog :)
    Just running some tests on this and noted that the boundary test fails for cases such as:

    poligono = [(-1,-1), (1,-1),(1,1), (-1,1)]
    lat= -1
    lon= -0.5

    Which loops through like:

    x,y = (-1, -0.5)
    p1 = (-1, -1),
    p2 = (1, -1),
    min(p1[0], p2[0]) = -1,
    max(p1[0], p2[0]) = 1
    p1[1] == p2[1]: True
    p1[1] == y: False
    x > min(p1[0], p2[0]): False
    x < max(p1[0], p2[0]): True


    x,y = (-1, -0.5)
    p1 = (-1, -1),
    p2 = (1, -1),
    min(p1[0], p2[0]) = -1,
    max(p1[0], p2[0]) = 1
    p1[1] == p2[1]: True
    p1[1] == y: False
    x > min(p1[0], p2[0]): False
    x < max(p1[0], p2[0]): True


    x,y = (-1, -0.5)
    p1 = (1, -1),
    p2 = (1, 1),
    min(p1[0], p2[0]) = 1,
    max(p1[0], p2[0]) = 1
    p1[1] == p2[1]: False
    p1[1] == y: False
    x > min(p1[0], p2[0]): False
    x < max(p1[0], p2[0]): True


    x,y = (-1, -0.5)
    p1 = (1, 1),
    p2 = (-1, 1),
    min(p1[0], p2[0]) = -1,
    max(p1[0], p2[0]) = 1
    p1[1] == p2[1]: True
    p1[1] == y: False
    x > min(p1[0], p2[0]): False
    x < max(p1[0], p2[0]): True


    OUT

    ReplyDelete
  14. No offense, but no one should use this code. The best solution appears to be implemented in Matplotlib : http://matplotlib.org/faq/howto_faq.html#test-whether-a-point-is-inside-a-polygon

    ReplyDelete
  15. Hi I was wonder the test on the edges works? Because I tested with my codes and it doesn't work.

    ReplyDelete