pandas 如果speed标签本身没有给出,我如何从gpx文件计算速度?

64jmpszr  于 2023-01-28  发布在  其他
关注(0)|答案(1)|浏览(131)

例如,我可能有这样的一些信息:

<trkpt lat="-33.8161780" lon="150.8710320">
 <ele>73.0</ele>
 <time>2017-07-08T22:05:45Z</time>
 <extensions>
  <power>0</power>
  <gpxtpx:TrackPointExtension>
    <gpxtpx:atemp>7</gpxtpx:atemp>
    <gpxtpx:hr>115</gpxtpx:hr>
    <gpxtpx:cad>27</gpxtpx:cad>
  </gpxtpx:TrackPointExtension>
 </extensions>
</trkpt>

我该如何从这些信息中计算速度呢?我已经使用了python etree模块来解析这个文件,并且把所有的信息都放在了Pandas数据库中。
有人提到我可能应该显示Pandas的数据框,它看起来像这样:

longitude   latitude   ele   temp
time                
2017-07-08 22:05:45 150.8710320 -33.8161780 73.0    7
2017-07-08 22:05:46 150.8710350 -33.8161500 73.0    7
2017-07-08 22:05:47 150.8710440 -33.8161170 73.0    7
2017-07-08 22:05:48 150.8710540 -33.8160820 73.0    7
2017-07-08 22:05:49 150.8710690 -33.8160430 73.0    7

等等。

gmxoilav

gmxoilav1#

speed等于distance / timelongitudelatitude大概代表地球表面的位置。如果我们接受半径为6371 km的球体作为地球的近似,那么我们可以很容易地将longitudelatitude转换为xyz坐标:

r = 6371000 # meters
df['theta'] = np.deg2rad(df['longitude'])
df['phi'] = np.deg2rad(df['latitude'])
df['x'] = r*np.cos(df['theta'])*np.sin(df['phi'])
df['y'] = r*np.sin(df['theta'])*np.sin(df['phi'])
df['z'] = r*np.cos(df['phi'])

现在计算连续点之间的距离并不困难:

df['x2'] = df['x'].shift()
df['y2'] = df['y'].shift()
df['z2'] = df['z'].shift()
df['distance'] = np.sqrt((df['x2']-df['x'])**2 + (df['y2']-df['y'])**2 + (df['z2']-df['z'])**2)

然而,这是一个弦长--球体表面上两点之间的直线距离。如果两点相距很远,弦将穿过地球表面。假定运动是在地球表面上进行的。因此,更精确的距离计算将使用arclength

df['central angle'] = np.arccos((df['x']*df['x2'] + df['y']*df['y2'] + df['z']*df['z2'])/r**2)
df['arclength'] = df['central angle']*r

中心角使用dot product formula
计算完弧长(距离)后,我们还必须计算连续观测(即数据框的行)之间的time间隔:

df['time'] = (df.index.to_series().diff() / pd.Timedelta(seconds=1))

因此,使用speed = distance / time

df['speed'] = df['arclength'] / df['time']  # in meters/second
import numpy as np
import pandas as pd

df = pd.DataFrame({'ele': [73.0, 73.0, 73.0, 73.0, 73.0], 'latitude': [-33.816178, -33.81615, -33.816117, -33.816082, -33.816043], 'longitude': [150.871032, 150.871035, 150.87104399999998, 150.87105400000002, 150.871069], 'temp': [7, 7, 7, 7, 7], 'time': ['2017-07-08 22:05:45', '2017-07-08 22:05:46', '2017-07-08 22:05:47', '2017-07-08 22:05:48', '2017-07-08 22:05:49']})
df['time'] = pd.to_datetime(df['time'])
df = df.set_index('time')
columns = df.columns.tolist()

r = 6371000 # radius of the Earth in meters
df['theta'] = np.deg2rad(df['longitude'])
df['phi'] = np.deg2rad(df['latitude'])
df['x'] = r*np.cos(df['theta'])*np.sin(df['phi'])
df['y'] = r*np.sin(df['theta'])*np.sin(df['phi'])
df['z'] = r*np.cos(df['phi'])
df['x2'] = df['x'].shift()
df['y2'] = df['y'].shift()
df['z2'] = df['z'].shift()
df['distance'] = np.sqrt((df['x2']-df['x'])**2 + (df['y2']-df['y'])**2 + (df['z2']-df['z'])**2)

df['central angle'] = np.arccos((df['x']*df['x2'] + df['y']*df['y2'] + df['z']*df['z2'])/r**2)
df['arclength'] = df['central angle']*r

df['time'] = (df.index.to_series().diff() / pd.Timedelta(seconds=1))
df['speed'] = df['arclength'] / df['time']  # in meters/second
df = df[columns + ['speed']]
print(df)

收益率

ele   latitude   longitude  temp     speed
time                                                            
2017-07-08 22:05:45  73.0 -33.816178  150.871032     7       NaN
2017-07-08 22:05:46  73.0 -33.816150  150.871035     7  3.119892
2017-07-08 22:05:47  73.0 -33.816117  150.871044     7  3.712201
2017-07-08 22:05:48  73.0 -33.816082  150.871054     7  3.940673
2017-07-08 22:05:49  73.0 -33.816043  150.871069     7  4.433590

如果您注解掉

df = df[columns + ['speed']]

并重新运行脚本,您将看到所有中间计算。您将注意到df['distance']非常接近df['arclength']。由于地球表面上的点相距不远,因此弦长是弧长的一个很好的近似值。因此,对于您发布的数据

df['speed'] = df['distance'] / df['time']

本来也可以工作得一样好。但是,arclength的计算更健壮一点,因为如果点相距较远,它会给出更准确的值。

相关问题