人工智能与地理大数据实验--出租车GPS数据—时空大数据Python处理基础（二）

环境：Windows 10 专业版 + Python 3.9.1 + anaconda 2020（ 4.8.2）

系列文章：

人工智能与地理大数据实验--出租车GPS数据—时空大数据Python处理基础（一）

人工智能与地理大数据实验--出租车GPS数据—时空大数据Python处理基础（二）

目录

五、源码流程分析

（一）、出租车数据的时间完整性评估

（二）、出租车数据的空间完整性评估

（三）、出租车订单出行特征分析

六、实验结果

（一）、出租车数据的时间完整性评估

（二）、出租车数据的空间完整性评估

（三）、出租车订单出行特征分析

五、源码流程分析

（一）、出租车数据的时间完整性评估

1.时间字段的处理

利用字符串str方法，把Series对象转换为pandas内置的StringMethods对象，并使用StringMethods对象的slice方法截去字符串切片，传入开始与结束位置的索引，提取前两位字符，从而获取小时信息。

# 把Stime列转换为StringMethods对象，进行切片操作，再赋值给Hour列
data['Hour'] = data['Stime'].str.slice(0,2)
data

2.数据量的小时集计

同一小时内产生的GPS数据，Hour字段中会拥有相同的值，可以根据这一字段对数据采用groupby方法进行集计操作，再利用count方法指定VehicleNum列统计每小时的数据量。

# 分组并统计各组数量
Hourcount = data.groupby('Hour')['VehicleNum'].count()
# 更改Series的列名，并通过reset_index将Series变成DataFrame
Hourcount = Hourcount.rename('count').reset_index()
Hourcount

3.数据量时间分布的折线图绘制

通过matplotlib库将数据绘制为简单的折线图与柱状图。

# 1.创建图表
# 导入包
import matplotlib.pyplot as plt
# 创建一个图，图的尺寸为8×4in，dpi为300
fig = plt.figure(1,(8,4),dpi = 300)
# 在图中创建子图
# 111分别表示：创建共一个子图，子图的布局为1行1列
ax = plt.subplot(111)
# 2.在图上画
# 绘制折线图，分别传入节点的x坐标与y坐标，'k-'定义了黑色实线
plt.plot(Hourcount['Hour'],Hourcount['count'],'k-')
# 绘制散点图，分别传入节点的x坐标与y坐标，'k.'定义了黑色散线
plt.plot(Hourcount['Hour'],Hourcount['count'],'k.')
# 绘制柱状图，分别传入节点的x坐标与y坐标
plt.bar(Hourcount['Hour'],Hourcount['count'])
# 3.调整图中元素
# 加y轴标题
plt.ylabel('Data volume')
# 加x轴标题
plt.xlabel('Hour')
# 调整x轴标签
plt.xticks(range(24),range(24))
# 加图标题
plt.title('Hourly data volume')
# 设置y轴范围
plt.ylim(0,30000)
# 显示图
plt.show()

（二）、出租车数据的空间完整性评估

1.出租车GPS数据空间分布栅格图

（1）研究区域栅格生成

用Python代码实现研究范围内栅格，并存储为GeoDataFrame形式。

①研究范围行政区划边界的读取

使用Geopandas读取Shapefile文件

# 导入Geopandas
import geopandas as gpd
# 读取shp格式地理数据文件
sz = gpd.read_file(r'D:\Project\Jupyter_Project\Data\data\sz\sz.shp',encoding = 'utf-8')
# 绘制地理数据文件
sz.plot()

查看读取进来的变量类型及内容

type(sz)

geopandas.geodataframe.GeoDataFrame

sz

②研究范围内栅格的生成

首先定义研究范围与栅格的大小，计算栅格的△lon与△lat。

# 导入math包
import math
# 划定栅格划分范围
lon1 = 113.75194
lon2 = 114.624187
lat1 = 22.447837
lat2 = 22.864748
# 取得左下角经纬度
latStart = min(lat1,lat2)
lonStart = min(lon1,lon2)
# 定义栅格大小，单位为米
accuracy = 500
# 计算栅格的经纬度增加量大小△Lon和△Lat，地球半径取6371004米
deltaLon = accuracy * 360 / (2 * math.pi * 6371004 * math.cos((lat1 + lat2) * math.pi / 360))
deltaLat = accuracy * 360 / (2 * math.pi * 6371004)
deltaLon,deltaLat

(0.004872614089207591, 0.004496605206422906)

然后，定义一个测试GPS点的经纬度，计算其对应的栅格编号与中心点经纬度。

