Clean polygons boundaries

[DanPatterson]

Start with a polygons boundary. 

Where are the unnecessary points?

• they can be points that are duplicate
• segments that overlap and/or reverse directions (squish a Z to a _ ... hard to see the middle 2 points ehh?)
• and other mishaps

Can you spot the redundant points?  Original image on the left, cleaned image on the right?

Perhaps if I make it more obvious.  A densified version of the inputs and the cleaned version.

Better?

The code

def _clean_segments_(a, tol=1e-06):
    """Remove overlaps and extra points on poly* segments. In `npg_geom_hlp`.

    Parameters
    ----------
    a : array
        The input array or a bit from a Geo array.
    tol : float
        The tolerance for determining whether a point deviates from a line.

    Notes
    -----
    - Segments along a straight line can overlap (a construction error).
          [[0,0], [5, 5], [2, 2], [7, 7]]  # points out of order
    - Extraneous points can exist along a segment.
          [[0,0], [2, 2], [5, 5], [7, 7]]  # extra points not needed for line.
    """
    cr, ba, bc = _bit_crossproduct_(a, extras=True)
    # -- avoid duplicating the 1st point (0).
    whr = [i for i in np.nonzero(np.abs(cr) > tol)[0] if i != 0]
    vals = np.concatenate((a[0][None, :], a[whr], a[-1][None, :]), axis=0)
    return vals

Which calls

def _bit_crossproduct_(a, is_closed=True, extras=False):
    """Cross product.

    Used by `is_convex`, `_angles_3pnt_` and `_clean_segments_`.

    .. note::

        np.cross for 2D arrays was deprecated in numpy 2.0, use `cross_2d`
    """
    def cross2d(x, y):
        return x[..., 0] * y[..., 1] - x[..., 1] * y[..., 0]

    if is_closed:
        if np.allclose(a[0], a[-1]):  # closed loop, remove dupl.
            a = a[:-1]
    ba = a - np.concatenate((a[-1][None, :], a[:-1]), axis=0)
    bc = a - np.concatenate((a[1:], a[0][None, :]), axis=0)
    # cr = np.cross(ba, bc) + 0.0  # deprecated
    cr = cross2d(ba, bc)
    if extras:
        return cr, ba, bc
    return cr

Now vals is a NumPy array which can be converted back to an arcpy shape using _arr_poly_ and _poly_to_array_ can be used to get the array(s) from arcpy geometry.

import arcpy
from arcpy import Array, Exists, Multipoint, Point, Polygon, Polyline

def _arr_poly_(arr, SR, as_type):
    """Slice the array where nan values appear, splitting them off."""
    aa = [Point(*pairs) for pairs in arr]
    if as_type.upper() == 'POLYGON':
        poly = Polygon(Array(aa), SR)
    elif as_type.upper() == 'POLYLINE':
        poly = Polyline(Array(aa), SR)
    return poly

def _poly_to_array_(polys):
    """Convert polyline or polygon shapes to arrays for use in numpy.

    Parameters
    ----------
    polys : tuple, list
        Polyline or polygons in a list/tuple
    """
    def _p2p_(poly):
        """Convert a single ``poly`` shape to numpy arrays or object."""
        sub = []
        pt = Point()  # arcpy.Point()
        for arr in poly:
            pnts = [[p.X, p.Y] for p in arr if pt]
            sub.append(np.asarray(pnts, dtype='O'))
        return sub
    # ----
    if not isinstance(polys, (list, tuple)):
        polys = [polys]
    out = []
    for poly in polys:
        out.extend(_p2p_(poly))  # or append, extend it is
    return out

Enough for now.  More on my github page numpy geometry