1、某市某条公交线路数据的爬取分析

1、创建属于自己的API的key值,我创建的应用类型为出行

Android 高德公交 高德地图公交出行方案_ci

2、高德地图开发文档的内容解析

1、定位到Web服务API

Android 高德公交 高德地图公交出行方案_数据_02


这里包含有第一步怎么获取Key值。

1、地理编码:

将详细的结构化地址转换为高德经纬度坐标。且支持对地标性名胜景区、建筑物名称解析为高德经纬度坐标。结构化地址举例:北京市朝阳区阜通东大街6号转换后经纬度116.480881,39.989410
地标性建筑举例:天安门转换后经纬度:116.397499,39.908722

2、逆地理编码:

将经纬度转换为详细结构化的地址,且返回附近周边的POI、AOI信息。例如:116.480881,39.989410 转换地址描述后:北京市朝阳区阜通东大街6号

Android 高德公交 高德地图公交出行方案_Android 高德公交_03

3、查询的URL地址

(https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=bc6bab115306dde0c6b07db965106b7a&ouput=json&city=南京&offset=1&keywords=88路&platfrom=JS)参数设置有待在研究

输入网址,获得如下内容,将内容拷贝到json在线解析网址

Android 高德公交 高德地图公交出行方案_数据_04


解析后的内容

Android 高德公交 高德地图公交出行方案_json_05

4、数据抓取
1、

此次我们需要抓取公交线路名称、公交线路数据、票价信息、公交站点名称、公交站点坐标

#1、获取南京市88路公交线路信息
url='https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=bc6bab115306dde0c6b07db965106b7a&ouput=json&city=南京&offset=1&keywords=88路&platfrom=JS'
url=re.get(url).text
data=json.loads(url)

dt={}
dt['line_name']=data['buslines'][0]['name']
dt['polyline']=data['buslines'][0]['polyline']
dt['total_price']=data['buslines'][0]['total_price']
da_name=[]
da_cor=[]
for i in data['buslines'][0]['busstops']: #由于busstops里面有很多条字典,循环访问每一条字典
    da_name.append(i['name']) #放入每条字典的公交站名
    da_cor.append(i['location'])#放入每条字典的公交站名位置
    
dt['bus_name=']=da_name
dt['bus_cor']=da_cor
dm=pd.DataFrame([dt])#将字典转化为一个二维数组
print(dm)

全国城市公交网数据,我们定位到南京市公交线路网页,我们按照数字查找进行数据爬取

Android 高德公交 高德地图公交出行方案_数据_06


在Xpath Helper中定位到我们要查村的线路名信息

Android 高德公交 高德地图公交出行方案_json_07

lst = []
for i in range(1,10):
    url = 'http://{}.gongjiao.com/lines_{}.html'.format('nanjing',i)#把城市和数字查找进行参数化
#         print(url)
    r = re.get(url).text        
    et = etree.HTML(r)
    line = et.xpath('//div[@class="list"]//a/text()')
    for l in line:
        lst.append(l.split('南京')[1]) #将每一条公交线路名称加入到结果列表里

这是提取的结果

Android 高德公交 高德地图公交出行方案_Android 高德公交_08

2、上述两个代码的整和(添加了异常处理操作)
#3 获取某个城市的所有公交线路数据
def get_lines(cityE,city):
    lst=[]
    for i in range(1,10):
        url='http://{}.gongjiao.com/lines_{}.html'.format(cityE,i)
#        print(url)
        r=re.get(url).text
        et=etree.HTML(r)
        line=et.xpath('//div[@class="list"]//a/text()')
#        print(line)
        for l in line:
            lst.append(l.split(city)[1])
    return lst

def get_roadinfo(city,line):
    url='https://restapi.amap.com/v3/bus/linename?s=rsv3&extensions=all&key=bc6bab115306dde0c6b07db965106b7a&ouput=json&city={}&offset=1&keywords={}&platfrom=JS'.format(city,line)
    url=re.get(url).text
    data=json.loads(url)
    #异常数据的处理
    try:
        if data['buslines']:
            print("data is avaliable!")
            if len(data['buslines'][0])!=0:    
                dt={}
                dt['line_name']=data['buslines'][0]['name']
                dt['polyline']=data['buslines'][0]['polyline']
                dt['total_price']=data['buslines'][0]['total_price']
                da_name=[]
                da_cor=[]
                for i in data['buslines'][0]['busstops']:
                    da_name.append(i['name'])
                    da_cor.append(i['location'])
                    
                dt['bus_name=']=da_name
                dt['bus_cor']=da_cor
                dm=pd.DataFrame([dt])#将字典转化为一个二维数组
                dm.to_csv('{}_lines.csv'.format(cityE),mode='a',header=False,encoding='gbk')
                
            else:
                print('nodata in list')
    except:
        print("There's an error here!")
        time.sleep(2) #休眠两秒后
        get_roadinfo(city,line) #再次执行程序
        

