Point in Polygon in Python

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.

Similar Posts:

    None Found

4 Comments

  1. Dan:

    I was looking for the exact same thing two days ago. I found this snippet on the Python list from somebody calling themselves Guido ;-)

    http://mail.python.org/pipermail/python-list/1999-September/012609.html

    Do a search for Paul Burke’s original page (the link in the mail is dead), as it explains the algorithms a lot more clearly.

  2. Justin Bronn:

    You could’ve just stopped at GeoDjango :) . It includes the original ctypes interface to GEOS, and has Point and Polygon objects that expose all of the geometric operations that GEOS provides. Good suggestion about having “point in polygon” phrase in the docs, and I’ll try to incorporate it into GeoDjango’s GEOS docs: http://geodjango.org/docs/geos.html. Here’s an example:

    >>> from django.contrib.gis.geos import Point, Polygon
    >>> pnt = Point(1, 1)
    >>> poly = Polygon( ((0, 0), (0, 2), (2, 2), (2, 0), (0, 0)) )
    >>> poly.contains(pnt)
    True
    >>> pnt.within(poly)
    True

  3. Neil:

    For other google searchers, there is a fast version of this code in the python plotting library, matplotlib:

    >>> from matplotlib.nxutils import points_inside_poly
    >>> help(points_inside_poly)
    points_inside_poly(…)
    mask = points_inside_poly(xypoints, xyverts)

    Return a boolean ndarray, True for points inside the polygon.

    *xypoints*
    a sequence of N x,y pairs.
    *xyverts*
    sequence of x,y vertices of the polygon.
    *mask* an ndarray of length N.

    A point on the boundary may be treated as inside or outside.
    See `pnpoly `_