I run into a use case where I needed to find if coordinates are inside a polygon, in other words the “point in polygon” test was called for. But I didn’t know to search for point in polygon at first, so this post is more of a description of the progress that lead me to the easiest solution for my use case.
What was surprising to me was that there don’t seem to be that many Python implementations out there just for point in polygon test. The first hit I got was Patrick Jordan’s function which was translated from the original C source. I run some tests with that, but for some reason I was getting errors that did not make sense, and decided to translate the other common algorithm from the C source by counting angles:
# See http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ import math class Point: def __init__(self, h, v): self.h = h self.v = v def InsidePolygon(polygon, p): angle = 0.0 n = len(polygon) for i, (h, v) in enumerate(polygon): p1 = Point(h - p.h, v - p.v) h, v = polygon[(i + 1) % n] p2 = Point(h - p.h, v - p.v) angle += Angle2D(p1.h, p1.v, p2.h, p2.v); if abs(angle) < math.pi: return False return True def Angle2D(x1, y1, x2, y2): theta1 = math.atan2(y1, x1) theta2 = math.atan2(y2, x2) dtheta = theta2 - theta1 while dtheta > math.pi: dtheta -= 2.0 * math.pi while dtheta < -math.pi: dtheta += 2.0 * math.pi return dtheta
It is very unpythonic, but does follow the original C implementation very closely. This algorithm worked for my use case, but then I noticed a typo that I had somehow put in with the first version, so in the end both worked for me.
Both of those algorithms have issues with points on the edge of the polygon. However, since my use case is really about polygons on the surface of the earth with relatively sparse points to test for, getting a point on an edge is pretty unlikely. Of course the curvature of the earth will also cause errors, but my polygons are pretty small and given that I don’t require 100% accuracy it doesn’t really matter.
After I had done all that I finally discovered GeoDjango. And reading through some of that documentation lead me to GIS-Python, which lead me to Shapely. And who would have thought it, Shapely has objects like Points and Polygons and the polygon has the method contains. Argh! Had I thought for a minute about the bigger picture to begin with, I would probably have come up with the right keywords to throw at Google and arrive at Shapely’s page directly and not waste all that time. On the other hand, it wouldn’t hurt for Shapely and these other options to actually contain the string “point in polygon” on some page to help other lost souls like me.
- None Found