1

So I have written some code to check whether a point [(x,y)] lies within a polygon [(x,y), (x,y), (x,y), (x,y)] but the code fails if the point lies on a boundary or a vertex, I need this to be classed as lying within the polygon. This is the code:

def areasign(poly):
    sumarea = 0
    for i in range(0, len(poly)-1):
        xi = poly[i][0]
        yi = poly[i][1]
        xj = poly[i+1][0]
        yj = poly[i+1][1]
        sumarea = sumarea + ((xi*yj)-(yi*xj))

    if sumarea == 0: return 0
    if sumarea < 0:  return -1
    if sumarea > 0:  return 1

def lineintersection(xy, ij):
    xyi = []
    xyj = []
    ijx = []
    ijy = []
    xyi.extend(xy)
    xyj.extend(xy)
    xyi.append(ij[0])
    xyj.append(ij[1])
    xyi.append(xy[0])
    xyj.append(xy[0])
    ijx.extend(ij)
    ijy.extend(ij)
    ijx.append(xy[0])
    ijy.append(xy[1])
    ijx.append(ij[0])
    ijy.append(ij[0])
    a = areasign(xyi)
    b = areasign(xyj)
    c = areasign(ijx)
    d = areasign(ijy)
    if (a!= b) and (c!=d):
        return True
    else:
        return False


def openpointfile(fname):
    p = []
    f = open(fname)
    for line in f:
        line = line.replace('POINT(', '')
        line = line.replace(')', '')
        vals = line.split(' ')
        res = map(float, vals)
        ctuple = tuple(res)
        p.append(ctuple)
        print ctuple
    f.close()
    return p

points =         openpointfile('C:/Users/Dan/Documents/Informatics/Data/Point_On_Line.txt')

def openpolygonfile(fname):
    p = []
    f = open(fname)
    for line in f:
        line = line.replace('POLYGON((', '')
        line = line.replace('))', '')
        s = line.split(',')
        poly = []
        for coord in s:
            vals = coord.split(' ')
            res = map(float, vals)
            ctuple = tuple(res)
            poly.append(ctuple)
        p.append(poly)      
    f.close()
    return p

polys = openpolygonfile('C:/Users/Dan/Documents/Informatics/Data/Polygon_On_Line.txt')

def pointinpoly(points, polys):
    infinity = tuple([100000000.0, 100000010.0])
    sum = 0
    testline = (points, infinity)
    for i in range(0, len(poly)-1):
        start = poly[i]
        end = poly[i+1]
        line = (start, end)
        l = lineintersection(line, testline)
        if l == True:
            sum = sum + 1
    if sum % 2 == 0:
        return  False
    else:
        return True

pointcount = 0
for point in points:
    pointcount = pointcount + 1
    polycount = 0
    for poly in polys:
        polycount = polycount + 1
        l = pointinpoly(point, poly)
        res = 'Point ' + str(pointcount) + ' lies within Polygon ' + str(polycount) + ' = ' + str(l)
        print l
        print (res)

Like i said before this doesnt work if the point lies on a boundary or a vertex, i have this function that i found on here that seems to work but i can't find a way to get it my files into the correct format for the function and to get it to iterate oer a list of polygons.

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

Any help on how to build that function into my existing code or any help on a new function that could be used to check if the point lies on the boundary of a polygon or on a vertex

3

whether a point [(x,y)] lies within a polygon [(x,y), (x,y), (x,y), (x,y)] but the code fails if the point lies on a boundary or a vertex

You need to understand the topologic rules:

Polygons do not contain their boundary (Geometry A contains Geometry B if no points of B lie in the exterior of A, and at least one point of the interior of B lies in the interior of A)

Contains Returns True if the object’s interior contains the boundary and interior of the other object and their boundaries do not touch at all.

Thus, it is normal that your code fails if the point lies on a boundary or a vertex

A polygon, its LinearRing (boundary) and a vertex point

enter image description here

1) Control of the predicates with the Java Topology Suite (same for polygon and a vertex point)

enter image description here

2) Control with Shapely

from shapely.geometry import Point, Polygon
polygon = Polygon([(0, 0), (1, 1), (1, 0)])
# boundary of the polygon = LinearRing
linearring = LinearRing(list(polygon.exterior.coords))
print linearring
LINEARRING (0 0, 1 1, 1 0, 0 0)
# contains
polygon.contains(linearring)
False
# a vertex 
point = Point(1,1)
polygon.contains(point)
False

But

polygon.touches(linearring)
True
polygon.touches(point)
True
polygon.intersect(linearring)
True
polygon.intersect(point)
True

And

linearring.contains(point)
True
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.