Discriminate convex vs concave vertices within each polygon feature

781
2
Jump to solution
12-05-2022 07:40 AM
Labels (2)
teakplantation
New Contributor III

Using the calculate geometry tool, one can calculate the number of vertices for each polygon in a feature class (https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/calculate-geometry-attribute....

A vertex can be concave or convex, depending on whether the inner angle at the vertex is larger or smaller than 180 degrees.

Does anybody know of a tool which can calculate how many vertices per polygon are convex and how many are concave?

 

Tags (1)
0 Kudos
1 Solution

Accepted Solutions
JohannesLindner
MVP Frequent Contributor

Starting from Dan's suggestion of convexHull(), which seemed good:

with arcpy.da.UpdateCursor("TestPolygons", ["SHAPE@", "IntegerField"]) as cursor:
    for shp, n in cursor:
        n = len(shp[0]) - len(shp.convexHull()[0])  # only for singlepart polygons
        cursor.updateRow([shp, n])

 

This works well for simple geometries:

JohannesLindner_0-1670325926043.png

 

But it starts to break down for more complex geometries, where convex vertices are not always on the convex hull:

JohannesLindner_1-1670326259996.png

 

I was ready to start fussing with angles between polylines when I realized that you can quite simply say whether a vertex is concave or convex: cretae a triangle of the vertex and its neighbours. If that triangle is inside the polygon, the vertex is convex, else it's concave.

JohannesLindner_2-1670329512390.png

 

def get_concave_vertex_count(polygon):
    n = 0
    for part in polygon:
        for i, vertex in enumerate(part):
            prev_vertex = part[i-1]
            try:
                next_vertex = part[i+1]
            except IndexError:
                next_vertex = part[0]
            triangle = arcpy.Polygon(arcpy.Array([prev_vertex, vertex, next_vertex]))
            n += triangle.touches(polygon)  # touch: only boundaries intersect
    return n


with arcpy.da.UpdateCursor("PolygonLayer", ["SHAPE@", "ConcaveVertices"]) as cursor:
    for shp, n in cursor:
        n = get_concave_vertex_count(shp)
        cursor.updateRow([shp, n])

Have a great day!
Johannes

View solution in original post

2 Replies
DanPatterson
MVP Esteemed Contributor

You could return the difference between the polygon and its convex hull (symmetric difference perhaps)

That will at least give you the start and end points where convexity changes.  You could convert to points for further work.  

If you need detailed work, you will have to convert the poly* features to point pairs

I use numpy generally but this idea could be easily translated to cursors.

import numpy as np
def polyline_angles(bits, fromNorth=False):
    """Polyline/segment angles.

    Parameters
    ----------
    bits : array-like
        XY coordinate pairs, usually a numpy ndarray or list of pairs.
    fromNorth : boolean
        True for North calculated angles.  False for x-axis calculations.
    """
    out = []
    for b in bits:
        dxy = b[1:] - b[:-1]  # -- sequential difference in pair coordinates
        ang = np.degrees(np.arctan2(dxy[:, 1], dxy[:, 0]))
        if fromNorth:
            ang = np.mod((450.0 - ang), 360.)
        out.append(ang)
    return out

 


... sort of retired...
JohannesLindner
MVP Frequent Contributor

Starting from Dan's suggestion of convexHull(), which seemed good:

with arcpy.da.UpdateCursor("TestPolygons", ["SHAPE@", "IntegerField"]) as cursor:
    for shp, n in cursor:
        n = len(shp[0]) - len(shp.convexHull()[0])  # only for singlepart polygons
        cursor.updateRow([shp, n])

 

This works well for simple geometries:

JohannesLindner_0-1670325926043.png

 

But it starts to break down for more complex geometries, where convex vertices are not always on the convex hull:

JohannesLindner_1-1670326259996.png

 

I was ready to start fussing with angles between polylines when I realized that you can quite simply say whether a vertex is concave or convex: cretae a triangle of the vertex and its neighbours. If that triangle is inside the polygon, the vertex is convex, else it's concave.

JohannesLindner_2-1670329512390.png

 

def get_concave_vertex_count(polygon):
    n = 0
    for part in polygon:
        for i, vertex in enumerate(part):
            prev_vertex = part[i-1]
            try:
                next_vertex = part[i+1]
            except IndexError:
                next_vertex = part[0]
            triangle = arcpy.Polygon(arcpy.Array([prev_vertex, vertex, next_vertex]))
            n += triangle.touches(polygon)  # touch: only boundaries intersect
    return n


with arcpy.da.UpdateCursor("PolygonLayer", ["SHAPE@", "ConcaveVertices"]) as cursor:
    for shp, n in cursor:
        n = get_concave_vertex_count(shp)
        cursor.updateRow([shp, n])

Have a great day!
Johannes