2

I'm looking for an automated way to select line features that have specified angle between vertexes. That in order to check water mains layer -steel pipes can have only 45 or 90 degrees turns and I need to check the digitizer work. ArcMap Data Reviewer have a tool for that, but i don't have this extension.

  • 3
    Do you have some python code that you have started to write? The solution is likely going to be a search cursor to loop over each geometry, a loop over each vertex, and the cosine of the dot product of each line segment to find the angle. – Jay Laura Jul 21 '13 at 18:24
4
+50

Something like the following should get you started at least. I think this is to what Jay Laura was referring. If anyone wants a better explanation of a particular part, please ask.

from math import acos
from numpy.linalg import norm
import arcpy

shpin = arcpy.GetParameterAsText(0)

def findangle(p1,p2,p3):
    """This function returns the angle p1p2p3 given tuple inputs of the coordinates."""

    #The vectors from p2->p1 and p2->p3
    v1 = [p1[0]-p2[0], p1[1]-p2[1]]
    v2 = [p3[0]-p2[0], p3[1]-p2[1]]

    #The sum of the dotproducts of each vector.    
    dotproduct = sum(p*q for p,q in zip(v1,v2))

    #The multiplication of the magnitudes of each vector.
    magnitude = norm(v1) * norm(v2)

    #Rounding here to account for number of decimal places in rads.
    return round(acos(dotproduct/magnitude), 5)

#Write OID and coordinates of every vertex to list of lists.
#OID field is used to break up coordinates by feature.
coords = [[row[0], row[1].firstPoint.X, row[1].firstPoint.Y] for row in
          arcpy.da.SearchCursor(shpin, ["OID@", "SHAPE@"], explode_to_points=True)]

#Get total count of features in input.
entries = int(arcpy.GetCount_management(shpin).getOutput(0))

#pi/4 and pi/2 radians, respectively (45 and 90 degrees)
rads = [0.78534, 1.57080]
angles = []

arcpy.SetProgressor("step", "Determining angles", 0, entries-1, 1)
for i in xrange(entries):
    arcpy.SetProgressorPosition(i)
    #Filter out segments not equal to current OID.
    featcoords = filter(lambda f: f[0] == i, coords)

    if len(featcoords) >= 3: #Only check segments with 3 or more vertices.        
        j = 0        
        while j < len(featcoords) - 2: 
            try: #Sometimes a math domain error is thrown.
                #If the angle is not 45 or 90 degrees, append the current OID.
                if findangle(featcoords[j][1:], featcoords[j+1][1:], featcoords[j+2][1:]) not in rads:
                    angles.append(i)
            except:
                pass
            j += 1

#Remove uniques as some segments may have multiple offenders.
arcpy.AddMessage(list(set(angles)))
|improve this answer|||||
  • +1 It is. I figured that if the OP wanted the tool they would respond to the question and/or provide some code. I would suggest that for a use everyday tool via ArcMap that it would be trivial to pack this logic as a Python add-in and let the tool select all the necessary OIDs (or FIDs) via a single click. – Jay Laura Jul 24 '13 at 4:20
  • Thanks for the education Paul and Jay. I'm not sure if I would have found this option on my own. This will come in handy for sure. – Dror Har Gil Jul 25 '13 at 13:27
1

It's not in python, but ET Geowizard (http://www.ian-ko.com/) can do it. Using the conversion tool (from Polylines to Points), you can have the angle at each vertex. I don't know if it's the bearing of the segment or the angle of the segment compared with the previous one, but then, it's easy enough to calculate.

|improve this answer|||||
  • Thanks fgCartigraphix, but ET GeoWizard tool gives the vertex angle in relation to the start or to the end of line. What i need is to compare each vertex angle to the previous vertex and to the next vertex- in order to figure his turning angle. – Dror Har Gil Jul 21 '13 at 5:17
  • But is line1 (vertex1 to vertex2) is 34° and line2 (vertex2 to vertex3) is 124°, the difference is 90°. It can be quite easy to calculate, no? – fgcartographix Jul 25 '13 at 19:55
0

The Split Line at Verticies (Data Management) tool may be of use:

http://resources.arcgis.com/en/help/main/10.1/index.html#//00170000003z000000

Then use the ET GeoWizard tool to obtain the angles you're after.

|improve this answer|||||
  • Thanks Rayner for your comment, but I'm looking for a fast and simple tool i can use on daily base on my geometric network without corrupt it. In addition, your proposal can provide the azimuth of the vertex or the angle in relation to start or end of line - and not the turning angle of the line at the vertex. Geog – Dror Har Gil Jul 23 '13 at 4:59

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.