inpolygon for Python - Examples of matplotlib.path.Path contains_points() method?

Question

I have been searching for a python alternative to MATLAB's inpolygon() and I have come across contains_points as a good option.

However, the docs are a little bare with no indication of what type of data contains_points expects:

contains_points(points, transform=None, radius=0.0)

Returns a bool array which is True if the path contains the corresponding point.

If transform is not None, the path will be transformed before performing the test.

radius allows the path to be made slightly larger or smaller.

I have the polygon stored as an n*2 numpy array (where n is quite large ~ 500). As far as I can see I need to call the Path() method on this data which seems to work OK:

poly_path = Path(poly_points)

At the moment I also have the points I wish to test stored as another n*2 numpy array (catalog_points).

Perhaps my problem lies here? As when I run:

in_poly = poly_path.contains_points(catalog_points)

I get back an ndarray containing False for every value no matter the set of points I use (I have tested this on arrays of points well within the polygon).

Answer 1

Often in these situations, I find the source to be illuminating...

We can see the source for path.contains_point accepts a container that has at least 2 elements. The source for contains_points is a bit harder to figure out since it calls through to a C function Py_points_in_path. It seems that this function accepts a iterable that yields elements that have a length 2:

>>> from matplotlib import path
>>> p = path.Path([(0,0), (0, 1), (1, 1), (1, 0)])  # square with legs length 1 and bottom left corner at the origin
>>> p.contains_points([(.5, .5)])
array([ True], dtype=bool)

Of course, we could use a numpy array of points as well:

>>> points = np.array([.5, .5]).reshape(1, 2)
>>> points
array([[ 0.5,  0.5]])
>>> p.contains_points(points)
array([ True], dtype=bool)

And just to check that we aren't always just getting True:

>>> points = np.array([.5, .5, 1, 1.5]).reshape(2, 2)
>>> points
array([[ 0.5,  0.5],
       [ 1. ,  1.5]])
>>> p.contains_points(points)
array([ True, False], dtype=bool)

Answer 2

Make sure that the vertices are ordered as wanted. Below vertices are ordered in a way that the resulting path is a pair of triangles rather than a rectangle. So, contains_points only returns True for points inside any of the triangles.

>>> p = path.Path(np.array([bfp1, bfp2, bfp4, bfp3]))
>>> p
Path([[ 5.53147871  0.78330843]
 [ 1.78330843  5.46852129]
 [ 0.53147871 -3.21669157]
 [-3.21669157  1.46852129]], None)
>>> IsPointInside = np.array([[1, 2], [1, 9]])
>>> IsPointInside
array([[1, 2],
       [1, 9]])
>>> p.contains_points(IsPointInside)
array([False, False], dtype=bool)
>>> 

The output for the first point would have been True if bfp3 and bfp4 were swapped.

Answer 3

I wrote this function to return a array as in matlab inpolygon function. But this will return only the points that are inside the given polygon. You can't find the points in the edge of the polygon with this function.

import numpy as np
from matplotlib import path

def inpolygon(xq, yq, xv, yv):
    shape = xq.shape
    xq = xq.reshape(-1)
    yq = yq.reshape(-1)
    xv = xv.reshape(-1)
    yv = yv.reshape(-1)
    q = [(xq[i], yq[i]) for i in range(xq.shape[0])]
    p = path.Path([(xv[i], yv[i]) for i in range(xv.shape[0])])
    return p.contains_points(q).reshape(shape)

You can call the function as:

xv = np.array([0.5,0.2,1.0,0,0.8,0.5])
yv = np.array([1.0,0.1,0.7,0.7,0.1,1])
xq = np.array([0.1,0.5,0.9,0.2,0.4,0.5,0.5,0.9,0.6,0.8,0.7,0.2])
yq = np.array([0.4,0.6,0.9,0.7,0.3,0.8,0.2,0.4,0.4,0.6,0.2,0.6])
print(inpolygon(xq, yq, xv, yv))

As in the matlab documentation this function,

returns in indicating if the query points specified by xq and yq are inside or on the edge of the polygon area defined by xv and yv.