【大数据课程设计】出租车轨迹数据分析

在实验中，采用的数据是上海市某天24小时内的出租车GPS定位数据。

数据探索

原始给定的数据共有4316个文件，也即表示共有4316辆出租车数据，这里面不乏一些报废的出租车数据，这里会在数据清洗中给出，下面我们就一辆出租车的数据进行分析，拿编号位105的出租车数据来说，其数据格式如图所示：

上图的数据向我们反映出，每辆出租车的数据中都有很多条目的数据，分别记录了2007年2月20号从0点开始到这一天结束出租车所在位置以及载客情况的一些信息，通过分析我们可以知道，第一列属性表示出租车的编号，第二列属性较长，表示日期及时间，时间使精确到秒，第三和第四列属性分别是经度和纬度信息，这对我们而言已经不陌生了，第四第五列属性对我们而言并没有太大价值，因为我们在研究速度的时候会自己定义速度，最后一个属性则是载客情况，其实并不只是0和1，还有非零的其他数字也即表示载客的数量，所以这里在后面判断载客情况的时候需要有所注意。

数据预处理

首先要对数据进行的处理便是删除含有缺失值的数据。对所有含有缺失值的行都进行删除，这样可避免其对其他的数据产生影响
其次，删除所有载客距离小于300的数据。这部分的核心便是如何计算出载客距离。因为我们得到的数据中已经含有了出租车在某一时刻的经纬度，所以我们可以计算相邻的两个点的距离，然后再把载客期间所有的距离加起来，这样我们便能得到我们所需要的载客距离，只要这个载客距离小于300，便将这一段载客记录删除。
第三步要做的便是将静止时间小于60秒（这个是可以更改的）的删除。在做这一步时需要对数据的时间的格式进行转换，转换成可以运算的格式，这里我们用到的函数是to_datetime()函数。然后需要做的便是从经纬度未发生变化的第一个数据开始，直到经纬度发生变化的时间结束，计算这一段内的时间，最后的结果是用秒来显示。当我的静止时间之和大于60秒时，便将这一段数据删除。
最后一步也是最难的一步，需要将所有不属于上海市的数据进行删除。这里我刚开始使用的方法是将经纬度反解析成地址，然后判断解析出来的地址是否含有上海市字样。如果不含有上海市，那么便将该数据进行删除，但是实际操作中发现该操作太慢了，进行了四个小时，却只进行了100个左右文件。于是便放弃了这个想法。然后我用到的方法是从网上找出上海市的经纬度的边界，这样边界便形成了一个框，只要是经纬度在这个框内的数据都留下，否则便都删除。
最后我们得到的结果是只有2290个出租车的文件。

import pandas as pd
from math import radians, cos, sin, asin, sqrt
from pathlib import * 
import requests
path='上海市出租车GPS数据/Taxi_070220/'
def sudu(wenjian):  #对于速度大于60km/h的异常数据进行处理
    for i in range(len(wenjian)):
        if  wenjian[4][i]>60:
            wenjian.drop([i],axis=0,inplace=True)
    return wenjian
def geodistance(lng1,lat1,lng2,lat2):   #由经纬度计算实际的距离
    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
    dlon=lng2-lng1
    dlat=lat2-lat1
    a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 
    distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径，6371km
    distance=round(distance/1000,3)
    return distance            #结果的单位为km
def zaikedis(wenjian):  #去掉载客距离小于300米的数据
    t=0
    de=[]
    b=0
    ti=[]
    dis=0
    for i in range(len(wenjian)):
        if wenjian[6][i]!=0 and t==0:
            t=1
            b=i
        if t==1:
            ti.append(i)
            if b!=i:
                d=geodistance(wenjian[2][i-1],wenjian[3][i-1],wenjian[2][i],wenjian[3][i])
                dis=dis+d
        if wenjian[6][i]==0 and t==1:
            t=0
            ti.remove(i)
            
            if dis<0.3:
                de.append(ti)
            ti=[]
            dis=0
    for i in range(len(de)):
        wenjian.drop(de[i],axis=0,inplace=True)
    return wenjian
