每一个GISer都应该了解的算法

前言：作为一个 GISer 需要掌握多方面的知识，哪怕不是开发方向，了解了解天天和我们打交道的 GIS 数据变化背后的原因也是很好的！除了社交、娱乐，人们外出使用最多的软件可能就是各种地图导航软件了，最常用的比如高德和百度地图。打开软件，双指滑动手机屏幕，放大缩小地图的时候，你可能会注意到这样的情况：放大地图，道路变得崎岖，弯来绕去；缩小地图，道路一下子变得笔直。就像下面这种情况，造成这种情况的原因就是 RDP 算法或者是相关的改进算法。同样在 GIS 软件中的抽稀和概化功能就是相应算法的实现。

广州市道路 截取自百度地图；左：道路较为笔直；右：道路较为曲折什么是 RDP 算法

RDP 算法的全称是拉默-道格拉斯-普克算法（英文：Ramer–Douglas–Peucker algorithm），分别是三位科学家的名字，又称道格拉斯-普克算法（DP），也称迭代端点拟合算法（英语：iterative end-point fit algorithm）。这是一种将曲线（折线）降采样为点数较少的类似曲线（折线）的算法，简单来说就是简化；是线状要素抽稀的经典算法（也可以叫概化），在保证几何形状基本不变的情况下，去除大量冗余折点，缩小体积。这对于网络地图是非常重要的，能减少加载时间，增强程序稳定性，保证终端设备的流畅，进而提升用户的体验。RDP 算法的应用非常广泛，特别是 GIS 领域。与 RDP 相似的另一个著名的算法是Visvalingam-Whyatt。大家可登陆 mapshaper.org 网站在线使用这两种算法，可调节参数，可导出处理后的 shp  数据。算法基本原理

将待处理曲线（折线）的首末端点连成一条直线，求所有中间点与直线的距离，并找出最大距离 max，用 max 与抽稀阈值 ε 相比较：若 max <= ε，这条曲线上的中间点全部舍去;若 max > ε，保留 max 对应的坐标点，并以该点为界，把曲线分为两部分，对这两部分曲线重复上述过程，直至所有的点都被处理完成；最后将首尾端点和保存下来的点相连，获得简化后的曲线（折现）。

原理演示 from wikipedia

原理演示（动态） from wikipedia没人看的代码实现

