[Numpy-discussion] how to tell if a point is inside a polygon

[Bill Baxter]

On Tue, Oct 14, 2008 at 8:46 AM, Pierre GM <pgmdevlist@gmail.com> wrote:
>> 2008/10/13 Mathew Yeates <myeates@jpl.nasa.gov>
>>
>> > Is there a routine in scipy for telling whether  a point is inside a
>> > convex 4 sided polygon?
>
> Mathew,
> You could use OGR (www.gdal.org)
>
> Example
> -------------
> import osgeo.ogr as ogr
>
> vert = [(0,0),(0,1),(1,1),(1,0)]
> listvert = [" %s %s" % (x,y) for (x,y) in vert]
> listvert.append(listvert[0])
> geo = ogr.CreateGeometryFromWkt("POLYGON ((%s))" % ','.join(listvert))
>
> querypoint = (0.5, 0.5)
> qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % querypoint)
>
> assert geo.Contains(qpt)
>
> querypoint = (0.5, 1.5)
> qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % querypoint)
>
> assert not geo.Contains(qpt)


If all you really need is a point in convex polygon test, it's
probably a little too trivial to be worth dragging in a dependency on
anything.  But if you may find yourself needing more geometric tests
then it might be a good idea to get a good geom lib now.

As for this test, all you need to do is check that the point is to the
left of each of the edges, taken counter-clockwise.

Here's some code I wrote a while back that does a more general
even-odd test, but should work for your case too:

def inside_shape(p, verts, edges=None):
    """Test whether the point p is inside the specified shape.
    The shape is specified by 'verts' and 'edges'
    Arguments:
    p - the 2d point
    verts - (N,2) array of points
    edges - (N,2) array of vert indices indicating edges
            If edges is None then assumed to be in order. I.e.
              [[0,1], [1,2], [2,3] ... [N-1,0]]

    Returns:
    - True/False based on result of in/out test.

    Uses the 'ray to infinity' even-odd test.
    Let the ray be the horizontal ray starting at p and going to +inf in x.
    """
    verts = npy.asarray(verts)
    if edges is None:
        N = verts.shape[0]
        edges = npy.column_stack([npy.c_[0:N],npy.c_[1:N,0]])

    inside = False
    x,y=p[0],p[1]
    for e in edges:
        v0,v1 = verts[e[0]],verts[e[1]]
        # Check if both verts to the left of ray
        if v0[0]<x and v1[0]<x:
            continue
        # check if both on the same side of ray
        if (v0[1]<y and v1[1]<y) or (v0[1]>y and v1[1]>y):
            continue
        #check for horizontal line - another horz line can't intersect it
        if (v0[1]==v1[1]):
            continue
        # compute x intersection value
        xisect = v0[0] + (v1[0]-v0[0])*((y-v0[1])/(v1[1]-v0[1]))
        if xisect >= x:
            inside = not inside
    return inside

License: public domain.

--bb