def parktime(wenjian):   #将静止状态超过60秒的数据删除
    de=[]
    ti=[]
    t=0
    t1=0
    t2=0
    for i in range(1,len(wenjian)):
        x1=wenjian.iloc[i][2]-wenjian.iloc[i-1][2]
        x2=wenjian.iloc[i][3]-wenjian.iloc[i-1][3]
        if x1==0 and x2==0 and t==0:
            t=1
            t1=wenjian[1][i]
        if t==1:
            ti.append(i)
        if t==1 and (x1!=0 or x1!=0):
            t=0
            ti.remove(i)
            t2=wenjian.iloc[i][1]
            t3=t2-t1
            if t3.seconds>120:
                de.append(ti)
            ti=[]
    for i in range(len(de)):
        wenjian.drop(de[i],axis=0,inplace=True)
    return wenjian
def zhuanhuan(s):
    s=str(s)
    s=list(s)
    for i in range(len(s)):
        if s[i]=="\\":
            s[i]='/'
    s="".join(s)
    return s
def transform(location):
    parameters = {'coordsys':'gps','locations': location, 'key': 'bc6dfa0e10696090a74ed56d7976f254'}
    base = 'http://restapi.amap.com/v3/assistant/coordinate/convert'
    response = requests.get(base, parameters)
    answer = response.json()
    return answer['locations']
def geocode(location):
        parameters = {'location': location, 'key': 'bc6dfa0e10696090a74ed56d7976f254'}
        base = 'http://restapi.amap.com/v3/geocode/regeo'
        response = requests.get(base, parameters)
        answer = response.json()
        return answer['regeocode']['formatted_address']
def jingwei(wenjian):
    d=[]
    for i in range(len(wenjian)):
        s=str(wenjian[2][i])+','+str(wenjian[3][i])
        d.append(s)
    return d
p=Path(path)
for x in p.iterdir():
    s=x
    s=zhuanhuan(s)
    wenjian=pd.read_csv(s,encoding='utf-8',header=None)
    wenjian.drop_duplicates(subset=[2,3,4,5,6],keep=False,inplace=True)
    wenjian.reset_index(drop=True,inplace=True)
    d=[]
    for i in range(len(wenjian)):
        if wenjian.iloc[i,2]>122.134497 or wenjian.iloc[i,2]<120.856692:
            d.append(i)
            continue
        if wenjian.iloc[i,3]>31.869483 or wenjian.iloc[i,3]<30.596680:
            d.append(i)
            continue
    wenjian.drop(d,axis=0,inplace=True)
    if len(wenjian) != 0:
        wenjian.to_csv('上海数据/'+s[29:],index=False)

出租车轨迹基本特征分析

出租车路径长度的分布

从上图我们可以分析出来，其实出租车路程的分布接近于一个正态分布，基本上每天跑的路程集中在300-400km，这个数据在现实生活中是很合理的。

出租车载客路径长度的分布

从2-3我们可以看出载客的出租车路程是总路程的一半，一半集中在200km左右，这里我们可以看到0-50km的区间很密，说明也有些出租车载客路程比较短，但是也不会低于1000米，这些应该是出现在2-2中那些本身路程就很短的出租车中。

出租车载客次数分布

从图中可以得知，其实没有接客的出租车数量也不在少数，可以看到0上的车辆数最多的，也就是说出租车有很多接不到客人的，但是这也并不影响一个整体去世的分布，其实大多数的接客次数在20-50次之间，这也符合一个正态分布。

出租车空载速度的分布

从上图我们可以看出如果司机空载时候的速度基本是在30km/h附近，这个速度在城市中来回穿行应该是大多数出租车司机进行寻找客人的速度。

载客平均速度出租车分布

相比较空载的速度，载客的速度就会慢一点，在25km/h附近，因此大多数司机还是在保护乘客安全的情况下，速度不会太快，至于出现在将近60km/h的两个数据，这里暂时还想不出是什么原因，很可能是在当时的情境下比较紧急，而且又因为当时是在过年，所以在市区不超速的情况下，尽量的快速行驶。

出租车数据的可视化