if __name__=='__main__':
    t1=time.time()
    cityE='chongqing'
    city='重庆'
    road_lines=get_lines(cityE,city)
#    print(len(road_lines))
    print('{}:有{}条公交线路'.format(city,len(road_lines)))
    for i,line in enumerate(road_lines):
        print('id{}:{}==================='.format(i,line))
        get_roadinfo(city,line)
    t2=time.time()
    print("耗时{}秒".format(t2-t1))
3、爬取的城市交通线路数据进行清洗与分析
'''
爬取的城市交通线路数据进行清洗与分析
'''
import pandas as pd
import numpy as np
import requests as re
import json
import time
from  lxml import etree
#1将爬取得数据进行数据加载,然后对公交站点数据进行清洗
#print(pd.read_csv('nanjing_lines.csv').head(0))
df_lines = pd.read_csv('nanjing_lines.csv',names=["id","line_name",'line',"price","station_name","station_cor"])
df_sts = df_lines[["station_cor","station_name"]]
df_sts['station_cor'] = df_sts['station_cor'].apply(lambda x:x.replace('[','').replace(']','').replace('\'','').split(', ')) #DataFrame中apply的用法,将str转化为列表
df_sts['station_name'] = df_sts['station_name'].apply(lambda x:x.replace('[','').replace(']','').replace('\'','').split(', '))
#np.hstack(df_sts['station_name'].repeat(list(map(len,df_sts['station_name']))))#把所有数据展成一行
st_all = pd.DataFrame(\
                     np.column_stack((\
                                     np.hstack(df_sts['station_cor'].repeat(list(map(len,df_sts['station_cor'])))),
                                     np.hstack(df_sts['station_name'].repeat(list(map(len,df_sts['station_name'])))),
                                     )),
                      columns=['station_cor', 'station_name']                     
                     )
#print(st_all)
da_all=st_all.drop_duplicates() #去除重复的公交站点
da_all.to_csv('南京公交站点清洗.csv')

#2、线路数据清洗
df=df_lines [["line_name",'line']]
df['line']=df['line'].apply(lambda x:x.split(';'))
print(df.head(2))
df.to_json("南京公交站线路数据清洗.json")
4、转换坐标系-创建矢量文件

 

sf = shapefile.Reader("shapefiles/blockgroups.shp")

2、读取Shapefile元数据

sf.shapeType

形状类型由shapefile规范定义的0到31之间的数字表示,并在下面列出。重要的是要注意,编号系统有几个尚未使用的保留编号,因此现有形状类型的编号不是连续的:
NULL = 0
点= 1
折线= 3
多边形= 5
多点= 8
POINTZ = 11
POLYLINEZ = 13
POLYGONZ = 15
MULTIPOINTZ = 18
点数= 21
折线= 23
POLYGONM = 25
多点= 28
MULTIPATCH = 31
3、添加点形状

>>> w = shapefile.Writer('shapefiles/test/point')
>>> w.field('name', 'C')

>>> w.point(122, 37) 
>>> w.record('point1')

>>> w.close()
添加LineString形状
>>> w = shapefile.Writer('shapefiles/test/line')
>>> w.field('name', 'C')

>>> w.line([
...			[[1,5],[5,5],[5,1],[3,3],[1,1]], # line 1
...			[[3,2],[2,6]] # line 2
...			])

>>> w.record('linestring1')

>>> w.close()
添加多边形形状
>>> w = shapefile.Writer('shapefiles/test/polygon')
>>> w.field('name', 'C')

>>> w.poly([
...	        [[122,37], [117,36], [115,32], [118,20], [113,24]], # poly 1
...	        [[15,2], [17,6], [22,7]], # hole 1
...         [[122,37], [117,36], [115,32]] # poly 2
...        ])
>>> w.record('polygon1')

>>> w.close()

完整代码如下

