본문 바로가기
Python/Qgis(PyQgis) & ArcGis(Arcpy)

Qgis&Arcgis) arcpy 벡터의 외적을 이용한 두 벡터의 교차점 구하기

by 유노파이 2021. 8. 10.


# 벡터의 외적
def cross_product(x0,y0,x1,y1, x2,y2,x3,y3):
    val = (x1-x0)*(y3-y2) - (y1-y0)*(x3-x2)
	return val

# 두 벡터의 교차점 
def intersection_point(x0,y0,x1,y1, x2,y2,x3,y3, val):
    
    pre  = x1*y0 - y1*x0
    post = x3*y2 - y3*x2
    
    x = (pre * (x3-x2) - (x1-x0) * post) / val
    y = (pre * (y3-y2) - (y1-y0) * post) / val

    return x,y
                 
if __name__ == '__main__':
    
    for idx, row in enumaerat(arcpy.da.SearchCursor("폴리라인","SHAPE@")):
        if idx == 0: # 첫번째 라인
            x0 = row[0].firstPoint.X ## 폴리라인 시작점
            y0 = row[0].firstPoint.Y
            x1 = row[0].lastPoint.X  ## 폴리라인 끝점
            y1 = row[0].lastPoint.Y
        else:
            x2 = row[0].firstPoint.X
            y2 = row[0].firstPoint.Y
            x3 = row[0].lastPoint.X
            y3 = row[0].lastPoint.Y
    
    val = cross_product(x0,y0,x1,y1, x2,y2,x3,y3)
    if val == 0:
    	arcpy.AddMessage('두 벡터는 평행 합니다.')
    	sys.exit()
    
    x,y = intersection_point(x0,y0,x1,y1, x2,y2,x3,y3, val)
    
    # 포인트에 추가 
    insertTotal = arcpy.da.InsertCursor ('포인트', "SHAPE@")   
    insertTotal.insertRow([(x,y)])
    del insertTotal

    arcpy.RefreshActiveView()

 

궁금 사항은 문의주세요~

댓글