nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

arcpy function to calculate angle between two lines

CrispinCooper

Does this exist?

Thanks, C.

Vind meer berichten getagd met

Accepted answers

ChrisSnyder

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()

All comments

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]

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!

ChrisSnyder

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

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.

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)

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.

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.