合作者：guyueCT

城市公交网络分析与可视化

数据爬取与处理

公交站点信息爬取

第一步，从青岛市公交线路的百度百科（https://baike.baidu.com/item/%E9%9D%92%E5%B2%9B%E5%85%AC%E4%BA%A4/5016038?fr=aladdin）中爬取青岛市区的各路公交线路名称，整理成txt文档，用于后续爬取线路信息工作。

第二步，申请高德地图的开发者key，选择web端的服务。
高德key

利用所得到的key和一些公交线路信息构成url，提交给高德地图，并爬取返回的与路线相关的信息，关键代码如下：

url = 'https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=559bdffe35eec8c8f4dae959451d705c&output=json&city={}&offset=1&keywords={}&platform=JS'.format(city,line)
r = requests.get(url).text
rt = json.loads(r)#转换为字典

其中，city代表要查询的城市，line代表要查询的公交线路；我们以青岛市的321路公交线路为例，查看爬取的结果。

爬取的结果已经转换成了字典，其中一些对我们地图可视化分析有帮助的有公交线路的名称、各公交站点的名称与坐标、各公交线路的轨迹坐标以及公交线路中各站点的顺序。
第三步，将爬取得到的相关信息进行整理、汇总成表格形式，以便后续后续可视化及分析。

import requests
import pandas as pd
import json
import re
f = open("raw_stop.txt","r")   #文件内容为上述统计的站点名称
str = f.read()     #将txt文件的所有内容读入到字符串str中
f.close()   #将文件关闭
def get_dt(city,line):
    url = 'https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=559bdffe35eec8c8f4dae959451d705c&output=json&city={}&offset=1&keywords={}&platform=JS'.format(city,line)
    r = requests.get(url).text
    rt = json.loads(r)
    print(rt)
    try:
        if rt['buslines']:
            if len(rt['buslines']) == 0:  #有名称没数据
                print('no data in list..')
            else:
                dt = {
    }
                dt['line_name'] = rt['buslines'][0]['name'] 
                dt['total_price'] = rt['buslines'][0]['total_price']

                st_name = []
                st_coords = []
                st_sequence = []
                for st in rt['buslines'][0]['busstops']:
                    st_name.append(st['name'])
                    st_coords.append(st['location'])
                    st_sequence.append(st['sequence'])

                dt['station_name'] = st_name
                dt['station_coords'] = st_coords
                dt['sequence'] = st_sequence
                dm = pd.DataFrame([dt])
                return dm
        else:
            pass
    except:
        print('error..try it again..')
        time.sleep(2)
        get_dt(city,line)
stations=re.split(',| |\n',str)
stations=list(filter(None,stations))
for sta in stations:
    if sta== '1路':
        data=get_dt('青岛',sta)
    else:
        dm=get_dt('青岛',sta)
        data=pd.concat([data,dm],ignore_index=True)
for i in range(data.shape[0]):
    coord_x=[]
    coord_y=[]
    for j in data['station_coords'][i]:
        coord_x.append(eval(re.split(',',j)[0]))
        coord_y.append(eval(re.split(',',j)[1]))
    a=[[data['line_name'][i]]*len(data['station_coords'][i]),coord_x,coord_y,data['station_name'][i],data['sequence'][i]]
    df=pd.DataFrame(a).T
    if i==0:
        df1=df
    else:
        df1=pd.concat([df1,df],ignore_index=True)
df1.columns=['line_name','coord_x','coord_y','station_name','sequence']
df1.to_csv('bus_inf.csv',index=None,encoding='utf_8_sig')

公交线路轨迹爬取

与上面的方法一致，爬取得到的信息字典中，有’polyline’项，项中存储的为各公交线路的轨迹坐标，同理将这些信息进行整理。

