利用Python对成都市POI数据进行探索性数据分析

---

        大家过年好！！！悠闲（枯燥）的假期刚刚开始，不知道各位在家里待得怎么样，祝大家鼠年吉祥，身体健康，万事如意！现在疫情严重，大家一定戴好口罩保护好自己！

       首先为大家介绍一下我开这个公众号的初衷（还不是因为假期闲的无聊）。原本立志假期要看n篇文献，但当我打开endnote，一股强烈的睡意袭来，我意识到我之前想多了。但是也不能浪费时光吧，于是想和大家分享一下自己的一些学习收获（十分浅薄的收获），于是就申请了这个公众号。开篇之作，就以成都市POI为例，展示一下如何利用Python对POI数据进行探索性数据分析(Exploratory Data Analysis)。(POI:Point Of Interest，兴趣点）

---

第一步，载入需要的各种包，由于涉及地理数据，用到了geopandas,shapely和伟大的fiona包。

import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point
from fiona.crs import from_epsg
plt.rcParams['font.sans-serif']=['Arial Unicode MS']  ###设置显示中文
plt.rcParams['axes.unicode_minus']=False   

然后读取csv格式的poi数据。

#### 读取poi相关csv文件
traffic = pd.read_csv('./CSV版本/交通设施服务.csv',encoding='gbk')
catering = pd.read_csv('./CSV版本/餐饮.csv',encoding='gbk')
company = pd.read_csv('./CSV版本/公司企业.csv',encoding='gbk')
residential = pd.read_csv('./CSV版本/商务住宅.csv',encoding='gbk')
shopping = pd.read_csv('./CSV版本/购物.csv',encoding='gbk')
government = pd.read_csv('./CSV版本/政府机构机关及社会团体.csv',encoding='gbk')

这里我用到了六类POI，分别是交通设施、餐饮、公司、住宅、购物和机关，先简单地查看一下数据的结构：

traffic.head(5)

可以看到数据中包含了POI的种类，地址，wgs84的经纬度信息（这个数据包含的范围很广，停车场、公交站、车站都在里面，实际使用的时候还是要细化一下分类）。

接着清洗一下数据，去除空值。

### 去除文件中缺失信息
poi = [traffic.copy(),catering.copy(),company.copy(),residential.copy(),shopping.copy(),government.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) #利用apply函数，转换点数据
    poi[i] = gpd.GeoDataFrame(poi[i])
    poi[i].crs = {'init':'epsg:4326'}     ### 设置坐标系为wgs84

读取成都市行政区的shp数据，并转换为投影坐标系。

chengdu = gpd.GeoDataFrame.from_file('./成都街道/chengdu.shp')   ### 读取成都行政区的shp数据
chengdu = chengdu.to_crs(from_epsg(4549))  ## 投影到CGCS2000坐标系
chengdu.plot(facecolor='None',edgecolor='black')  ### 粗略查看一下

还不错，有那个味了。

---

接着对POI点位进行同样的投影转换。

for i in range(len(poi)):
    poi[i] = poi[i].to_crs(from_epsg(4549)) ## 对poi同样投影

至此，POI和行政区的数据预处理工作就全部完成啦，下一步要做的就是将点(POI)与面(行政区)进行空间连接，在arcgis中可以用spatial join工具，相比之下python中的操作就很简单，只需要一行代码即可(geopandas的sjoin函数）。

for i in range(len(poi)):
    poi[i] = gpd.sjoin(poi[i],chengdu,op='within',how='left') ### 进行空间连接 将poi连接到相应街道
    poi[i].dropna(inplace=True)  ### 去除不在行政区范围内的数据

查看一下POI在空间上的分布情况。

name = ['交通设施服务','餐饮','公司企业','商务住宅','购物','政府机关及社会团体']
### 绘图，看看poi在空间上的分布情况
fig,ax = plt.subplots(2,3,figsize=(30,20))
fig.tight_layout(h_pad=-15,w_pad=0)
for i in range(0,2):
    for x in range(0,3):
        chengdu.plot(ax=ax[i][x],facecolor='None',edgecolor='black')
        poi[x+i*3].plot(ax=ax[i][x],color='red',markersize=1)
        ax[i][x].set_title(name[x+3*i],fontsize=25)
        ax[i][x].axis('off')
plt.savefig('poi.jpg',dpi=150)

当当，封面图就做出来啦。当然，只画个图是远远不够的，多少得整点定量吧，要不然怎么忽悠甲方和老师（不小心说漏了什么）。

---

说起定量分析，大家最熟悉的就是线性回归啦。作为最基础也是最广泛的有监督统计学习方式，线性回归深受广大科研工作者和程序员的喜爱。原本想用didi盖亚的数据做个POI密度与网约车出行量之间关系的回归分析，但是申请的数据还没通过审核，那就充分利用现有条件，用多元线性回归分析一下行政区内交通设施POI数量与其他种类POI数量的关系吧。

先统计一下各个行政区内各类poi点位的数量，用到了groupby函数。

### 开始统计行政区内各类poi的数量
stat = chengdu.loc[:,['NAME']]
for x in range(len(poi)):
        stat = pd.merge(stat,poi[x].groupby('NAME')['名称'].count().to_frame().reset_index(),on='NAME',how='outer')
stat.columns = ['NAME','交通设施服务','餐饮','公司企业','商务住宅','购物','政府机关及社会团体']
stat.fillna(0,inplace=True) ### 填充缺失值

然后开始调包，用到的是statsmodels这个模块，话说python关于统计的模块实在是不得行，被R吊打（虽然我还不会R），但statsmodels也能凑合用。

import statsmodels.api as sm
import statsmodels.formula.api as smf
#### 导入statsmodels模块进行分析（调包侠的日常）
results = smf.ols('交通设施服务 ~ 餐饮 + 公司企业 + 商务住宅 + 购物 + 政府机关及社会团体', data=stat).fit() ## 这浓浓的山寨风，怎么那么像R
print(results.summary())

上面就是输出的结果，r2竟然达到了0.9，餐饮、公司企业、商务住宅和购物数量都和交通设施数量都显著正相关，机关数量不显著。不过warnings提醒我们，自变量之间可能存在很强的共线性，statsmodels也可以统计方差膨胀因子(vif)（对不起我不该黑你，statsmodels还是很强大的），下面我们来看看。

from statsmodels.stats.outliers_influence import variance_inflation_factor
### 继续调包，查看各个变量的方差膨胀因子
chengdu_vif = stat.iloc[:,2:]  
vif = pd.DataFrame()
vif['feature'] = chengdu_vif.columns
chengdu_vif['const'] = 1     ### 添加截距，等价于sm.add_constant() 
vif['vif'] = 0
x = np.array(chengdu_vif)
for i in range(len(vif)):
    vif.loc[i,'vif'] = variance_inflation_factor(x,i)
vif

果然，餐饮和购物的vif都大于10，确实是有多重共线性存在。如果要深入分析的话就可以逐步调整选取的自变量，观察vif，以确定最终的自变量，这里就深入不演示了。

---

前面演示了最经典的有监督统计学习方式，不来一个无监督的怎么行。最经典的无监督统计学习方法当然就是K-means聚类啦，话不多说让我们操作起来。

继续调包，大名鼎鼎的sklearn。

from sklearn.cluster import KMeans

首先确定最佳的聚类数量。

## 通过手肘法确认最佳的聚类数
score = []
for i in range(2,13):
    n = i
    kmean = KMeans(n_clusters=n)
    kmean.fit(stat.iloc[:,1:])
    score.append(kmean.inertia_)
plt.figure(figsize=(8,4))
plt.plot(range(2,13),score,c='b')
plt.plot(range(2,13),score,c='b',marker='o')
plt.xticks(np.arange(2,13))
plt.xlabel('cluster',fontsize=13)
plt.ylabel('SSE(sum of the squared errors)',fontsize=13)
plt.title('手肘法确定最佳聚类数量',fontsize=15)
plt.grid()

可以看到，当cluster等于3和4的时候曲线拐得比较厉害，我这里就把聚类数量设为4。

下面就开始聚类分析。

n = 4
kmean = KMeans(n_clusters=n)
kmean.fit(stat.iloc[:,1:])
clus = kmean.predict(stat.iloc[:,1:])
stat['clus'] = clus   ### 将聚类结果添加到stat中

chengdu = pd.merge(chengdu,stat,on='NAME',how='outer')   ### 将聚类结果添加到shp中

chengdu.plot(figsize=(10,10),column=clus,categorical=True,legend=True)  ## 绘图观察一下各类分布

聚类的结果如上，可以看到，第一类只有一个行政区（红色），激起了我的好奇心，查看一下这个行政区是哪一个。

print(chengdu[chengdu['clus']==1]['NAME'])

wow，桂溪街道，这是哪里？百度了一下，原来是高新区，果然是高富帅，鹤立鸡群。

最后，让我们来看看各个类别的中心点具体特征。

### 绘制雷达图，观察各类的特征情况
labels = ['交通设施服务','餐饮','公司企业','商务住宅','购物','政府机关及社会团体']
#每个类别中心点数据
plot_data = kmean.cluster_centers_
#指定颜色
color = ['b', 'g', 'r','black']
# 设置角度
angles = np.linspace(0, 2*np.pi, 6, endpoint=False)
# 闭合
angles = np.concatenate((angles, [angles[0]]))
plot_data = np.concatenate((plot_data, plot_data[:,[0]]), axis=1)

fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111, polar=True) 
for i in range(len(plot_data)):
    ax.plot(angles, plot_data[i], 'o-', color = color[i], label = '类别'+str(i), linewidth=1.5,markersize=6)
ax.set_thetagrids(angles * 180/np.pi, labels, fontproperties="Arial Unicode MS",fontsize=13)
ax.set_title('聚类中心点特征',fontsize=20)
plt.legend(bbox_to_anchor=(0.983,1.1)) # 设置图例位置

果然，高富帅的公司企业数量和餐饮数量远超其他种类的社区。

至此，EDA就算完成啦，第一篇推文，希望大家喜欢，多多关注多多转发！