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

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