便携设备可以方便的获取经纬度信息,如果按照一定的时间间隔就能获取到具体的行动轨迹。比如:
- 网约车平台通过获取实时的司机位置判定司机是否偏航
- 运动记录平台用户的运动轨迹,比如跑步路径等。
如果记录时间间隔较大,则会丢失很多数据。如果记录时间间隔较小,则会记录很多距离较近的点,比如网约车堵车。所谓的经纬度轨迹的抽稀,其实就是删除哪些不重要的点。以下是整理的一些方案。
为了更好的讲解数据,我们这里采用高德的线路规划接口的数据进行模拟。(苏州站→苏州园区站)
# 获取线路规划 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)
可以看到在拐弯处抽稀后比较密集,直线处比较稀疏。
具体在经纬度抽稀时使用何种方法,还需要看具体的使用场景,如果是界面上线程线路,建议使用格拉斯-普克算法。如果是做轨迹的相似度检验,建议使用按距离抽稀。