经纬度范围划分

绘制Taix_105轨迹热图，发现轨迹中集中与分散的轨迹对比强烈。

统计出租车数据的经纬度范围

预处理后的出租车数据的经纬度范围：
经度：120.85°-122.00°，相差115个0.01°
纬度：30.59°-31.87°，相差128个0.01°
对出行线路起始点经纬度的处理，将上海的经纬度范围以0.01度为间隔划分，将落在同一个方格内的点按照同一个经纬度来处理，我们进行经纬度划分：
经度：120.85-122.00，共计115个划分单元。
纬度：30.59-31.87，共计128个划分单元。

客人上下车站点统计

建立二维矩阵（115*128）统计出租车数据的客人上下车数据，矩阵的值表示客人上下车次数，保存为taxi_on.csv（客人上车站点数据）和taxi_off.csv（客人下车站点数据）。

全时段上下车站点可视化

使用plotly库画出地理地图，经纬度按0.01°划分，使用taxi_on.csv数据，按照客人上下车权值选择相应的颜色，绘制客人上车的站点图。

客人上车的站点可视化

客人下车的站点可视化

分时段上下车站点可视化

按照之前绘制客人上下车站点图的方法，只是改变数据统计的时间范围。可视化时间段分为0–6点、6–9点、9–12点、12–15点、15–18点、18–24点。

分时段统计上下车站点数据

对于是上下车站点的数据，取出数据中的时间列，得到小时hour，并根据hour进行判断统计。输出每个时间段的最值：

为了可视化对比强烈，以taxi_on_24.max()==321.0为参考绘制站点图，站点颜色分为11级，以321为第11级，则应该对所有的值除32对应到相应的颜色。

分时段上车站点可视化

分时段下车站点可视化

出行线路可视化

把每一辆出租车的轨迹作为一条出行线路进行可视化。

可视化所有出租车的出行线路，共2289条轨迹，发现极其密集，已经覆盖了所有地图上的线路，并且轨迹图加载缓慢，我们选出其中100条进行可视化。

还是很密集，放大查看：

# 可视化
mapbox_access_token = (
    'pk.eyJ1IjoibHVrYXNtYXJ0aW5lbGxpIiwiYSI6ImNpem85dmhwazAy'
    'ajIyd284dGxhN2VxYnYifQ.HQCmyhEXZUTz3S98FMrVAQ'
) # 此处的写法只是为了排版，结果为连接在一起的字符串
layout = go.Layout(
    autosize=True,
    mapbox=dict(
        accesstoken=mapbox_access_token,
        center=dict(
            lat=31.23, #上海市北纬度
            lon=121.47 #上海市北经度
        ),
        pitch=0,
        zoom=7,  #显示地图区域大小
    ),
)

Color = ('#1ffae3','#0bf7f7','#3effc3','#76ff89','#b6ff47','#ebff12','#ffea00','#ffa600','#ff6509','#ff2c19','#ff0523') #颜色渐变
data = [] #绘制数据
marked = set()
cnt=0

for i in range(len(trace)):
    cnt=i%11
    data.append(
        go.Scattermapbox(
            lon=trace[i][0], #站台经度
            lat=trace[i][1], #站台纬度
            mode='lines',#markers：点模型，  lines：线模型，  'markers+lines'：点-线模型
            #name='trace'+str(i), #线路名称，显示在图例（legend）上
            #text='trace'+str(i), #各个点的名称，鼠标悬浮在点上时显示

            # 设置标记点的参数
            marker=go.scattermapbox.Marker(size=5, color=Color[cnt])
        )
    )

fig = dict(data=data, layout=layout)
#py.iplot(fig) #直接显示地图
py.offline.plot(fig, filename='lines_all.html') #生成html文件并打开

出行网络拓扑结构分析

