Find if point is located inside or on the border of a polygon (Python - Shapely)

Question

Well, I need to check if a point is within or located on a border of a polygon, but for this particular case I can´t get it working. First some functions and variables:

from matplotlib import pyplot as plt    
from shapely.geometry import Polygon,Point,LinearRing,LineString
from math import sqrt
from shapely import affinity
from descartes.patch import PolygonPatch

GM = (sqrt(5)-1.0)/2.0
W = 8.0
H = W*GM
SIZE = (W, H)

BLUE = '#6699cc'
GRAY = '#999999'
DARKGRAY = '#333333'
YELLOW = '#ffcc33'
GREEN = '#339933'
RED = '#ff3333'
BLACK = '#000000'

COLOR_ISVALID = {
    True: BLUE,
    False: RED,
}

def plot_line(ax, ob, color=GRAY, zorder=1, linewidth=3, alpha=1):
    x, y = ob.xy
    ax.plot(x, y, color=color, linewidth=linewidth, solid_capstyle='round', zorder=zorder, alpha=alpha)

def plot_coords(ax, ob, color=GRAY, zorder=1, alpha=1):
    x, y = ob.xy
    ax.plot(x, y, 'o', color=color, zorder=zorder, alpha=alpha)

def color_isvalid(ob, valid=BLUE, invalid=RED):
    if ob.is_valid:
        return valid
    else:
        return invalid

def color_issimple(ob, simple=BLUE, complex=YELLOW):
    if ob.is_simple:
        return simple
    else:
        return complex

def plot_line_isvalid(ax, ob, **kwargs):
    kwargs["color"] = color_isvalid(ob)
    plot_line(ax, ob, **kwargs)

def plot_line_issimple(ax, ob, **kwargs):
    kwargs["color"] = color_issimple(ob)
    plot_line(ax, ob, **kwargs)

def plot_bounds(ax, ob, zorder=1, alpha=1):
    x, y = zip(*list((p.x, p.y) for p in ob.boundary))
    ax.plot(x, y, 'o', color=BLACK, zorder=zorder, alpha=alpha)

def add_origin(ax, geom, origin):
    x, y = xy = affinity.interpret_origin(geom, origin, 2)
    ax.plot(x, y, 'o', color=GRAY, zorder=1)
    ax.annotate(str(xy), xy=xy, ha='center',
                textcoords='offset points', xytext=(0, 8))

def set_limits(ax, x0, xN, y0, yN):
    ax.set_xlim(x0, xN)
    ax.set_xticks(range(x0, xN+1))
    ax.set_ylim(y0, yN)
    ax.set_yticks(range(y0, yN+1))
    ax.set_aspect("equal")

Here is my code:

pts_tri = [(23.5, 40.703193977868615), (6.5, 11.258330249197702), (19.650283067421512,12.213053056663625)]
test = (8.258506615, 14.304153051823665)
poly = Polygon(pts_tri)
linearring = LinearRing(list(poly.exterior.coords))
fig = plt.figure(figsize=(10,5))
ax2 = plt.subplot(1,2,1)
for point in pts_tri:
    plt.plot(point[0],point[1],'o',color='black')
plot_coords(ax2, poly.exterior)
patch = PolygonPatch(poly, facecolor=color_isvalid(poly), edgecolor=color_isvalid(poly, valid=BLUE), alpha=0.5, zorder=2)
ax2.add_patch(patch)
plt.plot(test[0],test[1],'o', color='red')
plt.show()
print((poly.contains(Point(test))),(poly.touches(Point(test))),linearring.intersects(Point(test)),poly.overlaps(Point(test)))

Result of the plt.show()

and the print result is: False False False False

The point is located exactly at the polygon border, but all I can get is False. What should I do? Thanks in advance.

Please edit this to include minimum viable code to reproduce the error — Hayden Eastwood, Feb 14 '20 at 18:38
It was missing the plt import. pts_tri is a list of the polygon vertices and test is the test point. — Altanir Junior, Feb 14 '20 at 18:58
I've had a look at this and can't solve it - will return to it in the morning. Have you tried putting points more obviously in the polygon area? — Hayden Eastwood, Feb 14 '20 at 21:57
Thanks. I´m using this code to verify points on a rock classification diagram, using thousands of points. When the point is not on a border, the shapely test algorithm works as expected. — Altanir Junior, Feb 17 '20 at 10:33

score 1 · Answer 1 · answered Feb 26 '20 at 13:55

Shapely intersects can act unexpectedly due to floating point precision issues. If I round your test point...

test = (round(8.258506615, 14), round(14.304153051823665, 14))

Then I can get

poly.contains(Point(test)) to True and

poly.intersection(Point(test)) to POINT (8.258506615 14.30415305182366)

An even better method may be to just check the distance

epsilon = 1e-15
print(Point(test).distance(poly) < epsilon)
>>> True

A more technical explanation can be found here

Find if point is located inside or on the border of a polygon (Python - Shapely)

1 Answers1