转载：“可视化”的数据分析落伍了？

作为一名程序员兼职业余空间数据分析师，空间数据的处理一直是一抹挥之不去的乌云。虽然GIS软件可视化的操作界面、包罗万象的功能，已经能解决日常工作中几乎所有问题，但对于身为程序员的我来说，一切不能用键盘上26个字母解决的问题，都是不科学的~所以这些年我一直致力于思考一个问题，如何把Arcgis的功能，用代码的方式实现。

我了解到，Geopandas是一个极好的工具，今天我们就以一个案例来管窥一二。

数据准备&任务明确

我有一张上海的街镇地图，一份上海的房源数据和一份上海餐厅数据。

今天，我们需要制作每一个房源的1KM缓冲区，计算缓冲区内有多少家餐厅，并且计算上海每个街镇内有多少家餐厅，然后绘制在地图上展示。

这是一个理论上简单但操作繁琐的过程，有些同学可能已经知道GIS的操作：

导入地图数据——导入两份EXCEL——并关联经纬度——获得点数据——并保存点数据图层——然后对房源点数据进行缓冲区制作——保存缓冲区图层——分别使用空间关联把餐厅数据关联到上海街镇和房源缓冲区图层中——并分别保存图层……

好了好了，我已经说不下去了，我们来看看代码怎么解决吧~

导入数据

先导入各种包

import pandas as pdimport geopandas as gpdfrom shapely import geosfrom shapely.geometry import Pointimport fionaimport matplotlib.pyplot as pltfrom fiona.crs import from_epsg,from_string

我们来观察下python中用到的地理处理工具包的名字：Geopandas. 顾名思义，这个包可以让我们像使用pandas一样处理地理数据，大家可以想象一下这过程，一定如丝般顺滑！Geopandas其实是各种地理数据分析包的集大成者，包括shapely,Fiona等,当然还有数据分析相关的包numpy,pandas，所以，在这些提及的包中的功能，都可以混合使用，切换不留痕迹。这里就一并导入了。

导入地图shp

shanghai_map = gpd.GeoDataFrame.from_file('./上海街镇/上海街镇.shp') #读取shapfile数据为geodataframe格式

Geopandas提供了一种数据格式叫GeoDataFrame，用直白的话概括就是DataFrame加了一列数据，表达地理信息。导入功能Geopandas底层调用的是Fiona包，所以，一些基本参数和可以导入的数据格式，可以参考Fiona的说明文档。

查看导入后的数据

 

 

和DataFrame如出一辙，但多了一列Geometry来存放地理信息，绘图看一下：

 

嗯！是长这样！

接下来，导入房源和餐厅数据

 

*这两份数据是CSV格式，导入成DataFrame，我们发现数据中含有经纬度字段，我们可以根据这两个字段，也把数据转换成GeoDataFrame格式

 

通过经纬度转换点数据：

 

defpoint_to_geo(df,lon,lat):df['geometry'] = gpd.GeoSeries(list(zip(df[lon],df[lat]))).apply(Point) #识别经纬度，转换点数据 df = gpd.GeoDataFrame(df) #转换Geodataframe格式 df.crs = {'init':'epsg:4326'} #定义坐标系WGS84del df[lon]del df[lat]return dfhouse_data = point_to_geo(house_data,'lon_WGS','lat_WGS') #转换Geodataframe格式restaurant_data = point_to_geo(restaurant_data,'lon_WGS','lat_WGS') #转换Geodataframe格式

*官方文档请看这里：

 

绘图，看看成功与否——以上海街镇为底图，两份点数据叠加在底图上

base = shanghai_map.plot(color='lightyellow',edgecolor='black',figsize=(15, 15)) #画底图restaurant_data.plot(ax=base,marker='o', color='green', markersize=5) #在底图上叠加餐厅点数据house_data.plot(ax=base,marker='o', color='red', markersize=5) #在底图上叠加房源点数据plt.gca().xaxis.set_major_locator(plt.NullLocator())#去掉x轴刻度plt.gca().yaxis.set_major_locator(plt.NullLocator())#去掉y轴刻度plt.savefig('./map.png',dpi=400) #保存图片

*这里设置下图形的基本参数，画图功能是基于matplotlib，所以它的一些写法和功能是通用的。

显示OK，那么就进入核心环节~