import requests
import pandas as pd
import json
import re
f = open("raw_stop.txt","r")   #设置文件对象
str = f.read()     #将txt文件的所有内容读入到字符串str中
f.close()   #将文件关闭
def get_dt(city,line):
    url = 'https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=559bdffe35eec8c8f4dae959451d705c&output=json&city={}&offset=1&keywords={}&platform=JS'.format(city,line)
    r = requests.get(url).text
    rt = json.loads(r)
    try:
        if rt['buslines']:
            if len(rt['buslines']) == 0:  #有名称没数据
                print('no data in list..')
            else:
                dt = {
    }
                dt['line_name'] = rt['buslines'][0]['name'] 
                dt['polyline'] = rt['buslines'][0]['polyline']
                dm = pd.DataFrame([dt])
                return dm
        else:
            pass
    except:
        print('error..try it again..')
        time.sleep(2)
        get_dt(city,line)
stations=re.split(',| |\n',str)
stations=list(filter(None,stations))
for sta in stations:
    if sta== '1路':
        data=get_dt('青岛',sta)
    else:
        dm=get_dt('青岛',sta)
        data=pd.concat([data,dm],ignore_index=True)
data.to_csv('bus_inf_trace.csv',index=None,encoding='utf_8_sig')  

坐标转换

像谷歌地图等地图，使用的是wgs84坐标，实际坐标与地图坐标相同，而高德、腾讯地图等为了保障国家安全，在wgs84坐标的基础上对国内范围进行加密，称为火星坐标，百度地图为了商业目的，在火星坐标的基础上进行了二次加密。在这里，我们是通过高德地图采取的公交线路信息，因此需要对火星坐标进行转化，否则绘制的路线会发生偏移，出现公交线路跑到海边的问题。

# -*- coding: utf-8 -*-
import json
import math
import pandas as pd

x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 扁率

def gcj02towgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return lng, lat
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [lng * 2 - mglng, lat * 2 - mglat]


def transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
        0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
        0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    if lng < 72.004 or lng > 137.8347:
        return True
    if lat < 0.8293 or lat > 55.8271:
        return True
    return False

t_x=[]
t_y=[]
df=pd.read_csv('bus_inf.csv')
for i in range(len(list(df['coord_x']))):
    [X,Y]=gcj02towgs84(list(df['coord_x'])[i],list(df['coord_y'])[i])
    t_x.append(X)
    t_y.append(Y)
t_x=pd.DataFrame(t_x)
t_y=pd.DataFrame(t_y)
df['coord_x']=t_x
df['coord_y']=t_y
df.to_csv('bus_inf_wgs84.csv',index=None,encoding='utf_8_sig')    

城市公交网络可视化

利用ArcMap实现地图的可视化

参考这篇文章
知乎-城市路网实时路况爬取与ArcGIS可视化

利用plotly实现地图可视化

# 导入库及数据
import requests
import time
import numpy as np
import math
import plotly.offline as py
import plotly.graph_objs as go
from plotly.offline import init_notebook_mode
import pandas as pd
import re
import json

bus_inf_stations = pd.read_csv("bus_inf_wgs84.csv")
bus_inf_trace=pd.read_csv('bus_inf_trace.csv')
init_notebook_mode(connected=True) #如果使用jupyter需要加上这一句，把地图显示在页面上
# 绘制青岛公交轨迹 + 公交站点图
mapbox_access_token = (
    'pk.eyJ1IjoibHVrYXNtYXJ0aW5lbGxpIiwiYSI6ImNpem85dmhwazAy'
    'ajIyd284dGxhN2VxYnYifQ.HQCmyhEXZUTz3S98FMrVAQ'
) # 此处的写法只是为了排版，结果为连接在一起的字符串
layout = go.Layout(
    autosize=True,
    mapbox=dict(
        accesstoken=mapbox_access_token,
        center=dict(
            lat=36.08997714, #青岛市北纬度
            lon=120.3899010 #青岛市北经度
        ),
        pitch=0,
        zoom=11,
    ),
)

