便携设备可以方便的获取经纬度信息，如果按照一定的时间间隔就能获取到具体的行动轨迹。比如：
- 网约车平台通过获取实时的司机位置判定司机是否偏航
- 运动记录平台用户的运动轨迹，比如跑步路径等。
如果记录时间间隔较大，则会丢失很多数据。如果记录时间间隔较小，则会记录很多距离较近的点，比如网约车堵车。所谓的经纬度轨迹的抽稀，其实就是删除哪些不重要的点。以下是整理的一些方案。
为了更好的讲解数据，我们这里采用高德的线路规划接口的数据进行模拟。（苏州站→苏州园区站）
# 获取线路规划
import requests
def get_route(start_ , end_coordinate):
api_url = 'https://restapi.amap.com/v5/direction/driving'
payloads = {
"key": "",
"origin": start_coordinate, # 苏州站
"destination": end_coordinate, # 上海虹桥机场
"show_fields": "polyline"
}
r = requests.get(api_url, payloads)
data_json = r.json()
steps = data_json['route']['paths'][0]['steps']
return steps
steps = get_route(start_coordinate="120.610868,31.329679", end_coordinate="120.710574,31.341096")
print(steps)
线路规划返回的数据：
[{
'instruction': '向东行驶309米左转',
'orientation': '东',
'step_distance': '309',
'polyline': '120.610587,31.330743;120.611289,31.330878;120.611353,31.330904;120.611504,31.33099;120.611643,31.331012;120.611761,31.331033;120.612394,31.331146;120.612614,31.331194;120.612764,31.331226;120.61329,31.331307;120.613585,31.33136;120.613666,31.331398;120.613719,31.331451'
}, {
'instruction': '向北行驶118米右转',
'orientation': '北',
'step_distance': '118',
'polyline': '120.613719,31.331451;120.613735,31.331505;120.61373,31.331586;120.613585,31.332111;120.613499,31.332498'
}, {
'instruction': '沿苏站路向东行驶415米右转',
'orientation': '东',
'road_name': '苏站路',
'step_distance': '415',
'polyline': '120.613499,31.332498;120.61417,31.33261;120.614299,31.332632;120.615151,31.332782;120.615312,31.332809;120.616552,31.332996;120.616621,31.333007;120.617791,31.333195'
}, {
'instruction': '沿江乾路向南行驶490米左转',
'orientation': '南',
'road_name': '江乾路',
'step_distance': '490',
'polyline': '120.617791,31.333195;120.617871,31.333007;120.617936,31.332857;120.618343,31.332015;120.618483,31.33165;120.618633,31.331189;120.618837,31.330411;120.619196,31.328941'
}, {
'instruction': '沿西汇路途径东汇路向东行驶861米左转',
'orientation': '东',
'road_name': '西汇路',
'step_distance': '861',
'polyline': '120.619282,31.32885;120.62161,31.329258;120.622672,31.329456;120.62279,31.329461;120.623638,31.329655;120.624716,31.329874;120.625017,31.329906;120.625312,31.329912;120.626025,31.32988;120.626208,31.32988;120.626508,31.329912;120.627093,31.330008;120.627372,31.330073;120.628176,31.330261'
}, {
'instruction': '沿东汇路支路向北行驶165米右转',
'orientation': '北',
'road_name': '东汇路支路',
'step_distance': '165',
'polyline': '120.62815,31.330368;120.628037,31.330808;120.627967,31.331081;120.627946,31.331173;120.627844,31.33157;120.627822,31.331666;120.627779,31.331854'
}, {
'instruction': '沿北环东路向东行驶39米直行进入中间岔道',
'orientation': '东',
'road_name': '北环东路',
'step_distance': '39',
'polyline': '120.627779,31.331854;120.628182,31.331913'
}, {
'instruction': '沿北环快速路入口途径北环快速路向东行驶890米靠左沿主路行驶',
'orientation': '东',
'road_name': '北环快速路入口',
'step_distance': '890',
'polyline': '120.628182,31.331913;120.628568,31.331961;120.629185,31.332079;120.629818,31.332192;120.629882,31.332203;120.63265,31.332696;120.632795,31.33276;120.632935,31.332852;120.633428,31.332905;120.636314,31.333007;120.637392,31.333045'
}, {
'instruction': '沿北环快速路向东行驶1.1千米靠左沿主路行驶',
'orientation': '东',
'road_name': '北环快速路',
'step_distance': '1065',
'polyline': '120.637392,31.333045;120.640402,31.333173;120.64096,31.333206;120.641507,31.333243;120.644296,31.333426;120.645771,31.333495;120.646099,31.333517;120.646818,31.333522;120.647542,31.333565;120.648588,31.333667'
}, {
'instruction': '沿娄江快速路向东南行驶843米靠左沿主路行驶',
'orientation': '东南',
'road_name': '娄江快速路',
'step_distance': '843',
'polyline': '120.648588,31.333667;120.649312,31.333748;120.64972,31.33378;120.650305,31.333796;120.650653,31.333769;120.650964,31.33371;120.651227,31.33364;120.651597,31.33349;120.651855,31.333372;120.651946,31.333329;120.652343,31.333055;120.652407,31.333007;120.652713,31.332771;120.653019,31.332492;120.653523,31.331993;120.653936,31.331596;120.654323,31.33121;120.654467,31.33106;120.654709,31.330835;120.655315,31.330293;120.655685,31.329982'
}, {
'instruction': '沿娄江快速路向东行驶1.6千米靠左沿主路行驶',
'orientation': '东',
'road_name': '娄江快速路',
'step_distance': '1612',
'polyline': '120.655685,31.329982;120.656034,31.329756;120.656265,31.329628;120.656747,31.329391;120.657048,31.329263;120.65745,31.329123;120.657965,31.328984;120.658442,31.328887;120.658909,31.328828;120.659354,31.328796;120.659832,31.328791;120.660304,31.328812;120.660749,31.328866;120.661623,31.329011;120.665395,31.329762;120.665893,31.329874;120.666923,31.330084;120.670357,31.330776;120.671709,31.331044;120.67187,31.331076;120.672036,31.331114'
}, {
'instruction': '沿娄江快速路向东行驶2.3千米向右前方行驶进入匝道',
'orientation': '东',
'road_name': '娄江快速路',
'step_distance': '2285',
'polyline': '120.672036,31.331114;120.672299,31.331173;120.672567,31.331231;120.673034,31.331317;120.673828,31.331473;120.673844,31.331478;120.677111,31.332133;120.679583,31.332626;120.681332,31.332991;120.681402,31.333007;120.682995,31.333372;120.684996,31.333817;120.68534,31.333892;120.688006,31.334498;120.689733,31.334858;120.691798,31.335239;120.695425,31.335904'
}, {
'instruction': '沿玲珑街立交向东行驶201米靠左',
'orientation': '东',
'road_name': '玲珑街立交',
'step_distance': '201',
'polyline': '120.695425,31.335904;120.69572,31.335893;120.6958,31.335904;120.696562,31.336022;120.697157,31.336108;120.697511,31.336183'
}, {
'instruction': '沿玲珑街立交途径玲珑街向东行驶1.5千米到达目的地',
'orientation': '东',
'road_name': '玲珑街立交',
'step_distance': '1548',
'polyline': '120.697511,31.336183;120.697962,31.336258;120.69896,31.33644;120.700961,31.336821;120.701138,31.336864;120.701288,31.336923;120.701411,31.336993;120.701508,31.337057;120.701588,31.337138;120.701658,31.337229;120.701728,31.337347;120.701771,31.337454;120.701798,31.337578;120.701808,31.337712;120.701792,31.337846;120.701744,31.33806;120.701637,31.338463;120.701588,31.338838;120.701605,31.339021;120.701621,31.339128;120.701674,31.339294;120.701696,31.339348;120.701808,31.339541;120.70205,31.339783;120.70228,31.339938;120.702586,31.340078;120.703418,31.340308;120.703509,31.34034;120.703579,31.340383;120.70499,31.340705;120.706406,31.341027;120.70743,31.341247;120.70772,31.341279;120.70794,31.34129;120.708766,31.34122;120.709034,31.341231;120.709576,31.341311;120.710016,31.341403;120.710381,31.341472;120.710465,31.341489'
}]
线路规划返回的经纬度坐标点拼接后在地图上呈现如下：
points = []
for step in steps:
points.extend(step['polyline'].split(';'))
point_list = [(float(point.split(',')[1]), float(point.split(',')[0])) for point in points]
m = folium.Map(location=[31.363364, 120.638043], zoom_start=12, width=800, height=400, zoom_control='False',
tiles='https://map.geoq.cn/arcgis/rest/services/ChinaOnlineStreetGray/MapServer/tile/{z}/{y}/{x}',
attr='GeoQ')
folium.PolyLine(locations = point_list, color='green', weight=6, opacity=0.5).add_to(m)
for point in point_list:
folium.Circle(location=point,weight=1,color='green', fill=True,fill_opacity=0.8,opacity=1).add_to(m)
可以看到点之间非常的密集，有的甚至连成了直线。接下来要做的就是对这些点进行抽稀。
基于道路的抽稀
从上面的数据看，高德返回的数据是按照道路进行返回的，为了能得到抽稀的效果。可以采用选取每条道路起始点+最后道路的终点。代码示例：
points = []
for i, step in enumerate(steps):
if i == len(steps) - 1:
points.append(step['polyline'].split(';')[0])
points.append(step['polyline'].split(';')[-1])
else:
points.append(step['polyline'].split(';')[0])
m = folium.Map(location=[31.363364, 120.638043], zoom_start=12, width=800, height=400, zoom_control='False',
tiles='https://map.geoq.cn/arcgis/rest/services/ChinaOnlineStreetGray/MapServer/tile/{z}/{y}/{x}',
attr='GeoQ')
point_list = [(float(point.split(',')[1]), float(point.split(',')[0])) for point in points]
folium.PolyLine(locations = point_list, color='green', weight=6, opacity=0.5).add_to(m)
for point in point_list:
folium.Circle(location=point,weight=1,color='green', fill=True,fill_opacity=0.8,opacity=1).add_to(m)
由于按道路抽稀，线段的长短与行驶的道路长短相关，市内相较于其他区域更加密集，若涉及到高速可能会更长。
基于距离的抽稀
基于距离的抽稀比较好理解，比如每个500米选取一个点。大致的逻辑如下：
points = []
for step in steps:
points.extend(step['polyline'].split(';'))
point_list = [(float(point.split(',')[1]), float(point.split(',')[0])) for point in points]
import math
def haversine(latlon1, latlon2):
lat1, lon1 = latlon1
lat2, lon2 = latlon2
lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
c = 2 * math.asin(math.sqrt(a))
r = 6371 # 地球平均半径，单位为公里
return c * r
def distance_simplify(coords, epsilon):
result_list = [coords[0]]
for i in range(1, len(coords)):
a = result_list[-1]
b = coords[i]
if haversine(a, b) < epsilon:
pass
else:
result_list.append(coords[i])
return result_list
points_simplify = distance_simplify(point_list, 0.5)
m = folium.Map(location=[31.363364, 120.638043], zoom_start=12, width=800, height=400, zoom_control='False',
tiles='https://map.geoq.cn/arcgis/rest/services/ChinaOnlineStreetGray/MapServer/tile/{z}/{y}/{x}',
attr='GeoQ')
folium.PolyLine(locations = points_simplify, color='green', weight=6, opacity=0.5).add_to(m)
for point in points_simplify:
folium.Circle(location=point,weight=1,color='green', fill=True,fill_opacity=0.8,opacity=1).add_to(m)
抽稀后整体点与点之间的距离差不多。
基于Douglas-Peucker 算法的抽稀
道格拉斯-普克算法(Douglas–Peucker algorithm，亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点，并减少点的数量的一种算法。该算法的原始类型分别由乌尔斯·拉默（Urs Ramer）于1972年以及大卫·道格拉斯（David Douglas）和托马斯·普克（Thomas Peucker）于1973年提出，并在之后的数十年中由其他学者予以完善。
算法逻辑：
- 在轨迹曲线在曲线首尾两点 A，B 之间连接一条直线 AB，该直线为曲线的弦；
- 遍历曲线上其他所有点，求每个点到直线AB的距离，找到最大距离的点 C，最大距离记为 maxDistance
- 比较该距离 maxDistance 与预先定义的阈值 epsilon 大小，如果 maxDistance < epsilon，则将该直线 AB 作为曲线段的近似，舍去 AB 之间的所有点，曲线段处理完毕；
- 若 maxDistance >= epsilon，则使 C 点将曲线 AB 分为 AC和 CB 两段，并分别对这两段进行（1）~（3）步处理；
- 当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即为原始曲线的路径。
使用Python实现道格拉斯-普克算法：
import numpy as np
# 定义一个函数计算点到线的距离
def perpendicular_distance(point, line_start, line_end):
if np.all(np.equal(line_start, line_end)):
return np.linalg.norm(point - line_start)
return np.divide(
np.abs(np.linalg.norm(np.cross(line_end - line_start, line_start - point))),
np.linalg.norm(line_end - line_start)
)
# 实现道格拉斯-普克算法
def douglas_peucker(points, tolerance):
if len(points) < 3:
return points
# 初始化最远点的距离和索引
max_distance = 0
index_of_farthest = 0
# 查找最远点
for i in range(1, len(points) - 1):
distance = perpendicular_distance(np.array(points[i]), np.array(points[0]), np.array(points[-1]))
if distance > max_distance:
max_distance = distance
index_of_farthest = i
# 如果最大距离超过容忍度，递归地抽稀
if max_distance > tolerance:
# 递归处理线段两边的点
left_side = douglas_peucker(points[:index_of_farthest + 1], tolerance)
right_side = douglas_peucker(points[index_of_farthest:], tolerance)
# 合并结果
return left_side[:-1] + right_side
else:
return [points[0], points[-1]]
# 设置容忍度
tolerance = 0.0001
# 调用道格拉斯-普克算法
reduced_points = douglas_peucker(point_list, tolerance)
m = folium.Map(location=[31.363364, 120.638043], zoom_start=12, width=800, height=400, zoom_control='False',
tiles='https://map.geoq.cn/arcgis/rest/services/ChinaOnlineStreetGray/MapServer/tile/{z}/{y}/{x}',
attr='GeoQ')
folium.PolyLine(locations = reduced_points, color='green', weight=6, opacity=0.5).add_to(m)
for point in reduced_points:
folium.Circle(location=point,weight=1,color='green', fill=True,fill_opacity=0.8,opacity=1).add_to(m)
可以看到在拐弯处抽稀后比较密集，直线处比较稀疏。
具体在经纬度抽稀时使用何种方法，还需要看具体的使用场景，如果是界面上线程线路，建议使用格拉斯-普克算法。如果是做轨迹的相似度检验，建议使用按距离抽稀。
