import arcgisscripting import math, os, sys, traceback def rowgen(rows): row = rows.next() while row: yield row row = rows.next() def partgen(parts): numparts = parts.PartCount for numpart in range(numparts): part = parts.getPart(numpart) yield part def calculate_bend(pnt0, pnt1, pnt2): # ------------------------------------------------------------------------------ # Calculate angle at middle vertex between three points. # Negative bend is right bend # Positive bend is left bend # Algorithim by Don Sawatzky # ------------------------------------------------------------------------------ #Get trig angle of first line delx1 = pnt1.X - pnt0.X dely1 = pnt1.Y - pnt0.Y angle1 = math.atan2(dely1, delx1) #Get trig angle of next line delx2 = pnt2.X - pnt1.X dely2 = pnt2.Y - pnt1.Y angle2 = math.atan2(dely2, delx2) #Convert del angle to degrees bend = math.degrees(angle2 - angle1) if abs(bend) > 180: if bend > 0: bend -= 360 else: bend +=360 return bend def calculate_distance(pt1, pt2): # ------------------------------------------------------------------------------ # Calculate the distance between two points. # ------------------------------------------------------------------------------ dx = abs(pt1.X - pt2.X); dy = abs(pt1.Y - pt2.Y); distance = math.sqrt((dx*dx) + (dy*dy)) return distance def CharacterizePolylines(polylineIn, polylinesOut, targetLength, minAngle, maxAngle): # ------------------------------------------------------------------------------ # Calculates various statistics for length based subsets of polylines. # ------------------------------------------------------------------------------ try: #Create output polylines feature dataset outdir = os.path.dirname(polylinesOut) outbase = os.path.basename(polylinesOut) spatref = gp.describe(polylineIn).spatialreference gp.CreateFeatureClass_management(outdir, outbase, 'POLYLINE', '#', '#', 'ENABLED', spatref) gp.AddField_management(polylinesOut, "ANGLE_PREF", "FLOAT") gp.AddField_management(polylinesOut, "ANGLE_TOT", "FLOAT") gp.AddField_management(polylinesOut, "ANGLE_POS", "FLOAT") gp.AddField_management(polylinesOut, "ANGLE_NEG", "FLOAT") gp.AddField_management(polylinesOut, "LEN_TOT", "FLOAT") gp.AddField_management(polylinesOut, "LEN_STRT", "FLOAT") gp.AddField_management(polylinesOut, "LEN_TIGHT", "FLOAT") gp.AddField_management(polylinesOut, "STRT_RATIO", "FLOAT") gp.AddField_management(polylinesOut, "ANG_RATIO", "FLOAT") lineArray = gp.CreateObject("Array") min_angle = float(minAngle) max_angle = float(maxAngle) target_length = float(targetLength) line_num = 1 #Open output feature class outrows = gp.insertcursor(polylinesOut) # # ----- Loop through each input polyline # for inrow in rowgen(gp.searchcursor(polylineIn)): gp.addmessage("Processing polyline: %d"%line_num) line_num = line_num + 1 pline = inrow.Shape #Insert output row for linepart in partgen(pline): #gp.addmessage("No. pnts: %d"%linepart.Count) #Process output features for pt_num in range(linepart.Count - 1): lineArray.RemoveAll() pos_angle = 0.0 neg_angle = 0.0 offset_to_next = 1 length_total = 0.0 length_straight = 0.0 length_tight = 0.0 in_straight_part = False in_tight_part = False # # ----- loop until the end-of-line or the subset length is > target length # while (pt_num + offset_to_next < linepart.Count and length_total < target_length): pt1 = linepart.GetObject(pt_num + offset_to_next - 1) pt2 = linepart.GetObject(pt_num + offset_to_next) length_segment = calculate_distance(pt1, pt2) length_total = length_total + length_segment # # ----- manage angles when we have 3 or more points # if (offset_to_next > 1): lineArray.add(pt2) cur_angle = calculate_bend(ptLast, pt1, pt2) # # ----- keep track of positive and negative turns # if ((abs(cur_angle) >= min_angle) and (abs(cur_angle) <= max_angle)): if (cur_angle < 0.0): neg_angle = neg_angle + cur_angle else: pos_angle = pos_angle + cur_angle # # ----- keep track of the straight lengths # if (abs(cur_angle) < min_angle): if (not in_straight_part): length_straight = length_straight + calculate_distance(ptLast, pt1) length_straight = length_straight + length_segment in_straight_part = True else: in_straight_part = False # # ----- keep track of the tight lengths # if (abs(cur_angle) > max_angle): if (not in_tight_part): length_tight = length_tight + calculate_distance(ptLast, pt1) length_tight = length_tight + length_segment in_tight_part = True else: in_tight_part = False else: lineArray.add(pt1) lineArray.add(pt2) offset_to_next = offset_to_next + 1 ptLast = pt1 # # ----- create output row and assign values # if (length_total >= target_length): newrow = outrows.NewRow() newrow.ANGLE_POS = pos_angle newrow.ANGLE_NEG = neg_angle newrow.ANGLE_TOT = pos_angle + abs(neg_angle) newrow.LEN_TOT = length_total newrow.LEN_STRT = length_straight newrow.LEN_TIGHT = length_tight newrow.STRT_RATIO = length_straight / length_total * 100.0 if (pos_angle > abs(neg_angle)): newrow.ANGLE_PREF = pos_angle if ((pos_angle + abs(neg_angle)) * 100.0 > 0): newrow.ANG_RATIO = pos_angle / (pos_angle + abs(neg_angle)) * 100.0 else: newrow.ANGLE_PREF = abs(neg_angle) if ((pos_angle + abs(neg_angle)) * 100.0 > 0): newrow.ANG_RATIO = abs(neg_angle) / (pos_angle + abs(neg_angle)) * 100.0 newrow.shape = lineArray newrow.ID = inrow.FID outrows.InsertRow(newrow) del outrows #Like closes the output db file except Exception, msg: # get the traceback object tb = sys.exc_info()[2] # tbinfo contains the line number that the code failed on and the code from that line tbinfo = traceback.format_tb(tb)[0] # concatenate information together concerning the error into a message string pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n " + \ str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n" # generate a message string for any geoprocessing tool errors msgs = "GP ERRORS:\n" + gp.GetMessages(2) + "\n" gp.AddError(msgs) # return gp messages for use with a script tool gp.AddError(pymsg) # print messages for use in Python/PythonWin print pymsg print msgs if __name__ == "__main__": #create geoprocessor gp = arcgisscripting.create() #Get input and output dataset names polylineIn = gp.getparameterastext(0) desc = gp.describe(polylineIn) #gp.addmessage(desc.shapetype) assert desc.shapetype == 'Polyline', 'Input features: %s are not Polyline type' % polylineIn polylinesOut = gp.getparameterastext(1) targetLength = gp.getparameterastext(2) minAngle = gp.getparameterastext(3) maxAngle = gp.getparameterastext(4) CharacterizePolylines(polylineIn, polylinesOut, targetLength, minAngle, maxAngle) gp.SetParameterAsText(1, polylinesOut)