color = ('blue', 'green', 'yellow', 'purple', 'orange', 'red', 'violet', 
         'navy', 'crimson', 'cyan', 'magenta', 'maroon', 'peru','gray',
        'brown','darkgray','deeppink','greenyellow','lemonchiffon','lightblue',
        'limegreen','mediumaquamarine','mintcream','orangered','powderblue',
        'seagreen','springgreen','tomato','steelblue','skyblue') #可按需增加
data = [] #绘制数据
marked = set()
cnt = 0
i=0
while(i<len(bus_inf_trace)): #遍历每一条线路
    #plots_name = [] #站台名称
    coord_lon=[]
    coord_lat=[]
    line_name=bus_inf_trace['line_name'][i]
    poly=re.split(';',bus_inf_trace['polyline'][i]) #当前线路的轨迹
    j=0
    while(j<len(poly)):#遍历所有轨迹点
        (lon,lat)=gcj02towgs84(eval(re.split(',',poly[j])[0]), eval(re.split(',',poly[j])[1]))
        coord_lon.append(lon)  #站台经度集
        coord_lat.append(lat)  #站台纬度集
        j=j+1
        
    data.append(
        go.Scattermapbox(
            lon=coord_lon, #站台经度
            lat=coord_lat, #站台纬度
            mode='lines',#markers：点模型，  lines：线模型，  'markers+lines'：点-线模型
            name=line_name, #线路名称，显示在图例（legend）上
            text=plots_name, #各个点的名称，鼠标悬浮在点上时显示
            
            # 设置标记点的参数
            marker=go.scattermapbox.Marker(size=5, color=color[cnt])
        )
    )
    #marked.add(uid) #添加已绘制线路的uid
    #marked.add(line['pair_line_uid']) #添加已绘制线路反向线路的uid
    cnt = (cnt + 1) % len(color)
    i=i+1
        
while(i<len(bus_inf_stations)): # <8501
    plots_name = [] #站台名称
    coord_lon=[]
    coord_lat=[]
    line_name=bus_inf_stations['line_name'][i]
    while(line_name==bus_inf_stations['line_name'][i]):
        plots_name.append(bus_inf_stations['station_name'][i])
        coord_lon.append(bus_inf_stations['coord_x'][i])
        coord_lat.append(bus_inf_stations['coord_y'][i])
        if i==8500:
            break
        i=i+1
        
    data.append(
        go.Scattermapbox(
            lon=coord_lon, #站台经度
            lat=coord_lat, #站台纬度
            mode='markers',#markers：点模型，  lines：线模型，  'markers+lines'：点-线模型
            name=line_name, #线路名称，显示在图例（legend）上
            text=plots_name, #各个点的名称，鼠标悬浮在点上时显示
            
            # 设置标记点的参数
            marker=go.scattermapbox.Marker(size=5, color='black')
        )
    )
    #marked.add(uid) #添加已绘制线路的uid
    #marked.add(line['pair_line_uid']) #添加已绘制线路反向线路的uid
    cnt = (cnt + 1) % len(color)
    i=i+1
    

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

公交路线基本特征分析

公交线路的平均长度

在爬取的数据字典中，每条公交线路都含有键名为’ distance’的键值对，按照上面爬取公交线路的思路，爬取相关数据，得到各公交线路的长度，并计算线路平均长度。

经过计算，我们得到了线路的平均线路长度为
19.452331013986004km

公交线路的平均站点数

经过分析可以得出，要统计站点数只需要统计站点信息表格中相同’line_name’的行数就可以了。

经过计算，我们得到了线路的平均站点数大约为
30（30.42909090909091）

公交线路的平均站距

平均站距的计算，可以直接用平均路线长度除以平均站点数得到。
经过计算，可以得到了线路的平均站距大约为
0.6609898713513099km

公交线路的平均直线系数

根据公式：
直线系数=公交线路的长度/公交线路的直线长度

可以求出各路线的直线系数(排除环形路线)，其中，公交线路的直线长度通过公交线路的起始和终点站的地理坐标来求解得到，再对这些直线系数求平均，能够得到公交线路的平均直线系数。


经过计算，我们得到了公交线路的平均直线系数为
2.306061353782556

