How to identify polygons with (approximate) right-angles?

[HayleyGlover]

Hi, I have a shapefile with thousands of polygons of all sorts of shapes, and am trying to extract only those which are roughly rectilinear (more or less straight edged, with angles that are close to 90 degrees, and have four sides) and eliminate all the round/irregular polygons. I've attached an image of some of the polygons I've been generating.

Is there a way to automatically identify polygons in ArcMap which have interior angles of (for example) 85-95 degrees? Most of what I have seen only looks at finding the angle (ie. compass bearing) of an overall polygon, which is not what I want.

Thank you so much in advance!

[DanPatterson_Retired · Accepted solution]

as a general function tool with a single line code change you can use

# -*- coding: utf-8 -*-
import numpy as np

def angles_poly(a, inside=True, in_deg=True):
    """Sequential angles from a poly* shape
     NOTE: this line.... angle = np.sum(angles)
           can be changed to `np.min`, `np.max`, np.ptp or others
            depending on what needs to be returned
    """
    a = a.getPart()
    a =np.asarray([[i.X, i.Y] for j in a for i in j])
    if len(a) < 2:
        return None
    elif len(a) == 2:  # **** check
        ba = a[1] - a[0]
        return np.arctan2(*ba[::-1])
    else:
        angles = []
        if np.allclose(a[0], a[-1]):  # closed loop
            a = a[:-1]
            r = (-1,) + tuple(range(len(a))) + (0,)
        else:
            r = tuple(range(len(a)))
        for i in range(len(r)-2):
            p0, p1, p2 = a[r[i]], a[r[i+1]], a[r[i+2]]
            ba = p1 - p0
            bc = p1 - p2
            cr = np.cross(ba, bc)
            dt = np.dot(ba, bc)
            ang = np.arctan2(np.linalg.norm(cr), dt)
            angles.append(ang)
    if in_deg:
        angles = np.degrees(angles)  #---- function to change ------
    angle = np.sum(angles)
    return angle
#__esri_field_calculator_splitter__
#angles_poly(!Shape!)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

change line 34 to

angles = np.ptp(np.abs(angles))‍‍

now if you run this field calculation a square would have a deviation of 0 degrees, an equilateral triangle 0 etc.

So you might have to try a couple of other field calculations

angle = np.min(angles)
angle = np.max(angles
etc
‍‍‍

[HayleyGlover]

Thank you so much Dan, this is just what I was looking for.

[JoshuaBixby]

trying to extract only those which are roughly rectilinear (more or less straight edged, with angles that are close to 90 degrees, and have four sides)

What is described in parentheses is rectangular, which is a special case of rectilinear.  Are you interested in finding "roughly rectilinear" or "roughly rectangular?"

[HayleyGlover]

Oh sorry - yes rectangular especially.

[DarrenWiens2]

Without specifically coding each requirement, you could look into running Minimum Bounding Geometry (RECTANGLE_BY_AREA) to create a rectangle for each feature, and select those meeting some threshold (e.g. 99% of area). The more rectangular the feature, the closer it will match the minimum rectangle.

[HayleyGlover]

Thanks Darren - I'd been wondering if that function could be useful for me.