## 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
p2 = poly
else:
p1 = poly[i-1]
p2 = poly[i]
if p1 == p2 and p1 == y and x > min(p1, p2) and x < max(p1, p2):
return "IN"

n = len(poly)
inside = False

p1x,p1y = poly
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)
```

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)
lon_geo = float(coord_split)
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

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

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("|")]

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.

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

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

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

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

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

thanks

8. I love the internet. Thanks!

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 == p2 and p1 == y and x > min(p1, p2) and x < max(p1, p2):
AttributeError: _ShapeRecord instance has no attribute '__getitem__'

Any insight?

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.

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!

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)

12. Here's a formatted version:

http://pastebin.com/tZE57xfT

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.

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, p2) = -1,
max(p1, p2) = 1
p1 == p2: True
p1 == y: False
x > min(p1, p2): False
x < max(p1, p2): True

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

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

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

OUT

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

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

16. I am struggling to get this function to work reading and iterating over 2 wkt files, 1 containing the points and 1 containing the polygons. Any ideas ?

17. you miss one line in the boundaries detection

if p1 == p2 and p1 == x and y > min(p1, p2) and y < max(p1, p2):
return "IN"

check collinearity for both x and why

18. Your posts focus on important things.
Aussie flings