文章記錄軌跡數(shù)據(jù)處理過程，其中使用到python

gpx數(shù)據(jù)來源：

鏈接: https://pan.baidu.com/s/1IELrjwXK1UYIMyQRnnLblQ
提取碼: n3en

import gpxpy  # 讀出軌跡數(shù)據(jù)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文檔

獲取軌跡數(shù)據(jù)，整理數(shù)據(jù)

with open('../gpx/test-1.gpx') as fh:
    gpx_file = gpxpy.parse(fh)
# 轉換為Pandas數(shù)據(jù)結構
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"總數(shù)據(jù)量：{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

原始數(shù)據(jù)展示

# 在地圖上顯示數(shù)據(jù)
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數(shù)據(jù)
# 速度值圖表
plt.plot(coords['speed'])
plt.title("speed-data")

image.png

刪除速度異常的數(shù)據(jù)

理論基礎：

箱形圖可以用來觀察數(shù)據(jù)整體的分布情況，利用中位數(shù)，25/%分位數(shù)，75/%分位數(shù)，上邊界，下邊界等統(tǒng)計量來來描述數(shù)據(jù)的整體分布情況。通過計算這些統(tǒng)計量，生成一個箱體圖，箱體包含了大部分的正常數(shù)據(jù)，而在箱體上邊界和下邊界之外的，就是異常數(shù)據(jù)。其中上下邊界的計算公式如下：

1. UpperLimit=Q3+1.5IQR=75%分位數(shù)+（75%分位數(shù)-25%分位數(shù)）*K
2. LowerLimit=Q1－1.5IQR=25%分位數(shù)-（75%分位數(shù)-25%分位數(shù)）*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"刪除異常數(shù)據(jù)點，{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
'刪除異常數(shù)據(jù)點，56'
# ------------------------------------
# 在地圖上顯示數(shù)據(jù)
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

壓縮數(shù)據(jù)

參考：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）二人又獨立于拉默提出了該算法。我們知道，一條曲線上包含著無數(shù)個點，但是計算機在存儲曲線時只能存取有限個點，通常存儲的點越多，那么對該曲線的描述也就越精確。當我們要對原本用N個點描述的曲線進行壓縮表示時，即采用K（K<N）個點來描述曲線，為了盡可能保證原有曲線的形態(tài)不至有太大的改變，我們就需要一種算法，而道格拉斯-普克算法就是這樣一種將曲線近似表示為一系列點，并減少點的數(shù)量的一種算法。它的優(yōu)點是具有平移和旋轉不變性，給定曲線與閾值后，抽樣結果一定。

# 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("優(yōu)化前：{}  優(yōu)化后：{} points".format(cut_plot_coords.shape[0], simplified_coords.shape[0]))

Kalman濾波器--平滑數(shù)據(jù)

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

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

import pykalman
import time
import numpy as np
# 平滑軌跡數(shù)據(jù)
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) # 次數(shù)不能過多
(smoothed_state_means, smoothed_state_covariances) = kf1.smooth(measurements)

# 濾波后經(jīng)度的對比
plt.figure(figsize=(20,20))
times = range(measurements.shape[0])
plt.plot(times, measurements[:, 0], 'b-',           # 優(yōu)化前
         #times, measurements[:, 1], 'ro',
         times, smoothed_state_means[:, 0], 'r--')  # 優(yōu)化后
         #times, smoothed_state_means[:, 2], 'r--',)
plt.show()

image.png

數(shù)據(jù)處理前后的對比

# 結果對比
fig = plt.figure(figsize=(20,20))
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values, 'r-')  # 藍線：未優(yōu)化路線
plt.plot(cut_plot_coords['lon'].values, cut_plot_coords['lat'].values, 'go')  # 紅點：未優(yōu)化定位點
plt.plot(simplified_coords[:,0], simplified_coords[:,1], 'y*')          #  黃色星星：壓縮后數(shù)據(jù)定位點
plt.plot(smoothed_state_means[:, 0], smoothed_state_means[:, 2], 'bx')  # 藍色x：平滑后數(shù)據(jù)的位置點

image.png