#平均站点数
import pandas as pd
df=pd.read_csv('bus_inf_wgs84.csv')
dic1={
    }
line=df['line_name']
line=line.drop_duplicates()
for L in line:
    temp=df[df['line_name']==L]
    dic1[L]=temp.shape[0]
L = len(dic1)
S = sum(dic1.values())
A = S / L

#路线平均长度
import requests
import pandas as pd
import json
import re
f = open("raw_stop.txt","r")   #设置文件对象
str = f.read()     #将txt文件的所有内容读入到字符串str中
f.close()   #将文件关闭
def get_dt(city,line):
    url = 'https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=559bdffe35eec8c8f4dae959451d705c&output=json&city={}&offset=1&keywords={}&platform=JS'.format(city,line)
    r = requests.get(url).text
    rt = json.loads(r)
    try:
        if rt['buslines']:
            if len(rt['buslines']) == 0:  #有名称没数据
                print('no data in list..')
            else:
                dt = {
    }
                dt['line_name'] = rt['buslines'][0]['name'] 
                dt['distance'] = eval(rt['buslines'][0]['distance'])
                dm = pd.DataFrame([dt])
                return dm
        else:
            pass
    except:
        print('error..try it again..')
        time.sleep(2)
        get_dt(city,line)
stations=re.split(',| |\n',str)
stations=list(filter(None,stations))
for sta in stations:
    if sta== '1路':
        data1=get_dt('青岛',sta)
    else:
        dm=get_dt('青岛',sta)
        data1=pd.concat([data1,dm],ignore_index=True)
B=sum(data1['distance'])/data1.shape[0]

#平均站距
C=B/(A-1)
C

#平均直线系数
def get_dt_l(city,line):
    url = 'https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=559bdffe35eec8c8f4dae959451d705c&output=json&city={}&offset=1&keywords={}&platform=JS'.format(city,line)
    r = requests.get(url).text
    rt = json.loads(r)
    try:
        if rt['buslines']:
            if len(rt['buslines']) == 0:  #有名称没数据
                print('no data in list..')
            else:
                if rt['buslines'][0]['start_stop'] != rt['buslines'][0]['end_stop']:
                    dt={
    }
                    dt['line_name'] = rt['buslines'][0]['name']
                    dt['start_stop'] = rt['buslines'][0]['start_stop']
                    dt['end_stop'] = rt['buslines'][0]['end_stop']
                    for st in rt['buslines'][0]['busstops']:
                        if st['name']==dt['start_stop']:
                            dt['start_coord']=st['location']
                        if st['name']==dt['end_stop']:
                            dt['end_coord']=st['location']
                    dm = pd.DataFrame([dt])
                    return dm
        else:
            pass
    except:
        print('error..try it again..')
        time.sleep(2)
        get_dt_l(city,line)
stations=re.split(',| |\n',str)
stations=list(filter(None,stations))
for sta in stations:
    if sta== '1路':
        data2=get_dt_l('青岛',sta)
    else:
        dm=get_dt_l('青岛',sta)
        data2=pd.concat([data2,dm],ignore_index=True)
start_coord_x=[]
start_coord_y=[]
end_coord_x=[]
end_coord_y=[]
for i in range(data2.shape[0]):
    end_coord_x.append(eval(re.split(',',data2['end_coord'][i])[0]))
    start_coord_x.append(eval(re.split(',',data2['start_coord'][i])[0]))
    end_coord_y.append(eval(re.split(',',data2['end_coord'][i])[1]))
    start_coord_y.append(eval(re.split(',',data2['start_coord'][i])[1]))
a=[start_coord_x,start_coord_y,end_coord_x,end_coord_y]
df=pd.DataFrame(a).T
df.columns=['start_coord_x','start_coord_y','end_coord_x','end_coord_y']
D=pd.concat([data2,df],axis=1)
D=D.drop(['start_coord','end_coord'],axis=1)
D.to_csv('bus_inf_start_end.csv',index=None,encoding='utf_8_sig')

