Python point location of intersect between two lines

[Madi]

import arcpy
import math
import numpy

class Toolbox(object):
def __init__(self):
"""Define the toolbox (the name of the toolbox is the name of the
.pyt file)."""
self.label = "LogLocator"
self.alias = "LogLocator"

# List of tool classes associated with this toolbox
self.tools = [Locator1]

class Locator1(object):
def __init__(self):
"""Define the tool (tool name is the name of the class)."""
self.label = "Type 1 Locator"
self.description = "Parse the LogLocator input strings and create the location"
self.canRunInBackground = False

def getParameterInfo(self):
"""Define parameter definitions"""
param0 = arcpy.Parameter (
displayName="PLSS Feature Class",
name="PLSS",
datatype="GPFeatureLayer",
parameterType="Required",
direction="Input")

param1 = arcpy.Parameter (
displayName="County Feature Class",
name="CountyFC",
datatype="GPFeatureLayer",
parameterType="Required",
direction="Input")

param0.filter.list = ["Polygon"]
param1.filter.list = ["Polygon"]

param2 = arcpy.Parameter(
displayName="Sec-Township-Range",
name="STR",
datatype="String",
parameterType="Required",
direction="Input")
param3 = arcpy.Parameter(
displayName="Footage from Lines",
name="ftFL",
datatype="String",
parameterType="Optional",
direction="Input")
param4 = arcpy.Parameter(
displayName="County",
name="County",
datatype="String",
parameterType="Required",
direction="Input")
param5 = arcpy.Parameter(
displayName="State",
name="State",
datatype="String",
parameterType="Required",
direction="Input")

# #################################
#Test Data
param0.value = "PLSS_Sections"
param1.value = "County"
##
## param2.value = "54-T7N-R2W"
## param3.value = "100 FNL 1000 FEL"
## param4.value = "Adams"
## param5.value = "MS"
##
## param2.value = "3-T16N-R10E"
## param3.value = "1000 FNL 1000 FEL"
## param4.value = "Choctaw"
## param5.value = "MS"

param2.value = "13-T8S-R4W"
param3.value = "1000 FNL 1000 FEL"
param4.value = "Lafayette"
param5.value = "MS"

params = [param0, param1, param2, param3, param4, param5]
return params

def isLicensed(self):
"""Set whether tool is licensed to execute."""
return True

def updateParameters(self, parameters):
"""Modify the values and properties of parameters before internal
validation is performed. This method is called whenever a parameter
has been changed."""

# Enable text input dialog only if features selected
if parameters[0].altered and parameters[1].altered == True:
parameters[2].enabled = True
parameters[3].enabled = True
parameters[4].enabled = True
parameters[5].enabled = True
else:
parameters[2].enabled = False
parameters[3].enabled = False
parameters[4].enabled = False
parameters[5].enabled = False
return

def updateMessages(self, parameters):
"""Modify the messages created by internal validation for each tool
parameter. This method is called after internal validation."""
return

def execute(self, parameters, messages):

# Function to test for line intersections
def islineintersect(line1, line2):
## arcpy.AddMessage("{} {}".format(line1,line2))
i1 = [min(line1[0][0], line1[1][0]), max(line1[0][0], line1[1][0])]
i2 = [min(line2[0][0], line2[1][0]), max(line2[0][0], line2[1][0])]
ia = [max(i1[0], i2[0]), min(i1[1], i2[1])]
if max(line1[0][0], line1[1][0]) < min(line2[0][0], line2[1][0]):
return False
m1 = (line1[1][1] - line1[0][1]) * 1. / (line1[1][0] - line1[0][0]) * 1.
m2 = (line2[1][1] - line2[0][1]) * 1. / (line2[1][0] - line2[0][0]) * 1.
if m1 == m2:
return False
b1 = line1[0][1] - m1 * line1[0][0]
b2 = line2[0][1] - m2 * line2[0][0]
x1 = (b2 - b1) / (m1 - m2)
if (x1 < max(i1[0], i2[0])) or (x1 > min(i1[1], i2[1])):
return False
return True