在这里对出行线路起始点经纬度的处理，参照了上面对出行网络可视化的工作，将上海的经纬度范围以0.01度为间隔划分，将落在同一个方格内的点按照同一个经纬度来处理（但在后面绘制的拓扑结构图各个点整齐排列，显得不自然）。
经过计算与划分，我们得到经纬度范围为：
经度：120.85-122.00，共计115个划分单元。
纬度：30.59-31.87，共计128个划分单元。
由于数据量比较大，在下面对出行网络拓扑图的绘制过程中，以两小时为间隔来绘制拓扑图，共绘制了12张图，并根据图片进行分析；对于一些其他的特征，则是针对所有数据进行分析。

节点数

经过统计，节点数为2135个。

边数

经过统计，边数为35237条。

权重及其分布

从上面的图片中，我们可以看到，在DataFrame中的‘trip’属性列为各条边的权重，表示从起点到终点所拥有的路径数。

节点度的分布

接下来分析出行网络各节点的度数，因为图像为有向图,所以我们绘制了入度、出度和度的分布情况。

度

入度

出度

总体上来看，还是上海市地域范围内中心部分的出租车交通流量更大一些。

拓扑结构图的可视化

对于拓扑图中节点的分布问题，我们利用的是Gephi软件中的geo layout插件，会按照节点的经纬度来对节点进行规划布局。

社团划分

在上面绘制的拓扑结构图基础上，我们对网络进行了社区划分，采用的是经典的walktrap算法，实现工具使用的是python库igraph。

通过对节点进行社区划分，我们重新生成了存储节点的csv文件，在文件中增加了一列表示节点所属于的社区编号。

再利用Gephi软件绘制拓扑结构图，节点的分布使用了geo layout插件，相关的设置如下：

节点颜色按照社区归属进行分配，节点大小按照节点度数的多少进行设置。线的颜色采用黄色，方便观察。

## 出行网络
import pandas as pd
import numpy as np
import os
from scipy.sparse import csc_matrix
from datetime import datetime
path='code/上海出租车数据_杨/上海数据/'
files=os.listdir(path)
Max_longitude=0
Min_longitude=140
Max_latitude=0
Min_latitude=45
for f in files:
    df=pd.read_csv(path+f)
    if(max(df['2'])>Max_longitude):
        Max_longitude=max(df['2'])
    if(min(df['2'])<Min_longitude):
        Min_longitude=min(df['2'])
    if(max(df['3'])>Max_latitude):
        Max_latitude=max(df['3'])
    if(min(df['3'])<Min_latitude):
        Min_latitude=min(df['3'])
for step in range(0,24,2):
    OD=np.zeros((115*128,115*128))
    start=datetime(2007,2,20,step,0,0)
    end=datetime(2007,2,20,step+1,59,59)
    for f in files:
        df=pd.read_csv(path+f)
        flag=0
        t=0
        for i in range(df.shape[0]-1):
            if datetime.strptime(df.iloc[i,1],'%Y-%m-%d %H:%M:%S')>start and datetime.strptime(df.iloc[i,1],'%Y-%m-%d %H:%M:%S')<=end:
                if(flag==0 and df.iloc[i,6]==0 and df.iloc[i+1,6]==1):
                    flag=1
                    SNo_longitude=int((df.iloc[i,2]-120.85)/0.01)
                    SNo_latitude=int((df.iloc[i,3]-30.59)/0.01)
                    s=SNo_latitude*115+SNo_longitude
                    t=t+1
                if(flag==1 and df.iloc[i,6]==1 and df.iloc[i+1,6]==0):
                    flag=0
                    ENo_longitude=int((df.iloc[i,2]-120.85)/0.01)
                    ENo_latitude=int((df.iloc[i,3]-30.59)/0.01)
                    e=ENo_latitude*115+ENo_longitude
                    t=t+1
                if(t==2):
                    OD[s,e]=OD[s,e]+1.0
                    t=0
    csc=csc_matrix(OD)
    coo = csc.tocoo(copy=False)
    edge=pd.DataFrame({'start':coo.row,'end':coo.col,'trip':coo.data})[['start','end','trip']].sort_values(['start','end']).reset_index(drop=True)
    tag=[]
    for i in range(edge.shape[0]):
        if(abs(int(edge.iloc[i,0])-int(edge.iloc[i,1]))>1 and abs(int(edge.iloc[i,0])-int(edge.iloc[i,1]))!=116):
            tag.append(i)
    edge=edge.iloc[tag,:]
    edge.to_csv('edge_'+str(step)+'-'+str(step+2)+'.csv',index=False)
    print(step)

