topic Re: How to identify polygons with (approximate) right-angles? in Spatial Data Science Questions

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

HayleyGlover — Mon, 10 Sep 2018 03:43:31 GMT

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!

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

DanPatterson_Retired — Sun, 12 Dec 2021 07:57:25 GMT

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
‍‍‍

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

JoshuaBixby — Mon, 10 Sep 2018 13:27:52 GMT

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?"

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

DarrenWiens2 — Mon, 10 Sep 2018 17:15:59 GMT

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.

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

HayleyGlover — Mon, 10 Sep 2018 22:33:11 GMT

Oh sorry - yes rectangular especially.

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

HayleyGlover — Mon, 10 Sep 2018 22:33:54 GMT

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

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

HayleyGlover — Mon, 10 Sep 2018 22:34:42 GMT

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