1

I have a shapefile with polygons and want to write the Python script to eliminate collinear vertices in a polygon. I found a solution for the line feature class, but not for polygon: Deleting collinear vertex / vertices using ArcMap.

import arcpy
    
#----------------------------------------------------------------------
def are_collinear(p1, p2, p3, tolerance=0.5):
    """return True if 3 points are collinear.
    tolerance value will decide whether lines are collinear; may need
    to adjust it based on the XY tolerance value used for feature class"""
    x1, y1 = p1[0], p1[1]
    x2, y2 = p2[0], p2[1]
    x3, y3 = p3[0], p3[1]
    res = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    if -tolerance <= res <= tolerance:
        return True


#----------------------------------------------------------------------
def get_redundant_vertices(vertices):
    """get redundant vertices from a line shape vertices"""
    indexes_of_vertices_to_remove = []
    start_idx, middle_index, end_index = 0, 1, 2
    for i in range(len(vertices)):
        p1, p2, p3 = vertices[start_idx:end_index + 1]
        if are_collinear(p1, p2, p3):
            indexes_of_vertices_to_remove.append(middle_index)

        start_idx += 1
        middle_index += 1
        end_index += 1
        if end_index == len(vertices):
            break
    return indexes_of_vertices_to_remove


#----------------------------------------------------------------------
def clean_geometries(fc, densify_curves=False):
    """clean polyline features in the fc removing redundant vertices"""
    in_sr = arcpy.Describe(fc).spatialReference.factoryCode

    with arcpy.da.UpdateCursor(fc, ['SHAPE@', 'OID@']) as ucur:
        for row in ucur:
            print "OBJECTID", row[1]
            cleaned_parts = []
            shape = row[0]

            if 'curvePaths' in shape.JSON:
                if densify_curves:
                    shape = shape.densify('DISTANCE', 1, 1)
                else:
                    continue

            for part in range(shape.partCount):
                vertices = [(p.X, p.Y) for p in shape.getPart(part)]
                if len(vertices) < 3:  #polyline's part consists of 2 vertices
                    continue
                vertices_to_remove = get_redundant_vertices(vertices)
                vertices_to_keep = [
                    val for idx, val in enumerate(vertices)
                    if idx not in vertices_to_remove
                ]
                cleaned_part_as_array = arcpy.Array(
                    [arcpy.Point(*coords) for coords in vertices_to_keep])
                cleaned_parts.append(cleaned_part_as_array)

            if cleaned_parts:
                cleaned_shape = arcpy.Polyline(
                    arcpy.Array(cleaned_parts), in_sr)
                row[0] = cleaned_shape
                ucur.updateRow(row)


if __name__ == '__main__':
    fc = r'C:\GIS\Temp\ArcGISHomeFolder\Default.gdb\_SimpleRoads'
    #r'C:\GIS\Temp\ArcGISHomeFolder\empty.gdb\roads'
    clean_geometries(fc, densify_curves=False)

If I use this code for a polygon FC it will delete my polygons. I am new to Python and I don't know how to fix it in ArcPy.

Could you give me some tips or ideas how to achieve it? Unwanted vertices are marked with circles

2
  • Any Ramer-Douglas-Peucker implementation with a really tight threshold ought to accomplish this. For ArcGIS, with polygons, that could be arcpy.SimplifyPolygon_cartography.
    – Vince
    Feb 14, 2021 at 14:01
  • @Vince actually I tried RDP algorithm. If I use threshold like 0.5 it deletes some of the needed vertices, if I set 0.001 it deletes nothing. E.g.: s = df.simplify(tolerance=0.001, preserve_topology=True) Feb 14, 2021 at 14:10

1 Answer 1

2

I believe you only need to edit a few lines of code.

The lines I have edited/added are marked with ### at the end of the line. The key method is boundary(), it returns the polygon as a polyline and feeds that into the rest of the existing logic. The update cursor also now updates the geometry as a polygon.

So this code now expects the input to be a polygon dataset, not tested it with multi-parts either.

import arcpy

#----------------------------------------------------------------------
def are_collinear(p1, p2, p3, tolerance=0.5):
    """return True if 3 points are collinear.
    tolerance value will decide whether lines are collinear; may need
    to adjust it based on the XY tolerance value used for feature class"""
    x1, y1 = p1[0], p1[1]
    x2, y2 = p2[0], p2[1]
    x3, y3 = p3[0], p3[1]
    res = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    if -tolerance <= res <= tolerance:
        return True


#----------------------------------------------------------------------
def get_redundant_vertices(vertices):
    """get redundant vertices from a line shape vertices"""
    indexes_of_vertices_to_remove = []
    start_idx, middle_index, end_index = 0, 1, 2
    for i in range(len(vertices)):
        p1, p2, p3 = vertices[start_idx:end_index + 1]
        if are_collinear(p1, p2, p3):
            indexes_of_vertices_to_remove.append(middle_index)

        start_idx += 1
        middle_index += 1
        end_index += 1
        if end_index == len(vertices):
            break
    return indexes_of_vertices_to_remove


#----------------------------------------------------------------------
def clean_geometries(fc, densify_curves=False):
    """clean polyline features in the fc removing redundant vertices"""
    in_sr = arcpy.Describe(fc).spatialReference.factoryCode

    with arcpy.da.UpdateCursor(fc, ['SHAPE@', 'OID@']) as ucur:
        for row in ucur:
            print "OBJECTID", row[1]
            cleaned_parts = []
            poly = row[0] ###
            shape = poly.boundary() ###

            if 'curvePaths' in shape.JSON:
                if densify_curves:
                    shape = shape.densify('DISTANCE', 1, 1)
                else:
                    continue

            for part in range(shape.partCount):
                vertices = [(p.X, p.Y) for p in shape.getPart(part)]
                if len(vertices) < 3:  #polyline's part consists of 2 vertices
                    continue
                vertices_to_remove = get_redundant_vertices(vertices)
                vertices_to_keep = [
                    val for idx, val in enumerate(vertices)
                    if idx not in vertices_to_remove
                ]
                cleaned_part_as_array = arcpy.Array(
                    [arcpy.Point(*coords) for coords in vertices_to_keep])
                cleaned_parts.append(cleaned_part_as_array)

            if cleaned_parts:
                cleaned_shape = arcpy.Polygon(arcpy.Array(cleaned_parts), in_sr) ###
                row[0] = cleaned_shape
                ucur.updateRow(row)


if __name__ == '__main__':
    fc = r'C:\Scratch\fGDB_Scratch.gdb\testpoly_1'
    #r'C:\GIS\Temp\ArcGISHomeFolder\empty.gdb\roads'
    clean_geometries(fc, densify_curves=False)
    print("done!")
6
  • Thank you! I tried out and it deleted all my polygons :( Feb 15, 2021 at 12:28
  • I think you need to share some of your data so others can look at it? Maybe its to do with the tolerance value, there is an assumption that this data is NOT WGS84?
    – Hornbydd
    Feb 15, 2021 at 21:17
  • I guess, it is WGS84. Thank you so much! I attached the link in the bottom of my post. Feb 16, 2021 at 16:58
  • I think the issue is with the coordinate system, yours being in WGS84 and the tolerance setting for the function are_collinear. 0.5 would be 0.5 decimal degrees, not sure what that would be in metres. So you need to re-project your data to be using such code, Looking on SpatialReference website there is this coordinate system for Kazan?
    – Hornbydd
    Feb 16, 2021 at 18:47
  • Oh, right! Now I've got your point! I will try with the right one, thank you so much! :) Feb 16, 2021 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.