Calculate Polygon Angles at Vertices

[PeterMetzger1]

When I started this job someone had already had a contractor survey a bunch of our buildings at our installations and look at our cad drawings to make building footprints for real property management purposes.  I've been aligning floorplans to the these footprints, and much to my chagrin I find that a lot of these buildings that were entered are not square (apparently they thought maybe we buildt structures with 95° angles, needless to say we don't).  So now I'm writing a script to go through and calculate the angles at each vertex of the polygon so we can find out which buildings need to be resurveyed.  Right now I have the script so that it will calculate the vector of each line from the X axis.  But, what I need it to do is select the lines that interesect the vertex (the polygons are converted to polylines and then exploded, and points are created at the vertices) and I think I have that part down, but now I need to add the selected two lines vectors together to get the angle and put it in the field for angles at that vertex.  I then need to step through the process for each of the thousand some vertices with either a for or while loop.  I'm new to python but I think I can figure out the loop, but doing the math with the selected features is confusing me (actually any help with the loop would be nice too).

Thanks in advance, here's the code I have.

#
# This Script Validates if the surveyed building footprints are square enough to be acceptable#
#
#

                            # Import system modules
import arcpy
import math
from arcpy import env

# Set environment settings
arcpy.env.workspace = "C:\Users\peter.a.metzger\Desktop\Projects\CGN_Building_survey_validating\Building_Validate.gdb"
inFeatures = "Building_footprint"

#Turn Polygon into Polyline in scratch workspace
arcpy.PolygonToLine_management (inFeatures, "polygon2polyline")

#explode polyline into individual lines
arcpy.SplitLine_management ("polygon2polyline",  "explodedLine")

#convert polygon vertices to points
arcpy.FeatureVerticesToPoints_management (inFeatures, "Poly2Point")


# Set local variables

fieldName = "Angle"
expression = "GetAzimuthPolyline(!Shape!)"
codeblock = """def GetAzimuthPolyline(shape):
    radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
    degrees = radian * 180 / math.pi
    return degrees"""

# Execute AddField
arcpy.AddField_management("explodedLine", fieldName, "Double", 6)
arcpy.AddField_management("Poly2Point", fieldName, "Double", 6)

# Execute CalculateField 
arcpy.CalculateField_management("explodedLine", fieldName, expression, "PYTHON_9.3", codeblock)

# Make Poly2Point and explodedLine a layer file so that select by attribute can read it

arcpy.MakeFeatureLayer_management("Poly2Point", "Poly2Point_lyr")

arcpy.MakeFeatureLayer_management("explodedLine", "explodedLine_lyr")

#Step through each OBJECTID and select the polygons that share the point and subtract their values
##Select the vertex for the angle##

#NEED TO CHANGE TO A FOR LOOP#
#THIS IS WHERE I'M STUCK#
#NEED TO BE ABLE TO SUBTRACT THE TWO SELECTED LINES FROM ONE ANOTHER#
obid=1
arcpy.SelectLayerByAttribute_management ("Poly2Point_lyr", "NEW_SELECTION", "[OBJECTID] = 'obid'")

#select the two exploded lines from which to calculate the vertex angle

arcpy.SelectLayerByLocation_management ("explodedLine_lyr", "INTERSECT","Poly2Point_lyr", "", "NEW_SELECTION")

[markdenil]

You could get all the verticies in a building, and cycle through them, one at a time.
Use Near_analysis (in_features, near_features, "", "", ANGLE)
on each (as in_features) and the vertex before and vertex after it as near_features.
Setting the {angle} parameter to ANGLE reports the angle....

You can get the verticies as an array of points from the shape object.

[PeterMetzger1]

Thanks for the input, I thought that would work and would be much simpler, I would just tell it to find the angle between the point and the closest point and then find the angle between the point and then the second closest point.  However, a lot of our buildings are not rectangular and have a lot of things jutting out.  Take for example in the image, I want to find the angle of Vertex A.  It would be done using near to find the angle to point C from A, and B from A.  The problem is that D is closer to A and would calculate the wrong angle.[ATTACH=CONFIG]15973[/ATTACH].

For anyone else reading this post to clarify I am looking to solve for angle C in the following drawing, I have been able to calculate A and B
[ATTACH=CONFIG]15974[/ATTACH]

PS sorry about the MS Paint images

[SolomonPulapkura]

Once you have A and B isn't C = 180-(A+B)?

[markdenil]

I suggested cycling through the array of polygon points, not just getting angles to the closest two points.
The array will have each vertex point in order, running clockwise for an exterior ring and counter-clockwise for an internal ring.
See the page on Reading geometries for more information.

As well, unlike your diagram, the angles to nearest points 1 and 2 would be overlapping angles, not supplemental.
Both are measured from the x axis (east).
It would be the difference between the first and second angles you were looking to find.

Figuring if your calculated angle was internal or external is another issue, but if you are just flagging non-square corners, internal or external would not matter.

[PeterMetzger1]

