GIS编程：利用Arcpy实现道格拉斯-普克算法（核心代码已用类包装，复制粘贴即可用）

        道格拉斯-普克算法是一种通过递归实现折线&多边形抽稀简化的算法，原理很简单，大多数文章都有讲解，这里不多赘述，直接上代码。

# -*- coding:utf-8 -*-
import arcpy
import sys
from arcpy import env
import numpy as np
import math

reload(sys)
sys.setdefaultencoding('utf-8')
env.overwriteOutput = True

# 读取一个shapefile, 用户自定义
# input a shapefile customized by user
shp = r'C:\Users\Lenovo\Desktop\11111\polyline.shp'


# 从.shp文件中获取多段线顶点坐标列表
# access the nodes' coordinate of each single feature
def getPolyXY(shape):
    with arcpy.da.SearchCursor(shape, ["SHAPE@"]) as cursor:
        Plines = []
        for row in cursor:
            coordXY = row[0].getPart()
            record = []
            for v in range(row[0].pointCount - 1):
                pnt = coordXY.getObject(0).getObject(v)
                record.append([pnt.X, pnt.Y])
            Plines.append(record)
        return Plines


class Douglas(object):
    def __init__(self):
        self.pts_seq = list()
        self.flag = list()
        self.sigma = 20
        self.head = 0
        self.rear = 0

    def update(self, pts_seq):
        self.pts_seq = pts_seq
        length = len(pts_seq)
        self.rear = length - 1
        self.flag = [False] * length

    def max_dis(self, head, rear):
        x0, y0 = self.pts_seq[head][0], self.pts_seq[head][1]
        x1, y1 = self.pts_seq[rear][0], self.pts_seq[rear][1]
        a = y1 - y0
        b = x0 - x1
        c = x0 * (y0 - y1) + y0 * (x1 - x0)
        max_d = 0
        index = 0
        for i in range(head + 1, rear):
            pt = self.pts_seq[i]
            di = np.abs(a * pt[0] + b * pt[1] + c) / np.sqrt(a ** 2 + b ** 2)
            if di > max_d:
                max_d = di
                index = i
        return max_d, index

    def compress(self, head, rear):
        if rear - head > 1:
            max_d, index = self.max_dis(head, rear)
            if max_d <= self.sigma:
                self.flag[head] = True
                self.flag[rear] = True
            else:
                self.compress(head, index)
                self.compress(index, rear)
        else:
            self.flag[head] = True
            self.flag[rear] = True

    def draw_polyline(self):
        polyline = []
        for i, bool_value in enumerate(self.flag):
            if bool_value:
                polyline.append(self.pts_seq[i])
        return polyline


def main():
    Plines = getPolyXY(shp)
    dgls = Douglas()
    dgls.sigma = input('请输入容差/plz input the threshold:')
    pline_out = []
    for pl in Plines:
        dgls.update(pl)
        dgls.compress(dgls.head, dgls.rear)
        pline_out.append(dgls.draw_polyline())

    # customize output file path&name
    out_path = r'C:\Users\Lenovo\Desktop\11111'
    out_name = 'polyline_cps.shp'
    fc = out_path + '\\' + out_name
    arcpy.CreateFeatureclass_management(out_path,
                                        out_name,
                                        geometry_type='POLYLINE')
    cursor = arcpy.InsertCursor(fc)
    for pl in pline_out:
        array = arcpy.Array()
        point = arcpy.Point()
        for pt in pl:
            point.X = pt[0]
            point.Y = pt[1]
            array.add(point)
        new_pl = arcpy.Polyline(array)
        row = cursor.newRow()
        row.shape = new_pl
        cursor.insertRow(row)


if __name__ == '__main__':
    main()

如果要运行Arcpy第三方库，在安装有arcGIS的windows电脑中，可直接通过程序内python加载项运行代码（如下图1所示）。如果希望在IDE环境中调试，要注意将编译器定位到GIS文件夹中的python.exe文件（如下图2）。

图1

图2

运行结果如下图所示，所设阈值为25.