2

I am trying to create a tool for measuring the angle between two intersecting lines. I have created a buffer and clipped the lines to within this buffer and computed the angle between the lines using the formulas below:

dx = !X1a!-!X0a!

dy = !Y1a!-!Y0a!

dxa= !X1b!-!X0b!

dya= !Y1b!-!Y0b!

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((dxdya - dydxa))/(r*ra)) / math.pi * 180

((source How to find the angle between intersecting features in two featureclasses?, user: whuber))

This returns the angle between the two intersecting lines, but, I need the angle between the lines to always be between line1 counter-clockwise to line2. The tool will be used to calculate the angle between a road intersecting with a rail road.

Does anyone know the formula to do this?

The image below shows a sample of the data I am working with and the angle I would like to calculate. The diameter of the circle (buffer polygon) is 10 centimetres.

Description of the angle I want to calculate

2

There is a function in the math module that is much simpler to use:

deg = math.degrees(math.atan2(dy,dx))

Note the 2 on the end on tan that improves the original range from 0 - 90 degrees. This extended function will provide a 0 - 360 degree bearing for a slope (in radians). There are another pair of functions to convert between degrees and radians math.degrees() and math.radians(). No need to do that yourself.

You can then just subtract the full bearing of each line to get all the angles between two lines. The dy and dx can be negative values.

I can see that you might need to sort the end vertices to find out the direction of the lines. You do not explain why you need a particular angle. Is this because it is related to the existing line direction or do you always want the South facing angle?

All this is too fiddly to do in the expression calculator. The time has come to learn a bit of Python to put all the steps in a script and update the fields using a cursor.

   


#-------------------------------------------------------------------------------
# Name:        intersection_angle
# Purpose:     calculate intersection angles
#              should not be so hard!
# Author:      kimo
#
# Created:     03/08/2016
# Copyright:   (c) kimo 2016
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
# method
#       input
#             railway line, road lines
#       output crossing points with angles, segments with bearings
# should select only those intersecting to remove redundant lines
# 1 intersect railways and roads to get intersection points
# 2 buffer intersection points, retain intersection id
# 3 intersect buffers with original rail and roads, not clip to tag id
#
# 4 for each line
#      use geometry object to find ends and hence bearing
#      write out to two dictionaries
# 5 use dictionaries to calculate angles and store back in points

import math
import arcpy

arcpy.env.workspace = "j:/temp/sample.gdb"
arcpy.env.overwriteOutput = True
# simplified set that have intersecting lines, no junctions at crossings
# single part and generally error free
inFeatures = ["road", "rail"]

# Get intersections as points and add extra fields for output
arcpy.analysis.Intersect(inFeatures, "crossing1", "", 0.01, "point")
print arcpy.GetMessages()

# select only road/rail intersections not road/road or rail/rail
# but ESRI don't 'do' node topology any more
# start by deleting duplicate points
arcpy.DeleteIdentical_management(in_dataset="crossing1", fields="Shape", xy_tolerance="", z_tolerance="0")

# just select intersection points on railway lines, handle railway junctions later
arcpy.management.MakeFeatureLayer('crossing1','cross_lay')
arcpy.management.SelectLayerByLocation('cross_lay','INTERSECT','rail',1)
arcpy.management.CopyFeatures('cross_lay','crossing2')
arcpy.management.AddField("crossing2","id",'INTEGER')
arcpy.management.AddField("crossing2","angle1",'FLOAT')
arcpy.management.AddField("crossing2","angle2",'FLOAT')
arcpy.management.DeleteField("crossing2",['type'])

# Add a serial ID, use a cursor not field calculator if we are in python anyway
# easier to debug, trap errors, understand, just as fast
id = 1
with arcpy.da.UpdateCursor("crossing2",['id']) as cur:
    for row in cur:
        row[0] = id
        cur.updateRow(row)
        id+=1

# Buffer crossings by 10 meters to get local angles
arcpy.analysis.Buffer("crossing2", "crossing10", "10 meters")

# Intersect original merge (not clip because ID attributes lost)
arcpy.analysis.Intersect(["road", "crossing10"], "roadrail_clip1")
print arcpy.GetMessages()
arcpy.analysis.Intersect(["rail", "crossing10"], "roadrail_clip2")
print arcpy.GetMessages()
arcpy.management.Merge(["roadrail_clip1","roadrail_clip2"],"roadrail_clip")
# each road, rail line will inherit crossing id and a type tag

# calculate bearings from ends at buffer intersections
arcpy.management.AddField("roadrail_clip","bearing",'FLOAT')
dCross_road = {}
dCross_rail = {}

with arcpy.da.UpdateCursor("roadrail_clip",['id','bearing','SHAPE@','type']) as cur:
    for row in cur:
        dx = row[2].lastPoint.X - row[2].firstPoint.X
        dy = row[2].lastPoint.Y - row[2].firstPoint.Y
        bearing = math.degrees(math.atan2(dy,dx))
        # leave alone with original sense
        row[1] = bearing
        cur.updateRow(row)
        # to compare bearings by crossing ID
        # need to collect bearings by ID and type
        if row[3] == 'road':
            dCross_road[row[0]] = bearing
        elif row[3] == 'rail':
            dCross_rail[row[0]] = bearing
        else:
            pass
            # no interest

# angles can be calculated from the two dictionaries
# if both entries exist
i = 0
dAngle = {}
for id in dCross_rail:
    if dCross_rail.get(id,None) and dCross_road.get(id,None):
        b_rail = dCross_rail[id]
        b_road = dCross_road[id]
        # to decide which angle is required (anticlockwise from the road)
        # but still have to allow for two parts of the rail line
        if b_rail 

railway bearing reversed

  • Thanks Kimo. I'm not great at Python, I use it sometimes at work but that is all. To answer your question, I need to use the first section of rail road counter clockwise around the circle because that is the standard format for the people I am compiling the data for. They always want the angle between the road and the rail road to its right regardless of road or rail direction. I may just end up calculating the angles using the COGO report feature in the editor if I can't figure out how to do this. I tried some nested 'if' statements, but, couldn't get it to fully work. – jdavid05 Aug 3 '16 at 17:18
  • Ok, it is easy to get the required angle. Just fold all bearings to be positive to ignore the line direction. Then the subtraction of the road from rail always gives you the right answer. I will modify the script above. – kimo Aug 3 '16 at 21:11
  • Another tweak was required after I checked a few hundred real cases. If the lines have kinks inside the 10 m buffer then the result may be misleading. I also did not allow for the rail has two ends so a new test has allowed for the case when the rail bearing is less than the road bearing. Also 1.5 metre intersection results in offset intersections so I reduced the tolerance to 0.01 metres. – kimo Aug 3 '16 at 23:34
  • Thanks! Looks like this took a lot of work. This was really helpful. – jdavid05 Aug 4 '16 at 13:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.