Python 已经有实现其算法的第三方库 rdp，你可以通过以下命令安装。pip install rdp该模块非常简洁，仅专注于实现 rdp 算法，其 Python 源代码如下："""rdp~~~Python implementation of the Ramer-Douglas-Peucker algorithm.:copyright: 2014-2016 Fabian Hirschmann <[e-mail]>:license: MIT, see LICENSE.txt for more details."""from math import sqrtfrom functools import partialimport numpy as npimport sysif sys.version_info[0] >= 3:xrange = rangedef pldist(point, start, end):"""Calculates the distance from ``point`` to the line givenby the points ``start`` and ``end``.:param point: a point:type point: numpy array:param start: a point of the line:type start: numpy array:param end: another point of the line:type end: numpy array"""if np.all(np.equal(start, end)):return np.linalg.norm(point - start)return np.divide(np.abs(np.linalg.norm(np.cross(end - start, start - point))),np.linalg.norm(end - start))def rdp_rec(M, epsilon, dist=pldist):"""Simplifies a given array of points.Recursive version.:param M: an array:type M: numpy array:param epsilon: epsilon in the rdp algorithm:type epsilon: float:param dist: distance function:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`"""dmax = 0.0index = -1for i in xrange(1, M.shape[0]):d = dist(M[i], M[0], M[-1])if d > dmax:index = idmax = dif dmax > epsilon:r1 = rdp_rec(M[:index + 1], epsilon, dist)r2 = rdp_rec(M[index:], epsilon, dist)return np.vstack((r1[:-1], r2))else:return np.vstack((M[0], M[-1]))def _rdp_iter(M, start_index, last_index, epsilon, dist=pldist):stk = []stk.append([start_index, last_index])global_start_index = start_indexindices = np.ones(last_index - start_index + 1, dtype=bool)while stk:start_index, last_index = stk.pop()dmax = 0.0index = start_indexfor i in xrange(index + 1, last_index):if indices[i - global_start_index]:d = dist(M[i], M[start_index], M[last_index])if d > dmax:index = idmax = dif dmax > epsilon:stk.append([start_index, index])stk.append([index, last_index])else:for i in xrange(start_index + 1, last_index):indices[i - global_start_index] = Falsereturn indicesdef rdp_iter(M, epsilon, dist=pldist, return_mask=False):"""Simplifies a given array of points.Iterative version.:param M: an array:type M: numpy array:param epsilon: epsilon in the rdp algorithm:type epsilon: float:param dist: distance function:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`:param return_mask: return the mask of points to keep instead:type return_mask: bool"""mask = _rdp_iter(M, 0, len(M) - 1, epsilon, dist)if return_mask:return maskreturn M[mask]def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):"""Simplifies a given array of points using the Ramer-Douglas-Peuckeralgorithm.Example:# >>> from rdp import rdp# >>> rdp([[1, 1], [2, 2], [3, 3], [4, 4]])[[1, 1], [4, 4]]This is a convenience wrapper around both :func:`rdp.rdp_iter`and :func:`rdp.rdp_rec` that detects if the input is a numpy arrayin order to adapt the output accordingly. This means thatwhen it is called using a Python list as argument, a Pythonlist is returned, and in case of an invocation using a numpyarray, a NumPy array is returned.The parameter ``return_mask=True`` can be used in conjunctionwith ``algo="iter"`` to return only the mask of points to keep. Example:# >>> from rdp import rdp# >>> import numpy as np# >>> arr = np.array([1, 1, 2, 2, 3, 3, 4, 4]).reshape(4, 2)# >>> arrarray([[1, 1],[2, 2],[3, 3],[4, 4]])# >>> mask = rdp(arr, algo="iter", return_mask=True)# >>> maskarray([ True, False, False, True], dtype=bool)# >>> arr[mask]array([[1, 1],[4, 4]]):param M: a series of points:type M: numpy array with shape ``(n,d)`` where ``n`` is the number of points and ``d`` their dimension:param epsilon: epsilon in the rdp algorithm:type epsilon: float:param dist: distance function:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`:param algo: either ``iter`` for an iterative algorithm or ``rec`` for a recursive algorithm:type algo: string:param return_mask: return mask instead of simplified array:type return_mask: bool"""if algo == "iter":algo = partial(rdp_iter, return_mask=return_mask)elif algo == "rec":if return_mask:raise NotImplementedError("return_mask=True not supported with algo=\"rec\"")algo = rdp_recif "numpy" in str(type(M)):return algo(M, epsilon, dist)return algo(np.array(M), epsilon, dist).tolist()该源代码中，计算每个点到直线距离封装成了一个单独的函数：pldist 。然后根据迭代、或者递归这两种编程技巧（方式）分别实现了 rdp，供用户选择。默认选择是递归实现。同样阈值（epsilon）的设置非常重要，下面有三张在不同阈值情况下，进行曲线（折线）简化后的可视化对比：其中蓝色曲线是原始数据的曲线，红色是简化后的曲线。

阈值=0 原始点：77 简化后：50；当阈值为0时，仅会舍弃一条直线上的多余节点

阈值=1 原始点：77 简化后：16

阈值=4 原始点：77 简化后：4随着阈值从0到1再到4，可以直观的看到曲线（蓝色）从平滑慢慢变得菱角分明（红色），鱼与熊掌不可得兼，既要精简的数据量，又要尽可能的好看平滑，这就是算法开发人员的“白鲸”，梦想鱼与熊掌的最大兼得。随便一提，阈值为0可用于删减同一条直线上的多余的点。最后

新的一天从新的知识开始，今天你又学废了嘛？其实我还想写一个结合 OGR 库使用 RDP 算法处理 shp 数据的小案例，结果发现字可能挺多的，就以后有机会再写吧！荟GIS精粹，关注公众号：GIS荟欢迎交流，更多文章请使用搜索