# Geometry functions
def vectorize(inputline):
x1 = inputline[0][0]
x2 = inputline[1][0]
y1 = inputline[0][1]
y2 = inputline[1][1]
return [(x2-x1),(y2-y1)]
def length(v):
return math.sqrt(v[0]**2+v[1]**2)
def dot_product(v,w):
return v[0]*w[0]+v[1]*w[1]
def determinant(v,w):
return v[0]*w[1]-v[1]*w[0]
def inner_angle(v,w):
cosx=dot_product(v,w)/(length(v)*length(w))
rad=math.acos(cosx) # in radians
return rad # return rad*180/math.pi for angle in degrees
def angle_clockwise(A, B):
inner=inner_angle(A,B)
det = determinant(A,B)
if det<0: #this is a property of the det. If the det < 0 then B is clockwise of A
return inner
else: # if the det > 0 then A is immediately clockwise of B
return (2*math.pi - inner) # If in degrees 360-inner

def convert_points_to_line(pnt_array, SR):
# NOTE: The feature class name is hardwired here
# This should be a parameter in any final code
outPath = arcpy.Describe(parameters[0].value).path
inPtFC = arcpy.Describe(parameters[0].value).path + "\\SectPTS"
outFC = arcpy.Describe(parameters[0].value).path + "\\SectLINES"
for lyr in arcpy.mapping.ListLayers(mxd, "", dataFrame):
if lyr.name == "SectLINES":
arcpy.mapping.RemoveLayer(dataFrame,lyr)
if arcpy.Exists(outFC):
arcpy.Delete_management(outFC)
arcpy.PointsToLine_management(inPtFC,outFC,"idfield")
addLayer = arcpy.mapping.Layer(outFC)
arcpy.mapping.AddLayer(dataFrame, addLayer)

def output_point_location(pnt_array, SR):
# NOTE: The feature class name is hardwired here
# This should be a parameter in any final code
outFC = arcpy.Describe(parameters[0].value).path + "\\SectPTS"
for lyr in arcpy.mapping.ListLayers(mxd, "", dataFrame):
if lyr.name == "SectPTS":
arcpy.mapping.RemoveLayer(dataFrame,lyr)
if arcpy.Exists(outFC):
arcpy.Delete_management(outFC)
arcpy.da.NumPyArrayToFeatureClass(pnt_array, outFC, ["ID","XY"], SR)
addLayer = arcpy.mapping.Layer(outFC)
arcpy.mapping.AddLayer(dataFrame, addLayer)
#arcpy.SelectLayerByLocation_management("Location", "WITHIN", parameters[0].value, "", "NEW_SELECTION", "NOT_INVERT")
return

# ##############################
# Parse Input Data
# ##############################

# Get the Section-Township-Range
STR = parameters[2].valueAsText
STR = STR.strip()
name = STR.split('-')
if len(name) != 3:
arcpy.AddWarning("STR string ill-formed")
return
else:
Sec = name.pop(0)
twnshp = name.pop(0)
rnge = name.pop(0)

Township = twnshp[1:-1]
tdir = twnshp[-1]
if tdir != 'N' and tdir != 'S':
arcpy.AddWarning("Township direction not N or S")
Range = rnge[1:-1]
rdir = rnge[-1]
if rdir != 'E' and rdir != 'W':
arcpy.AddWarning("Range direction not E or W")

arcpy.AddMessage("\nSection {0}, Township {1}{2}, Range {3}{4}\n".format(Sec,Township,tdir,Range,rdir))