## 节点
for step in range(0,24,2):
    edge=pd.read_csv('edge_'+str(step)+'-'+str(step+2)+'.csv')
    point=list(edge['Source'])+list(edge['Target'])
    point=set(point)
    point=list(point)
    point.sort()
    c={'point':point}
    df=pd.DataFrame(c)
    longtitude=[]
    latitude=[]
    for i in range(df.shape[0]):
        latitude.append((df.iloc[i,0]//115)*0.01+30.59)
        longtitude.append((df.iloc[i,0]%115)*0.01+120.85)
    lat={'lat':latitude}
    long={'long':longtitude}
    lat=pd.DataFrame(lat)
    long=pd.DataFrame(long)
    df=pd.concat([df,long,lat],axis=1)
    df.to_csv('point_'+str(step)+'-'+str(step+2)+'.csv',index=False)

## ID转换
def get_keys(d, value):
    return [k for k,v in d.items() if v == value]
for step in range(0,24,2):
    edge=pd.read_csv('edge_'+str(step)+'-'+str(step+2)+'.csv')
    point=pd.read_csv('point_'+str(step)+'-'+str(step+2)+'.csv')
    dic={}
    for i in range(point.shape[0]):
        dic[i]=point.iloc[i,0]
    for i in range(edge.shape[0]):
        edge.iloc[i,0]=get_keys(dic,edge.iloc[i,0])[0]
        #print(get_keys(dic,edge.iloc[i,0])[0])
        edge.iloc[i,1]=get_keys(dic,edge.iloc[i,1])[0]
        #print(get_keys(dic,edge.iloc[i,1])[0])
    print(step)
    edge.to_csv('new/new_edge_'+str(step)+'-'+str(step+2)+'.csv',index=False)
    point['Id']=point.index
    point=point.drop(['point'],axis=1)
    point.to_csv('new/new_point_'+str(step)+'-'+str(step+2)+'.csv',index=False)
dic={}
for i in range(point.shape[0]):
    dic[i]=point.iloc[i,0]
for i in range(edge.shape[0]):
    edge.iloc[i,0]=get_keys(dic,edge.iloc[i,0])[0]
    #print(get_keys(dic,edge.iloc[i,0])[0])
    edge.iloc[i,1]=get_keys(dic,edge.iloc[i,1])[0]
    #print(get_keys(dic,edge.iloc[i,1])[0])
    print(i)
    point.to_csv('new_point.csv',index=False)

## 社区结构
for step in range(0,24,2):
    edge=pd.read_csv('new/new_edge_'+str(step)+'-'+str(step+2)+'.csv')
    point=pd.read_csv('new/new_point_'+str(step)+'-'+str(step+2)+'.csv')
    g=ig.Graph(1)
    g.add_vertices(point.shape[0]-1)
    weights=[]
    for i in range(edge.shape[0]):
        g.add_edge(edge.iloc[i,0],edge.iloc[i,1])
        weights.append(edge.iloc[i,2])
    se = g.community_walktrap(weights, steps=4)
    ss = se.as_clustering()
    result = [[]for i in range(ss.__len__())]
    for i in range(0, ss.__len__(), 1):
        result[i] = (ss.__getitem__(i))
    flag=0
    trip={}
    for club in result:
        for t in club:
            trip[t]=flag
        flag=flag+1
    trip_sort=sorted(trip.items(),key=lambda x:x[0])
    l=[]
    for i in trip_sort:
        l.append(i[1])
    L={'community':l}
    df=pd.DataFrame(L)
    df=pd.concat([point,df],axis=1)
    df.drop(['Id'],axis=1)
    df['Id']=df.index
    df.to_csv('new/point_community_'+str(step)+'-'+str(step+2)+'.csv',index=False)

文章中的代码并不是全部，只是部分关键代码，由于水平有限，可能写得代码并不是很好，欢迎大家指点批评。