【转载】利用python实现地理加权回归(GWR)与网约车订单数据挖掘

共计 11062 个字符，预计需要花费 28 分钟才能阅读完成。

利用 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 小时的分布情况，总体情况跟论文插图差不多，我们再画个跟论文插图一样的：

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()

这个足够以假乱真了吧，不难看出，总体趋势和 11 月 3 号是一致的，具体数值有小小的差别。

下面我们来统计空间分布，首先看看论文分析的区域：

文章把研究区域划成了 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_BW

area_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 大神及其弟子写的文章，十分硬核。这里我根据论文中的模型设置，选择了固定带宽，高斯核函数，AICc 法。（文章中没说他选的是固定带宽还是可变带宽，不过我看他最后给出了带宽的值为 2560，看样子应该是固定带宽）

结果如下：

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。