Python point location of intersect between two lines

Discussion in 'Python' started by Madi, Feb 28, 2018.

  1. Madi

    Madi

    Joined:
    Feb 28, 2018
    Messages:
    1
    Likes Received:
    0
    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!
     
    Madi, Feb 28, 2018
    #1
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.