# Get the ft FNL ft FEL
ftFromLine = parameters[3].valueAsText
# This is optional, so test if there is any input
if ftFromLine != None:
#ftFL = ftFromLine.rstrip()
name = ftFromLine.split() # split without an argument splits on whitespace
ftFU = int(name.pop(0))
FU = name.pop(0)
FU = FU[1:-1]
ftFS = int(name.pop(0))
FS = name.pop(0)
FS = FS[1:-1]
# arcpy.AddMessage("\n {0}ft F{1}L {2}ft F{3}L\n".format(ftFU,FU,ftFS,FS))
else:
ftFU = 0
ftFS = 0

# Get the County
County = parameters[4].valueAsText

# Get the State
State = parameters[5].valueAsText

# arcpy.AddMessage("\n {0} County, {1}\n".format(County,State))
#arcpy.AddMessage(outparams)

# ####################################################
# Find the correct S-T-R
# ####################################################

arcpy.SelectLayerByAttribute_management(parameters[1].value, "NEW_SELECTION", "CONAME = \'"+County+"\'")
arcpy.SelectLayerByLocation_management(parameters[0].value, "INTERSECT", parameters[1].value, "", "NEW_SELECTION", "NOT_INVERT")
arcpy.SelectLayerByAttribute_management(parameters[0].value, "SUBSET_SELECTION", "SECTION = "+Sec+" AND TOWNSHIP = "+Township+" AND TOWNSHIP_D = \'"+tdir+"\' AND RANGE = "+Range+" AND RANGE_D = \'"+rdir+"\'")
arcpy.SelectLayerByAttribute_management(parameters[1].value, "CLEAR_SELECTION", "")

# ####################################################
# Calculate footage displacements
# ####################################################

# Get the mapping parameters
mxd = arcpy.mapping.MapDocument("Current")
dataFrame = arcpy.mapping.ListDataFrames(mxd, "*")[0]
if dataFrame.mapUnits == 'Meters':
# conversion factor from feet to meters
ftFU = ftFU * 0.3048
ftFS = ftFS * 0.3048
SpatialRef = arcpy.Describe(parameters[0].value).spatialReference

centroid = arcpy.da.FeatureClassToNumPyArray(parameters[0].value,["SHAPE@XY"],explode_to_points=False)
centroid = centroid[0][0][0:2]

# Determine if there are any offsets
if ftFU == 0 or ftFS == 0:
# If no offsets center the location in the section
# In this case, we don't need to do anything else
# look at this for executing point location-once you find the correct input
pnt_location = [centroid[0],centroid[1]]
pnt_array = numpy.array([(1,(centroid[0],centroid[1]))],numpy.dtype([('idfield',numpy.int32),('XY', '<f8', 2)]))
output_point_location(pnt_array, SpatialRef)
return

# Otherwise we need to figure out how to interpret the offsets
pnt_array = arcpy.da.FeatureClassToNumPyArray(parameters[0].value,["OID@","SHAPE@XY"],explode_to_points=True)
pnt_array.dtype = [('idfield',numpy.int32),('XY', '<f8', 2)]
# Renumber the point ids
for indx in range(0,len(pnt_array)):
pnt_array[indx][0] = indx

# Determine orientations of all the line segments of the section polygon
LineTypeVect = []
for indx in range(0,len(pnt_array)-1):
mod = len(pnt_array)-1
#arcpy.AddMessage("({0} {1}) ({2} {3})".format((indx % mod),((indx+1) % mod),((indx+1) % mod),((indx+2) % mod)))
#Determine the clockwise angle between the line segment and the horizontal.
# Lines near horizontal are north or south, lines near vertical are east or west
# N-S or E-W determined by relationship to centroid of section
p1 = [pnt_array[indx % mod][1][0],pnt_array[indx % mod][1][1]]
p2 = [pnt_array[(indx+1) % mod][1][0],pnt_array[(indx+1) % mod][1][1]]
p3 = [pnt_array[(indx+1) % mod][1][0],pnt_array[(indx) % mod][1][1]] #Horizontal to p1
v1 = vectorize([p1,p2])
v2 = vectorize([p1,p3])
LineAngle = angle_clockwise(v1,v2) # Angles are returned in radians (not degrees)
if (abs(math.sin(LineAngle)) < math.sin(math.pi/4)):
if p1[1] > centroid[1] or p2[1] > centroid[1]:
LineType = "North"

