使用Python实现子区域数据分类统计
目录
一、前言
最近碰到一个需求，需要统计某省内的所有市的某数据分布情况信息。现有该省的数据分布情况以及该省的行政区划数据。我通过geopandas库实现了这一需求，在这里简单记录之，供需要的人借鉴。
二、geopandas简介
想必大家对pandas都不陌生，它是一个开源的强大的Python数据分析工具。pandas确实做到了灵活、快速、高效的进行数据处理，而geopandas是在pandas的基础上添加了对空间数据的支持，实现了读取空间数据以及对空间数据进行处理。关于其介绍和安装等请参考其github主页，本文不再赘述。
三、子区域数据分类统计
直接进入正题，现有某省的分类统计数据shp文件以及此省的行政区划数据shp文件。
3.1 引入geopandas
为了使用geopandas库，首先需要将其引入。代码如下：
from geopandas import *
3.2 读取此省分类统计数据及行政区划数据
然后从该省的分类统计数据shp文件中读出此数据。代码如下：
provincedata = GeoDataFrame.from_file(provinceshpfile)
代码很简单，只要给GeoDataFrame.from_file函数传入shp文件路径即可，其得到的是一个GeoDataFrame对象，类似于pandas中的DataFrame，区别会在下文讲到。同理，行政区划数据通过以下代码读入：
regiondata = GeoDataFrame.from_file(regionshppath)
3.3 投影转换
要使上述两个数据能够进行处理必须先要将其转换成相同的投影。geopandas进行投影转换很简单。代码如下：
gpd_new=gpd.to_crs(crs)
其中gpd表示原始数据，gpd_new为转换后的数据，crs表示需要转换的投影参数，在https://github.com/geopandas/geopandas/blob/46e50fe5a772944b325bc3c8806f4f96da76a0d8/doc/source/projections.rst中对此参数进行了详细描述，大家可以参照。
所以我们只需要将上述两种数据转换到同一投影即可，问题是假设我们不知道它们的投影类型，那么也很容易，只需要将其都转换成4326或其他投影即可，这样就能保证二者转换后为同一投影。但是问题又来了，如果该省分类数据特别大将会导致投影转换耗时过长。其实此处有个简单方法，我们只需要读出分类数据的crs并将行政区划数据转换成此投影即可，这样不仅代码简单而且能够节省大量时间。代码如下：
regiondata_crs = regiondata.to_crs(provincedata.crs)
3.4 对该省逐市进行数据分类提取
现在二者投影已经相同，我们不得不面对最核心的问题，如何能够从省的分类数据中提取出逐市分类数据情况。我们可以想到必须要将每个市的空间数据与该省的所有分类数据进行相交判断，判断哪些分类数据与该市相交，然后完成相应处理。
3.4.1 获取市的空间数据
首先对行政区划数据进行循环得到每个市的空间数据。代码如下：
for i in range(0, len(regiondata_crs)):
geo = regiondata_crs.geometry[i]
name = regiondata_crs.NAME[i]
其中geo就是取到的当前市的空间数据，可以看出GeoDataFrame与DataFrame的区别就在于多了一个geometry字段，它包含了数据的空间信息，可以对该字段进行空间操作。假设该shp文件还包含了一个NAME属性，那么我们就可以用“.NAME”的方式提取出当前市的name数据，其他属性同理。
3.4.2 获取此空间内的分类数据信息
有了geo之后就可以将其与省的分类数据中的每一个对象进行相交判断（循环判断），根据结果进行相应处理。这里我们假设统计不同种类数据的面积值，即每种类型的数据在该市所占面积大小。代码如下：
area = {}
for i in range(0, len(provincedata)):
geo_province = provincedata.geometry[i]
id = str(provincedata.ID[i])
if geo_province.intersects(geo) and geo_province.is_valid :
temp_geo = geo_province.intersection(geo)
temp_area = temp_geo.area
if area.get(id) == None:
area[id] = temp_area
else :
area[id] = area.get(id) + temp_area
其中area为字典对象，用来存储不同数据在该市所占面积。
首先通过provincedata.geometry[i]获取该数据的空间信息geo_province，然后使用provincedata.ID[i]取出该数据的编号值id。
geo_province.intersects(geo)用来判断该数据与当前市（geo）在空间上是否相交，geo_province.is_valid用来判断该数据是否合法，有无自相交环等。
如果相交则进行处理，首先通过geo_province.intersection(geo)来获取相交的部分temp_geo，然后通过temp_geo.area获取相交部分的面积temp_area。可以看出在geopandas中只需要对geometry对象使用area属性即可获取其面积。
最后将面积以id为key保存到area字典当中。
四、总结
这样就可以实现对该省的分类统计数据进行进一步细分，取出每个市的数据分类信息。当然并一定局限于省和市，比如全球和国家或者国家和省等。只要存在包含关系即可通过此种方式进行处理。这是鸡年的第一篇博客，愿所有人今年都能有个好的结果！
打个广告，年初我与朋友成立了一家公司，现长期招聘，如下：
公司名称：武汉一格空间科技有限公司
官网主页：http://www.phitrellis.com
公司简介：我们是一家初创公司，专注于地理信息系统，目前方向为研发一套地信大数据管理、可视化、分析平台
工作地址：甘肃兰州
招聘岗位：前端、后端均可
职位要求：没有特殊要求，只要你有能力，有拼搏进取的精神，能吃苦、肯学、肯上进，愿意与公司共同成长，我们都很欢迎
技术框架：前端主要用到React，后端主要用到Python、Scala，包含数据处理（numpy、pandas）、地信基础分析（geopandas、rasterio等）、Hadoop\Spark分布式框架。
薪资待遇：薪资面议，只要你足够优秀，我们愿意提供相应的薪资，公司各项福利优厚，工作环境宽松，不强制996，个人爱好和兴趣驱动。
联系方式：可以直接通过博客园联系，或Email：shoufengwei@phitrellis.com