分析过程

在制作缓冲区和空间关联之前，先需要对图形进行投影变换，这里仍然写一个函数：

def wgs84_to_CGCS2000(df,code):result = df.to_crs(from_epsg(code)) return resultshanghai_map_pcs = wgs84_to_CGCS2000(shanghai_map,4549)house_data_pcs = wgs84_to_CGCS2000(house_data,4549)restaurant_data_pcs = wgs84_to_CGCS2000(restaurant_data,4549)

*只要知道投影坐标系的ESPG代码，就可以任意转换，我们使用GSC2000坐标系。

接着，构建缓冲区图层，使用buffer()方法，参数为缓冲区的半径。

house_data_buffer = house_data_pcs.buffer(1000) #建立一公里缓冲区base = shanghai_map_pcs.plot(color='lightyellow',edgecolor='black',figsize=(15, 15))house_data_buffer.plot(ax=base,color='gray', markersize=5,alpha=0.5)plt.gca().xaxis.set_major_locator(plt.NullLocator()) #去掉x轴刻度plt.gca().yaxis.set_major_locator(plt.NullLocator()) #去年y轴刻度plt.savefig('./map2.png',dpi=400) #保存图片

构建完成，同样输出成地图查看效果：

增加了透明度，我们可以看到市中心的房源密度较高。生成的数据为GeoSeries，数据为缓冲区的地理信息。这里需要使用任意方法把地理信息和其他字段匹配在一起。

buffer_temp = house_data_pcs[['name','geometry']]buffer_temp['geometry'] = house_data_bufferhouse_data_buffer = buffer_temp

下面，我们使用空间关联，连接餐厅点数据，并获取餐厅的个数。

spacial_join_restaurant = gpd.sjoin(house_data_buffer,restaurant_data_pcs,how='left',op='contains') #空间连接spacial_join_restaurant = spacial_join_restaurant.groupby(['name']).count()['title'].to_frame().reset_index() #聚合计算个数spacial_join_restaurant.columns = ['name','restaurant_count'] #更改列名，方便操作buffer_result = pd.merge(house_data_pcs,spacial_join_restaurant,left_on='name',right_on='name',how='left') #字段匹配

*方法是sjoin()，功能和GIS相同，分相交，包含和被包含三种。我们使用包含进行连接，并使用pandas的groupby()方法分组计数。

 

这样我们就完成了空间关联和计算，我们可以使用地图显示，并添加图例：

base = shanghai_result.plot(column='restaurant_count', cmap='Oranges',scheme = 'fisher_jenks',legend=True,edgecolor='black',figsize=(15, 15)) #按个数多少叠加底色plt.gca().xaxis.set_major_locator(plt.NullLocator()) #去掉x轴刻度plt.gca().yaxis.set_major_locator(plt.NullLocator()) #去年y轴刻度plt.savefig('./map3.png',dpi=400) #保存图片

这样就可以绘制出以街镇为单位的点密度图：

最后，我们需要把数据导出保存。我们可以保存为shapfile格式，这里推荐大家直接保存为csv格式，方便分享。也可以根据需要，转成任意格式。

shanghai_result = wgs84_to_CGCS2000(shanghai_result,4326) #转换成地理坐标系shanghai_result.columns = ['town','region','geometry','restaurant_count'] #更改列名为英文，因为数据保存对中文支持不佳buffer_result = wgs84_to_CGCS2000(buffer_result,4326) #转换成地理坐标系shanghai_result.to_csv('./result/shanghai_result.csv') #保存成csvshanghai_result.to_file('./result/shanghai_result.shp') #保存成shapfilebuffer_result.to_csv('./result/buffer_result.csv') #保存成csvbuffer_result.to_file('./result/buffer_result.shp') #保存成shapfile

结果检查

我们把保存的数据导入Arcgis，可以正常使用~太棒了！

是不是很实用呢？所以，对于空间数据分析中一些重复操作，工程化的作业，建议大家可以考虑用代码实现，高效，便捷。

不用代码的工具有吗？

有！亿信华辰的酷屏！酷屏剔除了繁琐的更新步骤，提供丰富的二次开发接口，支持自行改进、增加组件，无需换包或打补丁，不重启服务器的情况下即可完成组件的更新，为实施人员提供便利。

快去使用吧！