废话不多说，直接开干！
最近在做一个车联网项目，有一个场景是车辆定时上报当前所在经纬度等位置信息上报给平台，平台通过web页面在高德地图上展示车辆行驶路径。
说明
道格拉斯-普克算法 (Douglas–Peucker algorithm，亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点，并减少点的数量的一种算法。它的优点是具有平移和旋转不变性，给定曲线与阈值后，抽样结果一定。—摘自百度百科
python 代码
安装模块Shapely
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple Shapely
# -*- coding:utf-8 -*-
"""
道格拉斯算法的实现
程序需要安装shapely模块
"""
import math
from shapely import wkt, geometry
import matplotlib.pyplot as plt
class Point:
"""点类"""
x = 0.0
y = 0.0
index = 0 # 点在线上的索引
def __init__(self, x, y, index):
self.x = x
self.y = y
self.index = index
class Douglas:
"""道格拉斯算法类"""
points = []
D = 1 # 容差
def readPoint(self):
"""生成点要素"""
g = wkt.loads("LINESTRING(1 4,2 3,4 2,6 6,7 7,8 6,9 5,10 10)")
coords = g.coords
for i in range(len(coords)):
self.points.append(Point(coords[i][0], coords[i][1], i))
def compress(self, p1, p2):
"""具体的抽稀算法"""
swichvalue = False
# 一般式直线方程系数 A*x+B*y+C=0,利用点斜式,分母可以省略约区
# A=(p1.y-p2.y)/math.sqrt(math.pow(p1.y-p2.y,2)+math.pow(p1.x-p2.x,2))
A = (p1.y - p2.y)
# B=(p2.x-p1.x)/math.sqrt(math.pow(p1.y-p2.y,2)+math.pow(p1.x-p2.x,2))
B = (p2.x - p1.x)
# C=(p1.x*p2.y-p2.x*p1.y)/math.sqrt(math.pow(p1.y-p2.y,2)+math.pow(p1.x-p2.x,2))
C = (p1.x * p2.y - p2.x * p1.y)
m = self.points.index(p1)
n = self.points.index(p2)
distance = []
middle = None
if (n == m + 1):
return
# 计算中间点到直线的距离
for i in range(m + 1, n):
d = abs(A * self.points[i].x + B * self.points[i].y + C) / math.sqrt(math.pow(A, 2) + math.pow(B, 2))
distance.append(d)
dmax = max(distance)
if dmax > self.D:
swichvalue = True
else:
swichvalue = False
if (not swichvalue):
for i in range(m + 1, n):
del self.points[i]
else:
for i in range(m + 1, n):
if (abs(A * self.points[i].x + B * self.points[i].y + C) / math.sqrt(
math.pow(A, 2) + math.pow(B, 2)) == dmax):
middle = self.points[i]
self.compress(p1, middle)
self.compress(middle, p2)
def printPoint(self):
"""打印数据点"""
for p in self.points:
print( "%d,%f,%f" % (p.index, p.x, p.y))
def main():
"""测试"""
d = Douglas()
d.readPoint()
# d.printPoint()
# 结果图形的绘制，抽稀之前绘制
fig = plt.figure()
a1 = fig.add_subplot(121)
dx = []
dy = []
for i in range(len(d.points)):
dx.append(d.points[i].x)
dy.append(d.points[i].y)
a1.plot(dx, dy, color='g', linestyle='-', marker='+')
d.compress(d.points[0], d.points[len(d.points) - 1]) #稀释后轨迹
# 抽稀之后绘制
dx1 = []
dy1 = []
a2 = fig.add_subplot(122)
for p in d.points:
print(p.x,p.y)
dx1.append(p.x)
dy1.append(p.y)
a2.plot(dx1, dy1, color='r', linestyle='-', marker='+')
plt.show()
if __name__ == '__main__':
main()
运行结果
php代码
/**
* 根据两个点求直线方程 ax+by+c=0
* @param $xy1 string 点1,例如"1,1"
* @param $xy2 string 点2,例如"2,2"
* @return array
*/
function getLineByPoint($xy1, $xy2)
{
$x1 = explode(",", $xy1)[0];//第一个点的横坐标
$y1 = explode(",", $xy1)[1];//第一个点的纵坐标
$x2 = explode(",", $xy2)[0];//第二个点的横坐标
$y2 = explode(",", $xy2)[1];//第二个点的横坐标
$a = $y2 - $y1;
$b = $x1 - $x2;
$c = ($y1 - $y2) * $x1 - $y1 * ($x1 - $x2);
return [$a, $b, $c];
}
/**
* 稀疏点
* @param $points array 参数为["1,2","2,3"]点集
* @param $max float 阈值,越大稀疏效果越好但是细节越差
* @return array
*/
function sparePoints($points, $max)
{
if (count($points) < 3) {
return $points;
}
// var_dump($points);die;
$xy1 = $points[0];//取第一个点
$xy1 = $xy1['glng'].','.$xy1['glat'];
$xy2 = end($points);//取最后一个点
$xy2 = $xy2['glng'].','.$xy2['glat'];
// var_dump(end($points)['glat']);die;
list($a, $b, $c) = getLineByPoint($xy1, $xy2);//获取直线方程的a,b,c值
$ret = [];//最后稀疏以后的点集
$dmax = 0;//记录点到直线的最大距离
$split = 0;//分割位置
for ($i = 1; $i < count($points) - 1; $i++) {
$d = getDistanceFromPointToLine($a, $b, $c, $points[$i]);
if ($d > $dmax) {
$split = $i;
$dmax = $d;
}
}
if ($dmax>$max) {
//如果存在点到首位点连成直线的距离大于max的,即需要再次划分
$child1 = sparePoints(array_slice($points, 0, $split + 1), $max);//按split分成左边一份，递归
$child2 = sparePoints(array_slice($points, $split), $max);//按split分成右边一份，递归
//因为child1的最后一个点和child2的第一个点,肯定是同一个（即为分割点）,合并的时候,需要注意一下
$ret = array_merge($ret, $child1, array_slice($child2, 1));
return $ret;
} else {
//如果不存在点到直线的距离大于阈值的,那么就直接是首尾点了
return [$points[0], end($points)];
}
}
$gps_list= [
"118.7727996753945,34.1246971539229",
"118.77485177448638,34.124702346094594",
"118.7831576932094,34.13298337245958",
"118.78533669945193,34.133335134266446",
];
$data = sparePoints($gps_list, "0.0001");//稀释万分之一
运行结果