# 定义一个GPS点测试栅格化
testlon = 114
testlat = 22.5
# 计算该GPS点对应的栅格编号
LONCOL = divmod(float(testlon) - (lonStart - deltaLon / 2), deltaLon)[0]
LATCOL = divmod(float(testlat) - (latStart - deltaLat / 2), deltaLat)[0]
# 计算该GPS点对应的栅格中心点经纬度
# 格子编号*格子宽+起始横坐标=格子中心横坐标
HBLON = LONCOL * deltaLon + lonStart
HBLAT = LATCOL * deltaLat + latStart
LONCOL,LATCOL,HBLON,HBLAT

(51.0, 12.0, 114.00044331854959, 22.501796262477075)

接下来计算每个栅格的四个顶点坐标，循环生成栅格，生成整个研究范围的栅格。

import geopandas as gpd
from shapely.geometry import Polygon
# 定义空的GeoDataFrame表，再往里加栅格
data = gpd.GeoDataFrame()
# 定义空的list，后面循环一次就往里面加东西
LONCOL_list = []
LATCOL_list = []
geometry_list = []
HBLON_list = []
HBLAT_list = []
# 计算行列要生成的栅格数量
# lon方向是lonsum个栅格
lonsnum = int((lon2-lon1) / deltaLon) + 1
# lat方向是latsum个栅格
latsnum = int((lat2-lat1) / deltaLat) + 1
for i in range(lonsnum):
    for j in range(latsnum):
        # 第i列，第j行的栅格中心点坐标
        HBLON = i * deltaLon + lonStart
        HBLAT = j * deltaLat + latStart
        # 用周围的栅格推算三个顶点的位置
        HBLON_1 = (i + 1 ) * deltaLon + lonStart
        HBLAT_1 = (j + 1) * deltaLat + latStart
        # 生成栅格的Polygon形状
        grid_ij = Polygon([
            (HBLON - deltaLon / 2, HBLAT - deltaLat / 2),
            (HBLON_1 - deltaLon / 2, HBLAT - deltaLat / 2),
            (HBLON_1 - deltaLon / 2, HBLAT_1 eltaLat / 2)])
        # 把生成的数据都加入到前面定义的空list里面
        LONCOL_list.append(i)
        LATCOL_list.append(j)
        HBLON_list.append(HBLON)
        HBLAT_list.append(HBLAT)
        geometry_list.append(grid_ij)
# 为GeoPandas文件的每一列赋值为刚刚的list
data['LONCOL'] = LONCOL_list
data['LATCOL'] = LATCOL_list
data['HBLON'] = HBLON_list
data['HBLAT'] = HBLAT_list
data['geometry'] = geometry_list
data

对创建的栅格绘制查看其空间分布。

data.plot(edgecolor = (0,0,0,1),linewidth = 0.2)

利用GeoDataFrame的intersects方法，输入行政区划的Polygon几何图形，将行政区划边界外的栅格剔除。

# 将前面读取的GeoDataFrame用unary_union方法合并为一个Polygon图形作为研究范围
roi = sz.unary_union
# 筛选出研究范围的栅格
grid_sz = data[data.intersects(roi)]
grid_sz.plot(edgecolor = (0,0,0,1),linewidth = 0.2)

将栅格数据保存至本地

grid_sz.to_file(r'D:\Project\Jupyter_Project\Data\data\grid_sz1.josn', driver = 'GeoJSON')

（2）GPS数据的栅格对应与集计

GPS数据以经纬度坐标的形式记录，计算得到每个GPS点对应的栅格编号。然后，再基于栅格标号统计每个栅格内GPS数据量。

import pandas as pd
# 读取数据
data = pd.read_csv('D:/Project/Jupyter_Project/Data/data/TaxiData-Sample.csv',header = None)
data.columns = ['VehicleNum','Stime','Lng','Lat','OpenStatus','Speed']
# 数据对应的栅格经纬度编号
data['LONCOL'] = ((data['Lng'] - (lonStart - deltaLon / 2))/deltaLon).astype('int')
data['LATCOL'] = ((data['Lat'] - (latStart - deltaLat / 2))/deltaLat).astype('int')
# 集计栅格数量
data_distribution = data.groupby(['LONCOL','LATCOL'])['VehicleNum'].count().rename('count').reset_index()
# 剔除不在研究范围内的OD记录
data_distribution = data_distribution[(data_distribution['LONCOL']>=0)&(data_distribution['LATCOL']>=0)&(data_distribution['LONCOL']=0)&(oddata['SLATCOL']>=0)&
                (oddata['ELONCOL']>=0)&(oddata['ELATCOL']>=0)&
                (oddata['SLONCOL']