So I went to the UC last week and had a very nice girl (Thanks Stephanie!) give me a hand on this, and we (mostly she) figured out how to get the angle at the vertices.  After the two lines at the vertex are selected the are exported to a feature class in memory and then a search cursor is used to get their angles, which are then put into a list and then subtracted from one another.  They are then put back into the layer file using an UpdateCursor.  This all worked great... when I was doing it on one vertex.  Once I tried to put it into a loop I keep getting an "IndexError: list index out of range".  Any help on what I have done wrong would be greatly appreciated.

#
#
# This Script Validates if the surveyed building footprints are square enough to be acceptable#
#
#

                            # Import system modules
import arcpy
import math
from arcpy import env

# Set environment settings

arcpy.env.workspace = "in_memory"

inFeatures = "C:/Users/xxxxxxxx/Desktop/Projects/CGN_Building_survey_validating/test.gdb/test_footprint"

#Turn Polygon into Polyline in scratch workspace
arcpy.PolygonToLine_management (inFeatures, "polygon2polyline")

#explode polyline into individual lines
arcpy.SplitLine_management ("polygon2polyline",  "explodedLine")

#convert polygon vertices to points
arcpy.FeatureVerticesToPoints_management (inFeatures, "Poly2Point")


# Set local variables

fieldName = "Angle"
expression = "GetAzimuthPolyline(!Shape!)"
codeblock = """def GetAzimuthPolyline(shape):
    radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
    degrees = radian * 180 / math.pi
    return degrees"""

# Execute AddField
arcpy.AddField_management("explodedLine", fieldName, "Double", 6)
arcpy.AddField_management("Poly2Point", fieldName, "Double", 6)

# Execute CalculateField 
arcpy.CalculateField_management("explodedLine", fieldName, expression, "PYTHON_9.3", codeblock)

# Make Poly2Point and explodedLine a layer file so that select by attribute can read it

arcpy.MakeFeatureLayer_management("Poly2Point", "Poly2Point_lyr")

arcpy.MakeFeatureLayer_management("explodedLine", "explodedLine_lyr")


#Identify how many vertices need to be calculated

numb_vertices = int(arcpy.GetCount_management("Poly2Point_lyr").getOutput(0)) 
print numb_vertices

#Loop through each vertex by OBJECTID

for i in range (0, numb_vertices):
    obj = "OBJECTID=%s" %(i)

    arcpy.SelectLayerByAttribute_management ("Poly2Point_lyr", "NEW_SELECTION", obj)

#select the two exploded lines from which to calculate the vertex angle

    arcpy.SelectLayerByLocation_management ("explodedLine_lyr", "INTERSECT","Poly2Point_lyr", "", "NEW_SELECTION")

#Export those lines to a new feature class

    arcpy.FeatureClassToFeatureClass_conversion ("explodedLine_lyr", "in_memory", "calc_values")

#Create List to do angle math

    Angle_math = []

#Return Angle_Values

    rows = arcpy.SearchCursor ("calc_values","","","Angle")

    current_Angle_Object = ""

#Iterate through the rows in the cursor and add them to the list

    for row in rows:
            if current_Angle_Object !=row.Angle:
                current_Angle_Object = row.Angle
            Angle_math.append(row.Angle)
#Calculate the angle at the vertex
    
    Angle_Final = 180 - Angle_math[0] - Angle_math[1]
    del row
    del rows

    print(Angle_Final)

#add data to the layer file

    rows = arcpy.UpdateCursor("Poly2Point_lyr")
    for row in rows:
        row.Angle = Angle_Final
        rows.updateRow(row)
    del row
    del rows
    del Angle_math
#Export all the final angles in Poly2Point_lyr to a feature class

angle_output_location = "C:/Users/xxxxxxxxx/Desktop/Projects/CGN_Building_survey_validating/test.gdb/Vertex_Angles"
arcpy.CopyFeatures_management("Poly2Point_lyr", angle_output_location)

[MarkMenzel]

Peter,

Did you ever get that script to work on a whole data set?

Thanks,
Mark

[CarlSunderman]

Peter,

Did you ever get that script to work on a whole data set?

Thanks,
Mark

I have some code for calculating angles of lines, you could add to convert polygons to lines, split at vertices and use the following code to calc angles. There could prob be some major improvements on it, bit it works.

with arcpy.da.UpdateCursor(<outFeatureclass>, [<put fields to calc in here>])as rows:
    count = 0
    for row in rows:
        # Print the current line ID
        while count == 0:
            prevm = 0
            prevX = 0
            prevY = 0
            break
        else:
            prevm = m
            prevX = endx
            prevY = endy            
        count = count+1
        print("Feature {0}:".format(row[0]))
        #Set start point
        startpt = row[1].firstPoint
        #Set Start coordinates
        startx = startpt.X
        #print startx
        starty = startpt.Y
        #print starty
        #Set end point
        endpt = row[1].lastPoint
        #Set End coordinates
        endx = endpt.X
        #print endx
        endy = endpt.Y
        #print endy

        # Get slope of line
        radian = math.atan((endx - startx)/(endy - starty))
        m = radian*180 / math.pi
        if (m<=0):
            m+=360.0
        slope = repr(m)
        change = prevm-m
        print slope
        print change

hope this helps and let me know if you have any questions