Polygons Intersection Points

[AhsanAbbas]

I want to calculate the intersection point of two polygon, the polygons are illustrate in picture. The picture clearly show that there are two points

of intersections.

I have write some code, but there is problem init, it's calculate more than two intersection point. Blow is my code...

def linesIntersection(A, B, C, D):

    x1 = A[0]

    y1 = A[1]

    x2 = B[0]

    y2 = B[1]

    x3 = C[0]

    y3 = C[1]

    x4 = D[0]

    y4 = D[1]

    #print x1, y1, x2, y2

    #print x3, y3, x4, y4

    denomenator = (((x1-x2)*(y3-y4)) - ((y1-y2)*(x3-x4)))

    #print denomenator

    if denomenator == 0:

        return denomenator

    else:

        P1 = (((x1*y2 - y1*x2)*(x3-x4)) - ((x1-x2)*(x3*y4-y3*x4)))

        Px = P1 / denomenator

        P2 = (((x1*y2 - y1*x2)*(y3-y4)) - ((y1-y2)*(x3*y4-y3*x4)))

        Py = P2 / denomenator

    print "The Intersection points are: ", Px, Py

polygon1 = [(1,0),(3,0),(3,2),(1,2),(1,0)]

polygon2 = [(0,1),(2,1),(2,3),(0,3),(0,1)]

poly1 = len(polygon1)

poly2 = len(polygon1)

#print poly1,poly2

for i in range(poly1-1):

    A = polygon1

    B = polygon1[i+1]

    for j in range(poly2-1):

        C = polygon2

        D = polygon2[j+1]

        linesIntersection(A, B, C, D)

I am using determinant for calculating the intersection of point. http://en.wikipedia.org/wiki/Line%E2%80%93line_intersection

There is a condition both lines should be similar planes, this is where problem begin, i don't know how to do it ? Pleas help me...

[DanPatterson_Retired]

There are more than one line intersections...

Why don't you try this demo script I use for education purposes...it uses "ray-casting"  It will give you some ideas, so it is quite verbose.

'''
PointInPolygonRayCasting.py
  Determine if a point is inside a given polygon or not
  Polygon is a list of [x,y] pairs, either as a list or tuple.
  A quick check is made to see if the pnt is one of the polygon's forming
  pnts.  Points that form a polygon are considered inside.  
  More complex analysis would check to see whether pnt, intersects the lines formed
  by the polygon points.
  
  This fuction returns True or False.  The algorithm is called
  "Ray Casting Method".
'''
def point_in_poly(pnt, poly, epsilon):
  '''Point in Polygon test using the Ray Casting Method.
     pnt = [x, y]
     poly = a list of pnts in clockwise order and the first and last are the same
     epsilon = the absolute value of the deviation to be considered if a point is
               within a polygon.  A small value such as 1e9 should suffice considering
               that is being used with real world coordinates.
     Sources: various but http://paulbourke.net/geometry/insidepoly/  or
              http://geospatialpython.com/2011_01_01_archive.html are good
  '''
  n = len(poly)
  inside = False
  x, y = pnt
  p1_X, p1_Y = poly[0]
  p2_X, p2_Y = poly[1]
  if pnt in poly:                 #Step 1 quick check to see if a pnt is one of the polygon pnts.
    p1_X, p1_Y = p2_X, p2_Y
    return True
  else:                           #Step 2  check for obvious cases as to whether points inside or outside
    for i in range(n + 1):            
      p2_X, p2_Y = poly[i % n]
      if y > min(p1_Y, p2_Y):
        if y <= max(p1_Y, p2_Y):
          if x <= max(p1_X, p2_X):
            if p1_Y != p2_Y:
              x_intercept = (y - p1_Y)*(p2_X - p1_X)/(p2_Y - p1_Y) + p1_X
              print "x_intercept: ", x_intercept
              inside = True
              return inside
            if p1_X == p2_X or x <= x_intercept:
              inside = not inside
      p1_X, p1_Y = p2_X, p2_Y
    return inside
if __name__ == "__main__":
  ''' Test
      Note point on boundary fails ie (0,10), floats don't matter
  '''
  polygon = [[0.0,0.0],[0.0,10.0],[10.0,10.0],[10.0,0.0],[0.0,0.0]]
  print "\nPoint in Polygon testing:  Note, points on boundary are not considered inside"
  print "\nPolygon boundary: ", polygon
  pnts = [[0,0],[0.001,0], [0.001,0.001],[1,5], [-10,-10],[9.9999999,9.9999999]]
  for pnt in pnts:## Call the fuction with the points and the polygon
    print "point: ", pnt, " in polygon? ", point_in_poly(pnt, polygon, 1e-9)

