Credit: SpikedMath.com |

*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.

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

ReplyDeleteWhe 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

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

ReplyDeleteDo this:

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

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"

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

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

thanks for the answer, but i have the same problem

ReplyDeletepoly_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

pairs = poly.split("|")

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

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

thanks, you really save my day.

ReplyDeleteit works perfect.

thanks

I love the internet. Thanks!

ReplyDeleteI 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.

ReplyDeleteI 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?

Jake,

DeleteI 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.

Hi, love this blog!

ReplyDeleteI 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!

Alison,

ReplyDeleteYou 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)

Here's a formatted version:

ReplyDeletehttp://pastebin.com/tZE57xfT

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.

DeleteHi Joel,

ReplyDeleteLove 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

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

ReplyDeleteI agree.

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

ReplyDeleteI 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 ?

ReplyDelete