from math import sin, asin, cos, radians, fabs, sqrt
EARTH_RADIUS = 6371  # 地球平均半径，6371km
def hav(theta):
    s = sin(theta / 2)
    return s * s
def get_distance_hav(lat0, lng0, lat1, lng1):
    "用haversine公式计算球面两点间的距离。"
    # 经纬度转换成弧度
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lng0 = radians(lng0)
    lng1 = radians(lng1)

    dlng = fabs(lng0 - lng1)
    dlat = fabs(lat0 - lat1)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))
    return distance
dic2={
    }
for i in range(len(D)):
    dic2[D['line_name'][i]]=get_distance_hav(D['start_coord_y'][i],D['start_coord_x'][i],D['end_coord_y'][i],D['end_coord_x'][i])
    if dic2[D['line_name'][i]]<0.5:
        del dic2[D['line_name'][i]]

dic3={
    }
for i in dic2.keys():
    dic3[i]=data1[data1['line_name']==i]['distance']/dic2[i]
for i in dic2.keys():
    dic3.values()   

S=0
num=0
for i in dic2.keys():
    S=S+list(dic3[i])[0]
    num=num+1
N=S/num

公交换乘网络搭建

换乘网络的搭建需要确定各公交站点的距离矩阵，距离矩阵如下：

矩阵中值为1，表示横纵坐标的站点之间能够通过一条公交线路到达；值为0，表示无法通过一条公交线路到达，需要换乘。
根据换乘矩阵，搭建包含所有的站点的换乘网络如下：

由于站点数量过多，无法很好的显示，我们只取十条公交线路进行绘制换乘网络，结果如下：

import pandas as pd
import numpy as np
df=pd.read_csv('bus_inf_wgs84.csv')
line=df['line_name']
line=line.drop_duplicates()
network=np.zeros((line.shape[0],line.shape[0]))
network=pd.DataFrame(network,index=line)
network.columns=line
for i in range(len(list(line))):
    for j in range(i,len(list(line))):
        temp1=df[df['line_name']==list(line)[i]]
        temp2=df[df['line_name']==list(line)[j]]
        T1=set(list(temp1['station_name']))
        T2=set(list(temp2['station_name']))
        if T1&T2 != set():
            network[list(line)[i]][list(line)[j]]=network[list(line)[j]][list(line)[i]]=1.0  
stage=[]
for i in list(line):
    for j in list(line):
        if network[j][i] == 1.0:
            stage.append((j,i))
import networkx as nx
import matplotlib.pyplot as plt
%matplotlib inline
#plt.figure(figsize=(8, 8.5))
G=nx.Graph()
G.add_nodes_from(list(line))
G.add_edges_from(stage)
nx.draw(G,with_labels=False,edge_color='b', node_color='g',node_size=100,font_family ='YouYuan',spring_layout=True,font_size=8)
plt.show()

公交换乘网络分析

节点数和边数

节点数为公交线路的长度为275，边数为距离矩阵中值‘1’的数量为10063。

节点的邻居数（度）

可以理解成为，距离矩阵每行或每列数值‘1’的个数，经过计算得到的结果为

节点度的分布（直方图）

网络的平均路径长度

给定一个网络G(n, e), n是节点数量，e是链接数量。对于网络中的任意一对节点source s和target t，我们可以计算它们之间的最短距离ds,t。将这些距离相加除以n(n−1)，就可以得到网络的平均最短路径长度。再这里，我们直接利用networkx的函数average_shortest_path_length()来直接进行计算，得到结果为
2.1946303180503923

参考资料

1. Plotly绘图库官网 https://plot.ly/
2. 使用networkx计算平均最短路径长度
   https://computational-communication.com/average-shortest-path-length/
3. python通过经纬度计算两点之间距离https://blog.csdn.net/zzmtkj/article/details/94591792
4. python复杂网络结构可视化——matplotlib+network
   https://zhuanlan.zhihu.com/p/36700425