arcpy function to calculate angle between two lines

[CrispinCooper]

Does this exist?

Thanks, C.

[ChrisMathers]

Not that I know of but you could get the geometry object to each, find their angles and calculate it that way. That would be a kind of  brute force way of doing it though.

[CrispinCooper]

Hi, thanks for the replies.

My second best option is to get the geometry for each and calculate manually.  Some variant on atan(dy1/dx1)-atan(dy2-dx2) will work for the planar case I imagine.  I wondered if arc would automatically handle spherical projections, great circles etc though - I'm not sure how to do that! (or is it simpler than I think?)

I'll try near_analysis (though note it isn't present in arcview which which I may need this script to be compatible).  Also am I right in thinking it requires that the features be (i) in a table and (ii) that that table doesn't already have a NEAR_ANGLE field?  I need to calculate multiple angles per feature (e.g. between each feature and n other objects)

[CrispinCooper]

I'm now using brute force...

def angle(pointlist):
    # calculates using cosine rule: A^2 = B^2 + C^2 - 2BC cos a
    # a = arccos ((B^2+C^2-A^2)/(2BC))
    # but multiply by 180/pi to convert from radians
    A_vector = pointlist[0] - pointlist[2]
    B_vector = pointlist[0] - pointlist[1]
    C_vector = pointlist[1] - pointlist[2]
    A2,B2,C2 = map(square_length_of_vector,(A_vector,B_vector,C_vector))
    B,C = map(numpy.sqrt,(B2,C2))
    side_ratio = (B2+C2-A2)/(2*B*C)
    if side_ratio > 1: # possible on sharp angle due to rounding errors
            return 0.
    elif side_ratio < -1: # possible on blunt angle due to rounding errors
            return 180.
    else:
            return 180/numpy.pi * numpy.arccos(side_ratio)

...shame if arc doesn't have a way of doing that that generalises to spherical geometry, though.

[ChrisSnyder · Accepted solution]

You can make rapid coordinate system conversions by using an appropriate spatial reference parameter in a cursor. For example, this code spits out the vertice coordinates of a township layer in WGS84 coordinates (not their native StatePlane coordinates).

import arcgisscripting
gp = arcgisscripting.create(9.3)
fc = r"\\snarf\am\div_lm\ds\gis\tools\sde_connections_read\ropa_gis_layer_user.sde\CADASTRE.TOWNSHIP_SV"
spatialRef = r"C:\Program Files\ArcGIS\Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
dsc = gp.describe(fc)
shapeFieldName = dsc.shapeFieldName
oidFieldName = dsc.oidfieldname
searchRows = gp.searchcursor(fc, "", spatialRef)
searchRow = searchRows.next()
while searchRow:
    oidFieldValue = searchRow.getvalue(oidFieldName)
    print oidFieldName + " = " + str(oidFieldValue)    
    shapeFieldObj = searchRow.getvalue(shapeFieldName)
    partCount = shapeFieldObj.partcount
    partIndex = 0
    while partIndex < partCount:
        print "Part " + str(partIndex)+ ":"
        partObj = shapeFieldObj.GetPart(partIndex)
        pointObj = partObj.next()
        pointCounter = 0
        while pointObj:
            print pointObj.x, pointObj.y
            pointObj = partObj.next()
            pointCounter += 1
            # If pointObj is null, either the part is finished or there is an interior ring
            if not pointObj: 
                pointObj = partObj.next()
                if pointObj:
                    print "Interior Ring:"
        partIndex += 1
searchRow = searchRows.next()

[CrispinCooper]

You can make rapid coordinate system conversions by using an appropriate spatial reference parameter in a cursor.

Ah, that's probably what I want then - thank you!

[DavidStrip]

If you are going to be computing this a lot and efficiency matters, you can do slightly better by observing that
   [INDENT]a · b = |a| |b| cos(theta)[/INDENT]
where · is the dot product, |a| is the length of a and theta is that angle between vectors a and b.
Rearrange and you get
[INDENT]theta = arccos((a · b)/(|a| |b|))[/INDENT]