import json
import urllib
import math
import pandas as pd
import numpy as np
import shapefile
# 其他坐标的转换方法见:GCJ02(火星坐标系)转GPS84
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方
def gcj02_to_wgs84(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:
    """
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)

def polit_to_wgs84(x):
    
    lng=float(x.split(',')[0])#数据需要转换为浮点型
    lat=float(x.split(',')[1])
    
    return gcj02_to_wgs84(lng, lat)

def lines_wgs84(x):
    lst = []
    for i in x:
#        print(i)
        lng = float(i.split(',')[0])
        lat = float(i.split(',')[1])
        lst.append(gcj02_to_wgs84(lng,lat))
    return lst

if __name__=='__main__':
   
    #1、加载公交站信息数据
    df_point=pd.read_csv("南京公交站点清洗.csv")
    #创建公交站点几何信息shape file
    df_point["point_to_wgs84"]=df_point['station_cor'].apply( polit_to_wgs84)
    #创建点矢量文件
    w=shapefile.Writer("bus_point.shp")
    w.field("name",'C')#创建字段
    for i in range(len(df_point['point_to_wgs84'])):
        w.point(df_point['point_to_wgs84'][i][0],df_point['point_to_wgs84'][i][1])
        w.record(df_point['station_name'][i].encode('gbk'))
    w.close()
#    
    #2、创建公交线路几何信息
    df_lines=pd.read_json("南京公交站线路数据清洗.json")
    df_lines['lines_wgs84'] = df_lines['line'].apply(lines_wgs84)
    print(df_lines)
    w = shapefile.Writer('bus_line.shp')
    w.field('name', 'C')
    for i in range(len(df_lines)):
        w.line([df_lines['lines_wgs84'][i]])
        w.record(df_lines['line_name'][i].encode('gbk'))
    w.close()
5、 道路数据与交通路线数据的关联矩阵

补充知识点空间连接:GIS数据除了可以通过基本属性表的字段进行关联,也可以通过地理空间位置的拓扑关系进行关联和连接。1.比如,一个点图层,需要添加每一个点所属的行政辖区信息;

Android 高德公交 高德地图公交出行方案_数据_09


Android 高德公交 高德地图公交出行方案_Android 高德公交_10

import arcpy
#道路网做字段融合,保证每条道路一条线;
#道路网络与公交线路做空间链接,按1对多建立关联矩阵;
path=r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata'
gdb_name="nanjing_analysis"
arcpy.CreateFileGDB_management(path, gdb_name) #创建地理数据库

arcpy.env.workspace=r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\day5\amap\basic_roads.gdb'#工作空间
lst=[fc for fc in arcpy.ListFeatureClasses()] #访问工作空间里的要素
arcpy.Merge_management(lst,r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb\all_roads')#道路要素进行合并
arcpy.Dissolve_management('all_roads',r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb\all_roads_dis','NAME')#道路数据按照NAME字段进行数据融合
arcpy.SpatialJoin_analysis("all_roads_dis","bus_line",out_feature_class=r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb\all_roads_join',
                           join_operation='JOIN_ONE_TO_MANY',match_option='WITHIN_A_DISTANCE',search_radius=0.0004)#道路数据和公交线路数据进行一对多空间连接,采用WITHIN_A_DISTANCE匹配策略

#提取道路数据和公交线路数据的关联矩阵.
arcpy.env.workspace=r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb'
fc=r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb\all_roads_join'
fds=["NAME","name_1"] #NAME是道路名称
with open(r'D:\BaiduNetdiskDownload\课程资料toStudents\课程资料toStudents\medata\nanjing_analysis.gdb\nanjing_data.csv','a') as f1:
    with arcpy.da.SearchCursor(fc,fds) as cursor:
        for row in cursor:
            if row[0]==None or row[1]==None:
                pass
            else:
                f1.write(row[0].encode('utf-8') + ','+ row[1].encode('utf-8') +'\n')
6、 拟合的交通路线数据展示与相关系数

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.pyplot import style
style.use('ggplot')
from scipy import optimize
import scipy.stats

'''
城市道路节点度显示
'''
#df = pd.read_csv('D:/BaiduNetdiskDownload/课程资料toStudents/课程资料toStudents/nanjing_data.csv',names=['roads','lines'],encoding='gbk')
#df = df[df['roads'] != ' '] #将roads的空值数据删除,目的是只保留有公交线路的道路数据
#dv = df.groupby(by='roads').count().sort_values(by='lines',ascending=False)#按道路数据进行分组,统计每条道路有多少条公交线路数据,采用降序
#dv= dv.reset_index()#创建数据索引
##print(dv.head(5))
#dv.to_csv("nanjing_gp.csv")
dk=pd.read_csv('nanjing_gp.csv')
dk = dk.reset_index()
dk['index'] = dk['index'] + 1 #让索引值从1开始

plt.figure(figsize=(10,5))
plt.scatter(dk['index'],dk['lines'])
plt.show()

plt.figure(figsize=(10,5))
x = dk['index']
y = dk['lines']
# 可能有更优的拟合曲线公式,为与参考论文一致改公式
def fm(x,a,b):
    return a*x**b

ft1,ft2 = optimize.curve_fit(fm,x,y)#ft1是参数拟合的a,b参数,ft2是参数的协方差矩阵
#print(ft1)
_y = [y for y in fm(x,ft1[0],ft1[1])] #拟合出来的y
s1, i, r, p, s2 = scipy.stats.linregress(y,_y)#为两组测量值计算线性最小二乘回归,s1回归线的斜率,i回归线截距,
#r相关系数,p利用检验统计量的t分布的Wald检验,对原假设为斜率为零的假设检验的双侧p值进行检验,s2估计梯度的标准误差。
# print(r**2)
# print(_y)
#画每条道路对应多条公交线路图
plt.scatter(x,y)
plt.plot(x,_y,color='k',lw=2)#绘制的拟合曲线
plt.text(1800,100,'${%.1f}x^{%.3f}$' % (ft1[0],ft1[1]),fontsize=12)
plt.text(1800,80,'$R^2={%.4f}$' % (r**2),fontsize=12)
plt.show()