轨迹数据处理

文章记录轨迹数据处理过程，其中使用到python

gpx数据来源：

链接: https://pan.baidu.com/s/1IELrjwXK1UYIMyQRnnLblQ
提取码: n3en

import gpxpy  # 读出轨迹数据GPX
import pandas as pd

import matplotlib.pyplot as plt
plt.rcParams['axes.xmargin'] = 0.1
plt.rcParams['axes.ymargin'] = 0.1
%matplotlib inline

import seaborn as sns
sns.set_style("whitegrid")
sns.set_context("talk")

%load_ext autoreload
%autoreload 2

读取GPX文档

获取轨迹数据，整理数据

with open('../gpx/test-1.gpx') as fh:
    gpx_file = gpxpy.parse(fh)
# 转换为Pandas数据结构
segment = gpx_file.tracks[0].segments[0]
coords = pd.DataFrame([{'lat': p.latitude, 
                        'lon': p.longitude, 
                        'time': p.time} for p in segment.points])
coords.set_index('time', drop=True, inplace=True)
print(f"总数据量：{len(coords)}")
coords.head()

time lat lon
2020-10-30 22:54:03+00:00 30.602108 119.897750
2020-10-30 22:54:03+00:00 30.602107 119.897749
2020-10-30 22:55:27+00:00 30.602098 119.897479
2020-10-30 22:56:27+00:00 30.602219 119.897779
2020-10-30 22:58:05+00:00 30.602341 119.898196

原始数据展示

# 在地图上显示数据
plot_coords = coords[['lat','lon']].reset_index()
plot_coords.drop('time', axis=1, inplace=True)
fig = plt.figure()
plt.plot(plot_coords['lon'].values, plot_coords['lat'].values)
plt.plot(plot_coords['lon'].values, plot_coords['lat'].values, 'ro')

image.png

删除异常点

# 获取开始和结束时间
start_time, end_time = segment.get_time_bounds()
duration = end_time - start_time
avg_speed = segment.length_2d() / duration.seconds
print(f"平均速度：{avg_speed:.2f} m/s , {avg_speed * 3.6:.2f} km/h")
# 处理没有速度的异常点
segment.points[0].speed = 0.0
segment.points[-1].speed = 0.0
gpx_file.add_missing_speeds()
coords['speed'] = [p.speed for p in segment.points]
coords['speed'] *= 3.6
coords['speed'].fillna(avg_speed, inplace=True) # 填充NaN数据
# 速度值图表
plt.plot(coords['speed'])
plt.title("speed-data")

image.png

删除速度异常的数据

理论基础：

箱形图可以用来观察数据整体的分布情况，利用中位数，25/%分位数，75/%分位数，上边界，下边界等统计量来来描述数据的整体分布情况。通过计算这些统计量，生成一个箱体图，箱体包含了大部分的正常数据，而在箱体上边界和下边界之外的，就是异常数据。其中上下边界的计算公式如下：

1. UpperLimit=Q3+1.5IQR=75%分位数+（75%分位数-25%分位数）*K
2. LowerLimit=Q1－1.5IQR=25%分位数-（75%分位数-25%分位数）*K

k=1.5时，计算出的是中度异常的范围。K=3计算出的是极度异常的范围

参考：https://www.zhihu.com/question/36172806

descri_speed = coords['speed'].describe()
print(coords['speed'].describe())
cut_coords = coords[coords['speed'] < descri_speed['75%'] + 3*(descri_speed['75%'] - descri_speed['25%']) ]
print(f"删除异常数据点，{len(coords['speed'])-len(cut_coords)}")
# --------------------------------------
# 输出结果
count     449.000000
mean       72.734144
std       345.554273
min         0.000000
25%         2.082639
50%         3.535751
75%         7.593871
max      4994.907107
Name: speed, dtype: float64
'删除异常数据点，56'
# ------------------------------------
# 在地图上显示数据
cut_plot_coords = cut_coords[['lat','lon']].reset_index()
cut_plot_coords.drop('time', axis=1, inplace=True)
cut_fig = plt.figure()
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values)
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values, 'ro')

image.png

压缩数据

参考：https://baimafujinji.blog.csdn.net/article/details/6475432

Ramer–Douglas–Peucker algorithm

道格拉斯-普克算法（Douglas–Peucker algorithm），亦称为拉默-道格拉斯-普克算法（Ramer–Douglas–Peucker algorithm），这个算法最初由拉默（Urs Ramer）于1972年提出，1973年道格拉斯（David Douglas）和普克（Thomas Peucker）二人又独立于拉默提出了该算法。我们知道，一条曲线上包含着无数个点，但是计算机在存储曲线时只能存取有限个点，通常存储的点越多，那么对该曲线的描述也就越精确。当我们要对原本用N个点描述的曲线进行压缩表示时，即采用K（K<N）个点来描述曲线，为了尽可能保证原有曲线的形态不至有太大的改变，我们就需要一种算法，而道格拉斯-普克算法就是这样一种将曲线近似表示为一系列点，并减少点的数量的一种算法。它的优点是具有平移和旋转不变性，给定曲线与阈值后，抽样结果一定。

# RDP algorithm as long as the rdp package is not iterative.
# See https://github.com/fhirschmann/rdp/issues/5
simplified_coords = rdp(cut_plot_coords[['lon', 'lat']].values, 0.0001)
print("优化前：{}  优化后：{} points".format(cut_plot_coords.shape[0], simplified_coords.shape[0]))

Kalman滤波器--平滑数据

JAVA例子： https://stackoverflow.com/questions/1134579/smooth-gps-data

原理讲解： https://zhuanlan.zhihu.com/p/306093273

import pykalman
import time
import numpy as np
# 平滑轨迹数据
measurements = np.array(simplified_coords)

initial_state_mean = [measurements[0, 0],
                      0,
                      measurements[0, 1],
                      0]

transition_matrix = [[1, 1, 0, 0],
                     [0, 1, 0, 0],
                     [0, 0, 1, 1],
                     [0, 0, 0, 1]]

observation_matrix = [[1, 0, 0, 0],
                      [0, 0, 1, 0]]

kf1 = pykalman.KalmanFilter(transition_matrices = transition_matrix,
                  observation_matrices = observation_matrix,
                  initial_state_mean = initial_state_mean)

kf1 = kf1.em(measurements, n_iter=30) # 次数不能过多
(smoothed_state_means, smoothed_state_covariances) = kf1.smooth(measurements)

# 滤波后经度的对比
plt.figure(figsize=(20,20))
times = range(measurements.shape[0])
plt.plot(times, measurements[:, 0], 'b-',           # 优化前
         #times, measurements[:, 1], 'ro',
         times, smoothed_state_means[:, 0], 'r--')  # 优化后
         #times, smoothed_state_means[:, 2], 'r--',)
plt.show()

image.png

数据处理前后的对比

# 结果对比
fig = plt.figure(figsize=(20,20))
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values, 'r-')  # 蓝线：未优化路线
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values, 'go')  # 红点：未优化定位点
plt.plot(simplified_coords[:,0], simplified_coords[:,1], 'y*')          #  黄色星星：压缩后数据定位点
plt.plot(smoothed_state_means[:, 0], smoothed_state_means[:, 2], 'bx')  # 蓝色x：平滑后数据的位置点

image.png