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...
Solved! Go to 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...
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)
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
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.)
Actually i don't understand how your point in polygon code helps me, you need to tell, what is problem in my code. ?
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...
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...