Hello 艾瑞巴蒂!
今天为大家带来的是本公众号第二篇文章,读完本文你将学会:
利用python进行网约车订单数据时空分布特性探索性挖掘
利用python进行空间自相关的检验并构建地理加权回归(GWR)模型
说到地理加权回归,相信大家肯定不会陌生。作为一种先进的空间数据分析技术,地理加权回归能够充分捕捉空间关系的非平稳性。举个简单的不恰当的例子,我们要对中国各个城市的奢侈品消费量与人均收入进行建模。正常的的理解是人均收入越高,奢侈品消费量就越大,在全国各个城市都应该是这种关系(这也正是全局模型的前提假设)。但事实真的是这样吗?现实情况可能是在一些比较张扬的地方(比如我们大东百),奢侈品销量与人均收入是负相关的,不管多穷说啥都得整一个LV;而在某些省份,大家赚得多买的也多,奢侈品销量与人均收入正相关。此时,全局模型不再适用,于是我们今天的主角GWR模型就闪亮登场。GWR估计的模型系数在空间上是变化的,拿上面的例子来说,GWR估计的系数在东百可能是负的,而在别的地区就是正的,明显更加科学,这也正是GWR的强大之处。
本期教程的研究对象是成都市滴滴的订单数据(来自盖亚数据公开计划),思路主要来源于下面这篇Sustainability杂志上的文章《揭示时空维度下城市建成环境对网约车出行的变化性影响——一个关于成都市的探索分析》(渣翻译,轻喷)。
分析中用到的地图数据以及该论文,大家向公众号后台发送关键字:GWR数据,即可获取。滴滴的数据需要自行申请,百度搜索滴滴盖亚计划即可。下面让我们正式开始。
一般来说,拿到时空数据的第一步就是看看其在时空上的分布,论文中分析的是11月3号的订单数据,由于我没拿到11月3号的数据,本文中由11月2号的数据替代。
首先是时间维度,论文中统计的结果如下:
可以看到,半夜出行量少,从5点开始一路飙增,中午13点达到最多。那么11月2号的情况是怎样呢,让我们来探索一下。
## 首先还是熟悉的调包import pandas as pdimport datetimeimport numpy as npimport matplotlib.pyplot as plt import pytz import mathimport geopandas as gpdfrom shapely.geometry import Point,Polygonfrom fiona.crs import from_epsgplt.rcParams['font.sans-serif']=['Arial Unicode MS']plt.rcParams['axes.unicode_minus']=False
## 导入订单数据df = pd.read_csv('order_20161102')df.columns = ['id','start_time','end_time','ori_lng','ori_lat','des_lng','des_lat']
看看数据长啥样:
可以看到,时间是unix时间戳,需要转换一下,此外坐标系并不是wgs,而是火星坐标系,后面我们也要转换一下。
进行数据的转换:
## 将unix时间戳转换为datetimedef time_trans(arr): return datetime.datetime.fromtimestamp(int(arr['start_time']))df['start_time'] = df.apply(time_trans,axis=1)### 提取每一单的开始时间的时段df['hour'] = df['start_time'].dt.hour ### 把经纬度转化为数字格式(原先是字符串)df['ori_lng'] = pd.to_numeric(df['ori_lng'],errors='coerce')df['ori_lat'] = pd.to_numeric(df['ori_lat'],errors='coerce')df['des_lng'] = pd.to_numeric(df['des_lng'],errors='coerce')df['des_lat'] = pd.to_numeric(df['des_lat'],errors='coerce')
画图,先看一下24小时的分布:
ridership = df.groupby('hour')['id'].count() ### 用groupby提取每一小时的出行量ridership24 = ridership.shift(-1)ridership24[23] = ridership[0] ### 把数据重新排序time = []for i in range(1,24,1): time.append('{0}:00'.format(i))plt.figure(dpi=100)plt.bar(np.arange(0,24,1),ridership24,color='red')plt.xticks(np.arange(0,24,1),time,rotation=90)plt.show()
上面是订单24小时的分布情况,总体情况跟论文插图差不多,我们再画个跟论文插图一样的&#x