else:
LineType = "South"
else:
if p1[0] > centroid[0] or p2[0] > centroid[0]:
LineType = "East"
else:
LineType = "West"
LineTypeVect.append(LineType)

arcpy.AddMessage("{:.0f} {:.0f}, {:.0f} {:.0f} - {} {} Line".format(p1[0],p1[1],p2[0],p2[1],LineAngle,LineType))

arcpy.SelectLayerByAttribute_management(parameters[0].value, "CLEAR_SELECTION", "")

# ####################################################
# Find the location point
# ####################################################

# Apply offsets
LocationPoint = [0,0]
PolyPointCnt = len(pnt_array)
line_array = numpy.empty((0,3),numpy.dtype([('idfield',numpy.int32),('P1', '<f8', 2),('P2', '<f8', 2)]))
for indx in range(0,len(pnt_array)-1):
p1 = [pnt_array[indx][1][0],pnt_array[indx][1][1]]
p2 = [pnt_array[indx+1][1][0],pnt_array[indx+1][1][1]]
t1 = (FS == 'E' and LineTypeVect[indx] == 'East')
t2 = (FS == 'W' and LineTypeVect[indx] == 'West')
t3 = (FU == 'N' and LineTypeVect[indx] == 'North')
t4 = (FU == 'S' and LineTypeVect[indx] == 'South')
if (t1 or t2 or t3 or t4):
# Shift the endpoints of the edge the offset distance perpendicular to orientation
p1x = p1[0] - ((ftFS/math.sqrt((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2))*(p1[1]-p2[1]))
p1y = p1[1] + ((ftFU/math.sqrt((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2))*(p1[0]-p2[0]))
p2x = p2[0] - ((ftFS/math.sqrt((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2))*(p1[1]-p2[1]))
p2y = p2[1] + ((ftFU/math.sqrt((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2))*(p1[0]-p2[0]))
lineitem = [indx,[p1x,p1y],[p2x,p2y]]
new_line = numpy.array([tuple(lineitem)],dtype=line_array.dtype)
line_array = numpy.append(line_array, new_line)

# Write the locations to a point feature
lineitem = [999+indx,[p1x,p1y]]
extr_point = numpy.array([tuple(lineitem)],dtype=pnt_array.dtype)
pnt_array = numpy.append(pnt_array, extr_point)
lineitem = [999+indx,[p2x,p2y]]
extr_point = numpy.array([tuple(lineitem)],dtype=pnt_array.dtype)
pnt_array = numpy.append(pnt_array, extr_point)

# Now we have the offset lines - which ones intersect?
for indx in range(0,len(line_array)-1):
for jndx in range(indx+1,len(line_array)):
Test = islineintersect([line_array[indx][1],line_array[indx][2]],[line_array[jndx][1],line_array[jndx][2]])
if Test:
arcpy.AddMessage("Line #{} intersects Line #{}\n".format(line_array[indx][0],line_array[jndx][0]))
else:
arcpy.AddMessage("Line #{} does not intersect Line #{}\n".format(line_array[indx][0],line_array[jndx][0]))

for indx in range(0,PolyPointCnt):
pnt_array = numpy.delete(pnt_array,0,0)
output_point_location(pnt_array, SpatialRef)

convert_points_to_line(pnt_array, SpatialRef)

# Find intersecting line location
# given two lines (pnt_array)[ID],[P1 P2][P3 P4], each point is an x,y
# for indx in range(0,len(line_array)-1):

THIS IS MY CODE SO FAR, i have to find the point location of the two lines that intersect and im not sure how to do that, anything helps!