[XanderBakker]

Just a silly thought... since we are at the GeoNet (ArcGIS) forum and talking about python code, why don't we use some arcpy "magic"...

import arcpy
lst_pol1 = [(1,0),(3,0),(3,2),(1,2),(1,0)]
lst_pol2 = [(0,1),(2,1),(2,3),(0,3),(0,1)]

polygon1 = arcpy.Polygon(arcpy.Array([arcpy.Point(a[0],a[1]) for a in lst_pol1]))
polygon2 = arcpy.Polygon(arcpy.Array([arcpy.Point(a[0],a[1]) for a in lst_pol2]))

dimension = 1
pnts = polygon1.intersect(polygon2, dimension)
for pnt in pnts:
    print pnt

returns:

2,0001220703125 2,0001220703125 NaN NaN

1,0001220703125 1,0001220703125 NaN NaN

[DanPatterson_Retired]

his previous posts Xander suggest that he wanted pure python implementations and not the useful "magic" shortcuts   left out numpy implementation for the moment

ps ... make sure you specify a spatial reference for the points while being accessed otherwise the results are not cast in the appropriate precision (I posted a number of threads on this issue and it has been seen by others).  It can affect some results like intersections etc when the geometries are close to within numeric precision.)

[AhsanAbbas]

Actually i don't understand how your point in polygon code helps me, you need to tell, what is problem in my code. ?

[XanderBakker · Accepted solution]

Hi,to answer our original question. The extra intersection points are caused by the code which checks for intersection of endless lines. The lines extent to both sides into infinity. To get the right intersections you should check if the resulting points are in the "domains" of the input lines. I'm sure there is an easier way to do this, but see below how I did this.

def linesIntersection(A, B, C, D):

    x1 = A[0]

    y1 = A[1]

    x2 = B[0]

    y2 = B[1]

    x3 = C[0]

    y3 = C[1]

    x4 = D[0]

    y4 = D[1]

    # print x1, y1, x2, y2, "     ", x3, y3, x4, y4

    # print x3, y3, x4, y4

    denomenator = (((x1-x2)*(y3-y4)) - ((y1-y2)*(x3-x4)))

    #print denomenator

    if denomenator == 0:

        return denomenator

    else:

        P1 = (((x1*y2 - y1*x2)*(x3-x4)) - ((x1-x2)*(x3*y4-y3*x4)))

        Px = P1 / denomenator

        P2 = (((x1*y2 - y1*x2)*(y3-y4)) - ((y1-y2)*(x3*y4-y3*x4)))

        Py = P2 / denomenator

##    print "Px: {0},  min(x)={1}, max(x)={2}".format(Px, min(x1, x2), max(x1, x2))

##    print "Px: {0},  min(x)={1}, max(x)={2}".format(Px, min(x3, x4), max(x3, x4))

##    print "Py: {0},  min(y)={1}, max(y)={2}".format(Py, min(y1, y2), max(y1, y2))

##    print "Py: {0},  min(y)={1}, max(y)={2}".format(Py, min(y3, y4), max(y3, y4))

    if Px >= min(x1, x2) and Px <= max(x1, x2):

        if Px >= min(x3, x4) and Px <= max(x3, x4):

            if Py >= min(y1, y2) and Py <= max(y1, y2):

                if Py >= min(y3, y4) and Py <= max(y3, y4):

                    print " - The Intersection points are: ", Px, Py

polygon1 = [(1,0),(3,0),(3,2),(1,2),(1,0)]

polygon2 = [(0,1),(2,1),(2,3),(0,3),(0,1)]

poly1 = len(polygon1)

poly2 = len(polygon1)

# print poly1,poly2

for i in range(0, poly1-1):

    A = polygon1

    B = polygon1[i+1]

    for j in range(0, poly2-1):

        C = polygon2

        D = polygon2[j+1]

        linesIntersection(A, B, C, D)

Kind regards, Xander

edit: for some strange reason the code isn't formatting...

[AhsanAbbas]

Yes it's works, It's true only for the case of both polygons have same no of points, other wise it's index is out of bound problem...

Well thanks for help...