Finding angle between intersecting features in two feature classes using ArcGIS Desktop and Python [closed]

Question

I have two intersecting line featureclasses. I want to find the angle at each point of intersection using ArcGIS 10 and Python.

Can anyone help?

Answer 1

There is a relatively simple workflow. It overcomes the potential problems that two features may intersect in more than one point. It does not require scripting (but can readily be turned into a script). It can be done primarily from the ArcGIS menu.

The idea is to exploit a layer of intersection points, one point for each distinct pair of intersecting polylines. You need to obtain a small piece of each intersecting polyline at these intersection points. Use the orientations of these pieces to compute their intersection angles.

Here are the steps:

1. Make sure each of the polyline features has a unique identifier within its attribute table. This will be used later to join some geometric attributes of the polylines to the intersection point table.

2. Geoprocessing|Intersect obtains the points (make sure to specify you want points for the output).

3. Geoprocessing|Buffer lets you buffer the points by a tiny amount. Make it really tiny so that the portion of each line within a buffer does not bend.

4. Geoprocessing|Clip (applied twice) limits the original polyline layers to just the buffers. Because this produces new datasets for its output, subsequent operations will not modify the original data (which is a good thing).

   

   Here is a schematic of what happens: two polyline layers, shown in light blue and light red, have produced dark intersection points. Around those points tiny buffers are shown in yellow. The darker blue and red segments show the results of clipping the original features to these buffers. The rest of the algorithm works with the dark segments. (You cannot see it here, but a tiny red polyline intersects two of the blue lines at a common point, producing what appears to be a buffer around two blue polylines. It's really two buffers around two overlapping points of red-blue intersection. Thus, this diagram displays five buffers in all.)

5. Use the AddField tool to create four new fields in each of these clipped layers: [X0], [Y0], [X1], and [Y1]. They will hold point coordinates, so make them doubles and give them lots of precision.

6. Calculate Geometry (invoked by right-clicking on each new field header) enables you to compute the x- and y- coordinates of the start and end points of each clipped polyline: put these into [X0], [Y0], [X1], and [Y1], respectively. This is done for each clipped layer, so 8 calculations are needed.

7. Use the AddField tool to create a new [Angle] field in the intersection point layer.

8. Join the clipped tables to the intersection point table based on common object identifiers. (Joins are performed by right-clicking on the layer name and selecting "Joins and Relates".)

   At this point the point intersection table has 9 new fields: two are named [X0], etc., and one is named [Angle]. Alias the [X0], [Y0], [X1], and [Y1] fields which belong to one of the joined tables. Let's call these (say) "X0a", "Y0a", "X1a", and "Y1a".

9. Use the Field Calculator to compute the angle in the intersection table. Here's a Python code block for the calculation:

   dx = !x1!-!x0!
   dy = !y1!-!y0!
   dxa = !x1a!-!x0a!
   dya = !y1a!-!y0a!
   r = math.sqrt(math.pow(dx,2) + math.pow(dy,2))
   ra = math.sqrt(math.pow(dxa,2) + math.pow(dya,2))
   c = math.asin(abs((dx*dya - dy*dxa))/(r*ra)) / math.pi * 180
   

   The field calculation expression is, of course, merely

   c
   

Despite the length of this code block, the math is simple: (dx,dy) is a direction vector for the first polyline and (dxa,dya) is a direction vector for the second. Their lengths, r and ra (computed via the Pythagorean Theorem), are used to normalize them to unit vectors. (There ought to be no problem with zero lengths, because clipping should produce features of positive length.) The size of their wedge product dxdya - dydxa (after division by r and ra) is the sine of the angle. (Using the wedge product rather than the usual inner product should provide better numerical precision for near-zero angles.) Finally, the angle is converted from radians to degrees. The result will lie between 0 and 90. Note the avoidance of trigonometry until the very end: this approach tends to produce reliable and easily computed results.

Some points may appear multiple times in the intersection layer. If so, they will get multiple angles associated with them.

Buffering and clipping in this solution are relatively expensive (steps 3 and 4): you don't want to do it this way when millions of intersection points are involved. I have recommended it because (a) it simplifies the process of finding two successive points along each polyline within the neighborhood of its intersection point and (b) buffering is so basic it is easy to do in any GIS--no additional licensing is needed above the basic ArcMap level--and usually produces correct results. (Other "geoprocessing" operations might not be so reliable.)

Answer 2

I believe you need to create python script.

You can do it using geoprocessing tools and arcpy.

Here is the main tools and ideas that can be useful for you:

1. Make intersection of your two polyline (lets call them PLINE_FC1, PLINE_FC2) featureclasses (you need point features as a result - POINT_FC) using tool Intersect. You will have IDs from PLINE_FC1, PLINE_FC2 in points POINT_FC.
2. Split PLINE_FC1 by POINT_FC using tool Split Line At Point. By result you will have splitted polylines - main advantage of it that you can just take first/last vertex of such line compare it to the next/previous vertex (coordinate difference) and calculate angle. So, you will have angle of your line in the point of intersect. There is one problem here - you have to run this tool manually several times to realize how output is written. I mean if it takes polyline, split it, write two result polylines to the output and then proceed to the next polyline and repeat. Or may be this part (result of splitting) are written to different memory feature classes, then appended to the output. This is the main problem - to realize how output is written to be able to filter only first part of each polyline after splitting. Another possible solution is to loop throught all result splitted polylines with SearchCursor and take only first encountered (by ID of source polylines PLINE_FC1).
3. In order to get angle you will need to access result polyline's verteces using arcpy. Write resulting angles to points (POINT_FC).
4. Repeat steps 2-3 for PLINE_FC2.
5. Substract angle attributes (in POINT_FC) and get result.

May be it will be very difficult to code step 2 (also some tools require ArcInfo license). Then you can also try to analyse verteces of every polyline (grouping them by ID after intersection).

Here is the way to do it:

1. Take first intersection point POINT_FC. Get its coordinates (point_x, point_y)
2. By its ID take respective source polyline from PLINE_FC1.
3. Take first (vert0_x, vert0_y) and second (vert1_x, vert1_y) verteces of it.
4. For the first vertex calculate tangent of line between this vertex and intersection point: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
5. Calculate same thing for second vertex: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
6. If tan1 is equal to tan2, then you have found two verteces of your line which have intersection point in between and you can calculate angle of intersection for this line. Otherwise you have to proceed to the next pair of verteces (second, third) and so on.
7. Repeat steps 1-6 for every intersection point.
8. Repeat steps 1-7 for second polyline featureclass PLINE_FC2.
9. Substract angle attributes from PLINE_FC1 and PLINE_FC2 and get result.

Answer 3

Recently I was trying to do it on my own.

My clue feature is based on circular points around the intersections of lines as well as points located at one-meter distance from intersections. The output is a polyline feature class that has attributes of angles' number on intersections and angle.

Note that lines should be planarized in order to find intersections and spatial reference has to be set with correct line length display (mine is WGS_1984_Web_Mercator_Auxiliary_Sphere).

Running in ArcMap console but easily can be turned to a script in toolbox. This script uses only line layer in TOC, nothing more.

import arcpy
import time

mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame


line = ' * YOUR POLYLINE FEATURE LAYER * ' # paste the name of line layer here    

def crossing_cors(line_layer):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = arcpy.Describe(line_layer).spatialReference
    
    dict_cors = {}
    dang_list = []

    with arcpy.da.UpdateCursor(line_layer, ['SHAPE@', 'OID@']) as uc:
        for row in uc:
            if row[0] is None:
                uc.deleteRow()

    with arcpy.da.UpdateCursor(line_layer, 'SHAPE@', spatial_reference = sr) as uc:
        for row in uc:
            line = row[0].getPart(0)
            for cor in line:
                coord = (cor.X, cor.Y)
                try:
                    dict_cors[coord] += 1
                except:
                    dict_cors[coord] = 1
    cors_only = [f for f in dict_cors if dict_cors[f]!=1]
    cors_layer = arcpy.CreateFeatureclass_management('in_memory', 'cross_pnt', "POINT", spatial_reference = sr)
    arcpy.AddField_management(cors_layer[0], 'ANGLE_NUM', 'LONG')
    with arcpy.da.InsertCursor(cors_layer[0], ['SHAPE@', 'ANGLE_NUM']) as ic:
        for x in cors_only:
            pnt_geom = arcpy.PointGeometry(arcpy.Point(x[0], x[1]), sr)
            ic.insertRow([pnt_geom, dict_cors[x]])
    return cors_layer

def one_meter_dist(line_layer):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = arcpy.Describe(line_layer).spatialReference
    
    dict_cors = {}
    dang_list = []
    cors_list = []
    with arcpy.da.UpdateCursor(line_layer, 'SHAPE@', spatial_reference = sr) as uc:
        for row in uc:
            line = row[0]
            length_line = line.length 
            if length_line > 2.0:
                pnt1 = line.positionAlongLine(1.0)
                pnt2 = line.positionAlongLine(length_line - 1.0)
                cors_list.append(pnt1)
                cors_list.append(pnt2)
            else:
                pnt = line.positionAlongLine(0.5, True)
    cors_layer = arcpy.CreateFeatureclass_management('in_memory', 'cross_one_meter', "POINT", spatial_reference = sr)
    ic = arcpy.da.InsertCursor(cors_layer[0], 'SHAPE@')
    for x in cors_list:
        ic.insertRow([x])
    return cors_layer

def circles(pnts):

    import math
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.env.overwriteOutput = True
    sr = df.spatialReference
    
    circle_layer = arcpy.CreateFeatureclass_management('in_memory', 'circles', "POINT", spatial_reference = sr)


    ic = arcpy.da.InsertCursor(circle_layer[0], 'SHAPE@')
    with arcpy.da.SearchCursor(pnts, 'SHAPE@', spatial_reference = sr) as sc:
        for row in sc:
            fp = row[0].centroid
            list_circle =[]
            for i in xrange(0,36):
                an = math.radians(i * 10)
                np_x = fp.X + (1* math.sin(an))
                np_y = fp.Y + (1* math.cos(an))
                pnt_new = arcpy.PointGeometry(arcpy.Point(np_x,np_y), sr)
                
                ic.insertRow([pnt_new])
    del ic 
    return circle_layer

def angles(centers, pnts, rnd):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    sr = df.spatialReference

    line_lyr = arcpy.CreateFeatureclass_management('in_memory', 'line_angles', "POLYLINE", spatial_reference = sr)
    arcpy.AddField_management(line_lyr[0], 'ANGLE', "DOUBLE")
    arcpy.AddField_management(line_lyr[0], 'ANGLE_COUNT', "LONG")

    ic = arcpy.da.InsertCursor(line_lyr[0], ['SHAPE@', 'ANGLE', 'ANGLE_COUNT'])

    arcpy.AddField_management(pnts, 'ID_CENT', "LONG")
    arcpy.AddField_management(pnts, 'CENT_X', "DOUBLE")
    arcpy.AddField_management(pnts, 'CENT_Y', "DOUBLE")
    arcpy.Near_analysis(pnts, centers,'',"LOCATION") 

    with arcpy.da.UpdateCursor(line, ['SHAPE@', 'OID@']) as uc:
        for row in uc:
            if row[0] is None:
                uc.deleteRow()

    with arcpy.da.UpdateCursor(pnts, [u'ID_CENT', u'CENT_X', u'CENT_Y', u'NEAR_FID', u'NEAR_DIST', u'NEAR_X', u'NEAR_Y'], spatial_reference = sr) as uc:
        for row in uc:
            row[0] = row[3]
            row[1] = row[5]
            row[2] = row[6]
            uc.updateRow(row)
            if row[4] > 1.1:
                uc.deleteRow()


    arcpy.Near_analysis(pnts, rnd,'',"LOCATION")     
 
    list_id_cent = []
    with arcpy.da.UpdateCursor(pnts, [u'ID_CENT', u'CENT_X', u'CENT_Y', u'NEAR_FID', u'NEAR_DIST', u'NEAR_X', u'NEAR_Y', 'SHAPE@'], spatial_reference = sr) as uc:
        for row in uc:
            pnt_init = (row[-1].centroid.X, row[-1].centroid.Y)
            list_id_cent.append([(row[1], row[2]), row[3], pnt_init])

    list_id_cent.sort()
    values = set(map(lambda x:x[0], list_id_cent))
    newlist = [[y for y in list_id_cent if y[0]==x] for x in values]

    dict_cent_angle = {}

    for comp in newlist:
        dict_ang = {}
        for i, val in enumerate(comp):

            curr_pnt = comp[i][2]
            prev_p = comp[i-1][2]
            init_p = comp[i][0]


            angle_prev = math.degrees(math.atan2(prev_p[1]-init_p[1], prev_p[0]-init_p[0]))
            angle_next = math.degrees(math.atan2(curr_pnt[1]-init_p[1], curr_pnt[0]-init_p[0]))

            diff = abs(angle_next-angle_prev)%180
            

            vec1 = [(curr_pnt[0] - init_p[0]), (curr_pnt[1] - init_p[1])]
            vec2 = [(prev_p[0] - init_p[0]), (prev_p[1] - init_p[1])]

            ab = (vec1[0] * vec2[0]) + (vec1[1] * vec2[1]) 
            mod_ab = math.sqrt(math.pow(vec1[0], 2) + math.pow(vec1[1], 2)) * math.sqrt(math.pow(vec2[0], 2) + math.pow(vec2[1], 2))
            cos_a = round(ab/mod_ab, 2)

            diff = math.degrees(math.acos(cos_a))

            pnt1 = arcpy.Point(prev_p[0], prev_p[1])
            pnt2 = arcpy.Point(init_p[0], init_p[1])
            pnt3 = arcpy.Point(curr_pnt[0], curr_pnt[1])
           

            line_ar = arcpy.Array([pnt1, pnt2, pnt3])
            line_geom = arcpy.Polyline(line_ar, sr)

            ic.insertRow([line_geom , diff, len(comp)])
    del ic

    lyr_lst = [f.name for f in arcpy.mapping.ListLayers(mxd)]
    if 'line_angles' not in lyr_lst:
        arcpy.mapping.AddLayer(df, arcpy.mapping.Layer(line_lyr[0]))


centers = crossing_cors(line)

pnts = one_meter_dist(line)

rnd = circles(centers)

angle_dict = angles(centers, pnts, rnd)