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.
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-p2, p1-p2] v2 = [p3-p2, p3-p2] #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, row.firstPoint.X, row.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 == 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)))
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.
The Split Line at Verticies (Data Management) tool may be of use:
Then use the ET GeoWizard tool to obtain the angles you're after.