利用Python爬取轨道交通站点周边POI数据及可视化

Hello！大家好，最近在龟速写论文，一直没有更新，今天为大家带来一篇利用百度API爬取地铁站周边POI数据的教程，基本原理与可达圈那篇一致，可以举一反三。话不多说，下面正式开始。

---

这次分析的地铁站还是宽窄巷子站，因为比较懒，就用了之前osm那篇文章的代码。先调包+把底图画出来：

import osmnx as ox
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point 
import requests
import json
import time
import matplotlib.pyplot as plt
pd.set_option('display.max_columns', None)
pd.set_option('display.max_columns', None)
plt.rcParams['font.sans-serif']=['Arial Unicode MS']
plt.rcParams['axes.unicode_minus']=False

ditie = gpd.read_file('./成都地铁/成都地铁/成都地铁_busstops.shp')
route = gpd.read_file('./成都地铁/成都地铁/成都地铁_buslines.shp')
### 宽窄巷子经纬度 POINT (104.04820 30.66783)
kuanzhai = gpd.GeoDataFrame([Point(104.04820,30.66783)])
kuanzhai.columns = ['geometry']
G = ox.graph_from_point(center_point=(30.66783,104.04820),distance=1000,network_type='drive')
G_gdf = ox.graph_to_gdfs(G)
base = G_gdf[1].plot(figsize=(10,10),edgecolor='grey')
G_gdf[0].plot(ax=base,color='blue')
kuanzhai.plot(marker='*',ax=base,color='red',markersize=200)
west,east = base.get_xlim()   ### 这里获取一下底图的东西南北边界，为下一步爬取边界内的POI做准备
south,north = base.get_ylim()
plt.show()

得到我们本次要分析的区域，跟上篇文章中的一样：

下面开始爬POI数据:

首先看看百度的官方文档怎么说的：

可以看到，必填项有POI的类别，矩形区域左下和右上的坐标，还有ak，值得注意的是，由于我们输入的是WGS84坐标，所以还要加上coord_type=1，否则会默认我们输入的坐标为百度坐标。还有一点比较重要，API最多支持查询400条POI，所以我们的区域不能选的太大，另外请求每次最多返回20条，需要我们识别一下是否需要翻页，下面我们简单的试试：

ak = '你的ak'
url = 'http://api.map.baidu.com/place/v2/search?query=美食&coord_type=1&bounds={},{},{},{}&page_size=20&page_num=0&output=json&ak={}'.format(south,west,north,east,ak)
data = requests.get(url).json()

然后得到返回的数据长这个样子：

是字典的格式，其中‘status'键值对代表查询的状态，'total'代表查询的区域内共有多少个POI，如果大于20的话，我们就需要翻页查询了，'result'里就是我们查询的POI的结果，有地址，坐标，甚至还有电话，这里的detail我没有展开，展开的话还有打分评价之类的。按照这个思路，我们就可以查询地铁站周围的各个类型的POI了，具体操作很简单，写个循环就行，唯一需要注意的就是要判别需不需要翻页，此外，百度api限制了并发时间，每次查询完我们需要time.sleep()一下，具体的代码大家如果需要的话可以后台私信我，有偿提供。这里我一共查了'美食','酒店','购物','生活服务','休闲娱乐','公司企业'这六类POI，结果如下：

再计算一下各类的数量：

df['category'].value_counts()

不愧是景区，4平方公里95个酒店，不过估计这个冬天很大一部分都歇菜了，共克时艰共克时艰。

然后我们把location解析一下，变成GeoDataFrame。

def get_lat(arr):
    return arr['location']['lat']
def get_lng(arr):
    return arr['location']['lng']
df['lat'] = df.apply(get_lat,axis=1)
df['lng'] = df.apply(get_lng,axis=1)
def point(arr):
    return Point(arr['lng'],arr['lat'])
df['geometry'] = df.apply(point,axis=1)
df = gpd.GeoDataFrame(df)

然后plot一下：

base = G_gdf[1].plot(figsize=(10,10),edgecolor='grey')
G_gdf[0].plot(ax=base,color='blue')
df.plot(column='category',legend=True,ax=base,markersize=4,cmap='Paired')
kuanzhai.plot(marker='*',ax=base,color='red',markersize=200)
plt.show()

我擦嘞，怎么偏移了。在想了很久为什么之后，我在官方文档里找到了答案：

原来输出的是百度坐标系。。。。。。（这里吐槽一下混乱的坐标系统，什么gcj坐标，bd坐标，说是为了加密，但也都被破解了，安安静静用84不好吗)。然后我又在github上找到一个老哥写的bd坐标转84的代码，也不知道对不对，试试效果：

for i in df.index:
    df.loc[i,'wgs_lng'] = bd09_to_wgs84(df.loc[i,'lng'],df.loc[i,'lat'])[0]
    df.loc[i,'wgs_lat'] = bd09_to_wgs84(df.loc[i,'lng'],df.loc[i,'lat'])[1]

def point(arr):
    return Point(arr['wgs_lng'],arr['wgs_lat'])
df['geometry'] = df.apply(point,axis=1)

base = G_gdf[1].plot(figsize=(10,10),edgecolor='grey')
df.plot(column='category',legend=True,ax=base,markersize=8,cmap='Paired')
kuanzhai.plot(marker='*',ax=base,color='red',markersize=200)
plt.show()

这回看样子是对了，不过右上角这一坨是什么鬼，于是我打开了百度地图：

原来如此，这里开了一家百盛，所以说，这个POI数据有时候直接放进模型里，问题大大的，会造成很大的bias，不过国内土地利用的数据又是机密，只能用POI这个瘸腿的玩意，把实证论文写在祖国大地上实在是有点困难啊。

下面我们来画一下POI分布的热力图：

import seaborn as sns
meishi = df[df['category']=='美食']
base = G_gdf[1].plot(figsize=(10,10),edgecolor='blue',alpha=0.3)
sns.kdeplot(meishi['wgs_lng'],meishi['wgs_lat'],shade=True,shade_lowest=False,cmap='Greys',n_levels=20,alpha=0.8,legend=False)
meishi.plot(ax=base,color='red',markersize=5)
#kuanzhai.plot(marker='*',ax=base,color='red',markersize=200)
plt.xlim(west,east)
plt.ylim(south,north)
plt.title('美食POI空间分布核密度图',fontsize=20)
plt.show()

嘎嘎，美食POI核密度分布就画好了，还是比较符合实际的，颜色最深的那个地方正好是奎星楼街那片，的确是美食聚集地，什么冒椒火辣，糖油果子，泡椒鸡爪子，冰粉，想的我都流口水了，赶紧开学吧，我要回成都

至此本篇文章就结束啦，感谢大家的支持，最近更新的几篇文章都比较基础，希望大家都能用得到，下次可能会更新一些更有挑战的内容，希望大家多多转发多多点赞，大家的支持是我坚持下去的动力（如果要是能打赏就更好了），谢谢各位老板！

往期内容: