交通轨迹抽稀-Douglas-peucker算法

Douglas-peucker算法(道格拉斯-普克算法)是一种经典的轨迹抽稀算法。

算法原理

• 将轨迹的起点(假如为A)和终点(假如为C)连接成一条直线；

• 遍历轨迹中的所有点(起终点除外)，计算各点到起终点连成的直线的距离，并选出最大的那个距离及其所对应的点(假如为B)；

• 如果最大距离小于所设定的偏移量Lmax，则该轨迹简化为起终点相连的直线；

• 如果最大距离大于等于所设定的偏移量Lmax，则以最大距离对应的点B为分界，将轨迹分为A-B和B-C两段，重复上述步骤；

• 处理完毕之后，将各个分界点相连，即可获得抽稀后的轨迹。

动画演示：

引自维基百科

代码演示：

笔者手边并没有可用的gps轨迹，因此使用高德-路径规划API生成一些gps轨迹，代码如下，

import requests
import folium
import warnings

warnings.filterwarnings("ignore")

def get_trace(origin, destination, key):
    """生成轨迹
    :param origin: 起点 格式: 120.638043,31.363364
    :param destination: 终点 格式: 120.710574,31.341096
    :param key: 用户key
    """
    base_url = f'''https://restapi.amap.com/v3/direction/driving?origin={origin}&destination={destination}&extensions=all&output=json&key={key}'''
    response = requests.get(base_url).json()
    steps = response['route']['paths'][0]['steps']
    
    points = []
    for step in steps:
        points.extend(step['polyline'].split(';'))

    data = pd.DataFrame(points, columns=['location'])
    data['lon'] = data['location'].apply(lambda x: float(x.split(',')[0]))
    data['lat'] = data['location'].apply(lambda x: float(x.split(',')[1]))
    
    gdata = gpd.GeoDataFrame(data, geometry=gpd.points_from_xy(data['lon'], data['lat']))
    return gdata

gdf = get_trace(origin='121.489663,31.141301', destination='121.320081,31.193964', key='你的key')
gdf.to_file('./save_data/轨迹.geojson', driver='GeoJSON')

接下来可以使用folium库可视化一下轨迹，

def plot_trace():
    gdf = gpd.read_file('./save_data/轨迹.geojson')
    
    tiles_gd = 'http://wprd04.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}'
    map = folium.Map(location=[31.193964, 121.320081], zoom_start=12, tiles=tiles_gd, attr='default')
    
    # 创建标记点样式
    marker = folium.Circle(radius=5, color='black', fill=True, fill_color='black')
    folium.GeoJson(gdf, marker=marker).add_to(map)
    map.save('./save_data/轨迹.html')

plot_trace()

可以看出轨迹点之间较为密集。

下面进行Douglas-peucker算法演示

先计算点到直线的距离，这里使用欧式距离近似计算，感兴趣的同学可以尝试用Haversine公式来进行计算

def get_distance(point, start_point, end_point):
    """计算点到直线的距离"""
    if (start_point.x == end_point.x) & (start_point.y == end_point.y):
        # 如果起终点相同，则取任一个点为垂足
        return np.sqrt(np.power(point.x - start_point.x, 2) + np.power(point.y - start_point.y, 2))
    
    # 根据向量算点到直线的距离
    start_point = np.array([start_point.x, start_point.y])
    end_point = np.array([end_point.x, end_point.y])
    point = np.array([point.x, point.y])
    
    # np.dot用来求向量点积
    # np.linalg.norm用来求向量的模
    return (np.dot(np.array(point) - np.array(start_point), np.array(end_point) - np.array(start_point)) / 
            np.linalg.norm(np.array(end_point) - np.array(start_point)))

Douglas-peucker算法计算过程如下，

def base_douglas_peucker(points, width):
    """基于Douglas-Peucker算法的轨迹抽稀 根据算法思想, 可能看成一个递归问题"""
    # 给出递归终止条件
    if len(points) < 3:
        return points
    
    # 计算距离最远的点
    max_distance = 0       
    max_distance_index = 0   # 最大距离对应的点在points中的索引
    
    for i in range(1, len(points) - 2):
        distance = get_distance(points[i], points[0], points[-1])
        if distance > max_distance:
            max_distance = distance
            max_distance_index = i
    
    # 判断最大距离是否大于所设定的阈值width
    if max_distance > width:
        # 则递归调用, 分别计算左右轨迹
        left_trace = base_douglas_peucker(points[:max_distance_index + 1], width)
        right_trace = base_douglas_peucker(points[max_distance_index:], width)
        return left_trace + right_trace[1: ]
    else:
        return [points[0], points[-1]]

trace = base_douglas_peucker(points, 0.01)

# 转成GeoDataFrame
gdata = gpd.GeoDataFrame(a, columns=['geometry'])
gdata.to_file('轨迹_douglas_peucker.geojson', driver='GeoJSON')

# 可视化
def plot_trace1():
    gdf = gpd.read_file('轨迹_douglas_peucker.geojson')
    
    tiles_gd = 'http://wprd04.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}'
    map = folium.Map(location=[31.193964, 121.320081], zoom_start=12, tiles=tiles_gd, attr='default')
    
    # 创建标记点样式
    marker = folium.Circle(radius=5, color='black', fill=True, fill_color='black')
    folium.GeoJson(gdf, marker=marker).add_to(map)
    map.save('轨迹_douglas_peucker.html')
plot_trace1()