利用python实现地理加权回归(GWR)与网约车订单数据挖掘
Hello 艾瑞巴蒂！
今天为大家带来的是本公众号第二篇文章，读完本文你将学会：
利用python进行网约车订单数据时空分布特性探索性挖掘
利用python进行空间自相关的检验并构建地理加权回归(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小时的分布情况，总体情况跟论文插图差不多，我们再画个跟论文插图一样的：
plt.figure(dpi=100)ridership_sel = []for i in range(1,24,2): ridership_sel.append(ridership24[i])time = []for i in range(1,24,2): time.append('{0}:00-{1}:00'.format(i,i+1))plt.bar(np.arange(0,24,2),ridership_sel,color='red',width=1.5)plt.xticks(np.arange(0,24,2),time,rotation=45)plt.title('自己画的图')plt.ylabel('Ridership')plt.xlabel('Time')plt.show()这个足够以假乱真了吧
下面我们来统计空间分布，首先看看论文分析的区域：
文章把研究区域划成了17*17的500m方格。那我们也把它画出来，首先明确一下研究的区域，对着b图，我在osm上大致确定了总体的位置：
大方框的经纬度如左上所示，下面让我们用geopandas把它画出来。
### 生成网格（这种方式很粗略，要想精确的还是到arcgis里操作比较好）y = np.linspace(30.6521,30.7290,18)x = np.linspace(104.0422,104.1284,18)grid = []for m,i in enumerate(x[:17]): for k,j in enumerate(y[:17]): cord = np.array([i,j,x[m+1],j,x[m+1],y[k+1],i,y[k+1]]).reshape(4,2) grid.append(Polygon(cord))## 创建网格的geodataframe，具体的方法第一篇文章都讲过area = gpd.GeoDataFrame(grid)area = area.reset_index()area.columns=['index','geometry']area.crs={'init':'epsg:4326'}chengdu = gpd.GeoDataFrame.from_file('./成都街道/chengdu.shp')base = chengdu.plot(figsize=(10,10),facecolor='None',edgecolor='Black')area.plot(ax=base,facecolor='None',edgecolor='red')plt.ylim(30.63,30.75)plt.xlim(104.02,104.15)plt.show()结果如上，差不多是那个意思。下一步就是统计每个格子里的出行量，画下面这张图：
首先要对网约车数据的坐标进行转换，这里用到转换的代码是我百度的，感谢这位csdn上的博主。
### 进行火星坐标系的转换x_pi = 3.14159265358979324 * 3000.0 / 180.0pi = 3.1415926535897932384626 # πa = float(6378245.0) # 长半轴ee = 0.00669342162296594323 # 扁率def gcj02towgs84(lng, lat): """ GCJ02(火星坐标系)转GPS84 :param lng:火星坐标系的经度 :param lat:火星坐标系纬度 :return: """ 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## ————————————————## 版权声明：本文为CSDN博主「GeoDoer」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。## 原文链接：https://blog.csdn.net/summer_dew/article/details/8072343### 应用转换函数把火星坐标转换为wgs84df13 = df[df['hour']==13] ### 筛选13点的出行订单df13 = df13.reset_index().drop(columns='index')for i in range(len(df13)): df13.loc[i,'ori_lng'] = gcj02towgs84(df13.loc[i,'ori_lng'],df13.loc[i,'ori_lat'])[0] df13.loc[i,'ori_lat'] = gcj02towgs84(df13.loc[i,'ori_lng'],df13.loc[i,'ori_lat'])[1]### 创建 geodataframedef point(arr): return Point(arr['ori_lng'],arr['ori_lat']) df13['geometry'] = df13.apply(point,axis=1)df13 = gpd.GeoDataFrame(df13) ### 用到在上一篇文章中提到的方法，df13.crs = {'init':'epsg:4326'}## 画图查看一下上车点在空间上的分布base = chengdu.plot(figsize=(10,10),facecolor='None',edgecolor='Black')area.plot(ax=base,facecolor='None',edgecolor='red')df13.plot(ax=base,markersize=1,color='blue')plt.ylim(30.63,30.75)plt.xlim(104.02,104.15)plt.show()可以看到，左下区域还是比较密（春熙路一带），接着我们集计一下。
### 将上车点与网格进行空间连接，并统计每个网格内的出行量join_ridership = gpd.sjoin(area,df13,op='contains',how='left') join_ridership.dropna(inplace=True)area = pd.merge(area,join_ridership.groupby('index')['id'].count().to_frame().reset_index(),on='index',how='outer')area['id'] = area['id'].fillna(0)area.columns = ['index','geometry','ridership']area.plot(figsize=(10,10),column='ridership',cmap='YlOrRd',edgecolor='Black',legend=True,legend_kwds={'shrink':0.72,'pad':-0.02,'label':'Ridership'})plt.axis('off')plt.show()总体上的跟插图差不多，我画的是13点这个时段出行的，文中是全天的，并且我画的是连续的热力图，文章里的做了断点处理分成了6个等级。
然后我们也统计一下各类POI在各个网格内的分布情况，这里我偷了个懒，文章中用了逐步回归，把不显著的变量都剔除了，最后的情况如下：
在13:00-14:00时间段，用到的解释变量有土地利用混合度、公交站、居民区、餐饮、购物和公司的POI。
下面我们导入相关poi数据：
bus = pd.read_csv('./CSV版本/交通设施服务.csv',encoding='gbk')residential = pd.read_csv('./CSV版本/商务住宅.csv',encoding='gbk')catering = pd.read_csv('./CSV版本/餐饮.csv',encoding='gbk')shopping = pd.read_csv('./CSV版本/购物.csv',encoding='gbk')company = pd.read_csv('./CSV版本/公司企业.csv',encoding='gbk')entertainment = pd.read_csv('./CSV版本/体育休闲服务.csv',encoding='gbk')residential = residential[(residential['中类']=='住宅区')]company = company[(company['中类']=='公司')|(company['中类']=='公司企业')]bus = bus[bus['中类']=='公交车站']### 去除文件中缺失信息poi = [bus.copy(),catering.copy(),company.copy(),residential.copy(),shopping.copy(),entertainment.copy()]for i in range(len(poi)): poi[i].dropna(inplace=True)def point(arr): return Point(arr['WGS84_经度'],arr['WGS84_纬度']) ## 定义转换函数for i in range(len(poi)): poi[i]['geometry'] = poi[i].apply(point,axis=1) #识别经纬度，转换点数据 poi[i] = gpd.GeoDataFrame(poi[i]) poi[i].crs = {'init':'epsg:4326'} ### 设置坐标系为wgs84for i in range(len(poi)): poi[i] = gpd.sjoin(poi[i],area,op='within',how='left') poi[i].dropna(inplace=True) ### 进行空间连接 将poi连接到相应格子### 开始计算格子内各类poi的数量stat = area.loc[:,['index']]for x in range(len(poi)): stat = pd.merge(stat,poi[x].groupby('index')['名称'].count().to_frame().reset_index(),on='index',how='outer')stat.columns = ['index','bus','catering','company','residential','shopping','entertainment']stat.fillna(0,inplace=True) ### 填充缺失值上面的代码都是第一篇文章介绍过的，这里直接拿来用。
还有一项指标是土地利用混合度，计算公式如下:
下面是计算的代码：
### 计算土地利用混合度stat['landuse_mix'] = 0for i in range(len(stat)): total = 0 for x in range(1,7): total += stat.iloc[i,x] for x in range(1,7): if total == 0: ### 有的网格里面一个poi也没有，把他的混合度定为1 stat.loc[i,'landuse_mix'] = 1 else: stat.loc[i,'landuse_mix'] += (stat.iloc[i,x]/total)**2stat.describe().T看一下统计的结果
上面是我们计算的，下面是论文中计算的：
大体上差不多，但土地利用混合度区别较大，原因在于文章中计算这个混合度用到了所有种类的poi，而我为了省事只用了建模时使用的6种，导致我们计算的土地利用混合度和论文中的有较大不同。
接着我们看看各类指标在空间上是如何分布的：
area = pd.merge(area,stat,on='index',how='outer')fig,ax = plt.subplots(2,4,figsize=(30,20))fig.tight_layout(h_pad=-50,w_pad=-2)name = ['ridership','bus','catering','company','residential','shopping','entertainment','landuse_mix']for i in range(0,2): for x in range(0,4): area.plot(ax=ax[i][x],column=name[x+4*i],cmap='YlOrRd',edgecolor='Black',legend=True,legend_kwds={'shrink':0.39,'pad':-0.02,'label':name[x+3*i]}) ax[i][x].set_title(name[x+4*i],fontsize=25) ax[i][x].axis('off')封面图就做出来了。
为了检验变量的空间自相关，论文统计了Moran指数。这时候有同学会问，arcgis能统计莫兰指数，python行不行啊？答案当然是可以的，下面上代码。
这里我们用到了pysal这个包，全名python spatial analysis library，计量经济学方面的空间分析工具这个包里面基本都有，后面的GWR我们也是根据这个包实现的，有兴趣的同学可以自己研究一下别的功能。
### 计算莫兰指数from pysal.lib.weights import Queenfrom pysal.explore.esda.moran import Moranw_queen = Queen.from_dataframe(area) ### 构建权重矩阵 queen矩阵moran = pd.DataFrame(index=name)for i in name: moran.loc[i,'Moran_I'] = Moran(np.array(area[i]),w_queen).I moran.loc[i,'Z_score'] = round(Moran(np.array(area[i]),w_queen).z_norm,3) moran.loc[i,'P_value'] = round(Moran(np.array(area[i]),w_queen).p_norm,3)moran可以看到，所有变量均显著，并且莫兰指数大于0，说明有空间集聚的情况。下面是论文中关于莫兰指数的统计：
跟我们计算的对比，总体上是那个意思，但是数值有的差了不少，应该是因为poi用的不一样，我这个poi好像是18年的，用在16年的滴滴数据上十分不合适
下面就是重头戏了，变量都有了，我们开始用GWR建模。
先调包，从pysal里导入GWR相关函数：
from pysal.model.mgwr.gwr import GWR from pysal.model.mgwr.sel_bw import Sel_BWarea_proj = area.to_crs(from_epsg(32648)) #### 投影转换u = area_proj.centroid.x ## 获取x投影坐标v = area_proj.centroid.y ## 获取y投影坐标g_y = area['ridership'].values.reshape((-1,1))g_X = area[['landuse_mix','bus','residential','catering','shopping','company']].valuesg_coords = list(zip(u,v))g_X = (g_X - g_X.mean(axis=0)) / g_X.std(axis=0) ## 标准化g_y = g_y.reshape((-1,1))g_y = (g_y - g_y.mean(axis=0)) / g_y.std(axis=0) ## 标准化上面首先进行一下投影，这样得到的结果会比较准确，然后对x和y进行标准化。
开始建模：
gwr_selector = Sel_BW(g_coords, g_y, g_X,fixed=True,kernel='gaussian')gwr_bw = gwr_selector.search(criterion='AICc')
gwr_results = GWR(g_coords, g_y, g_X,gwr_bw,fixed=True,kernel='gaussian').fit()解释一下上面的几个参数，fixed代表核函数的带宽是可变还是固定，kernel代表核函数的类型，主要有gaussian（高斯核函数）和bisquare（双平方）核函数，criterion是选择带宽的方式，主要有CV（交叉验证）和AICc（赤池信息准则），如果想深入了解的话大家可以看看Fotheringham大神及其弟子写的文章，十分硬核
结果如下：
pysal会同时给出全局OLS模型与GWR模型的估计参数，这里我们得到的带宽为8031，即每个点周围会建立半径为8公里的圆，其中的点会附上不同的权重，半径以外的点的权重非常小基本可以忽略不计。
下面我们来画一下各个变量的系数在空间上的分布（所谓的空间不平稳性）：
name = ['intercept','landuse_mix','bus','residential','catering','shopping','company']for i in range(7): area['gwr_'+name[i]] = gwr_results.params[:,i]color = ['Blues','Greens','YlOrBr','RdPu','Reds','Purples']fig,ax = plt.subplots(2,3,figsize=(30,20))name = ['landuse_mix','bus','residential','catering','shopping','company']for i in range(0,2): for x in range(0,3): area.plot(ax=ax[i][x],column='gwr_{}'.format(name[i*3+x]),cmap=color[i*3+x],k=6,scheme='fisher_jenks',edgecolor='Grey',legend=True) ax[i][x].set_title('gwr_{}'.format(name[i*3+x]),fontsize=25) ax[i][x].axis('off')可以看到，不同空间位置的系数确实不太一样。
下面我们来看看论文中的模型结果：
咦，细心的你会发现，怎么我们模型结果跟论文完全不一样啊，论文里GWR模型r方达到了0.82，较全局模型提高了0.02；而我们计算的模型里GWR和全局模型r2都一样，别的诊断信息也基本都一样，根本没什么提升；并且参数的分布，我们计算的GWR和论文中的GWR结果完全不同，这是为什么？
要想回答这个问题，我们首先要看一下我们估计的全局模型的信息：
其中x1（土地利用混合度）和x2（公交站点）这两个变量并不显著，放在GWR模型里并不合适，因此我们计算的GWR模型实际上是有问题的。我们为了省事直接用了论文的结论，选了6个变量。但由于我们的poi数据和论文中使用的并不一致，论文中六个变量并在全局模型中都是显著的，而我们的变量有的并不显著，导致了我们计算的GWR存在严重的问题。此外还有可能我们数据因变量与自变量之间的关系并不是空间不平稳的，此时用OLS全局模型即可，不需要GWR。
至此，本篇文章就结束啦，感谢大家的支持，我们下篇文章再见