目錄
一、前言
最近碰到一個需求，需要統計某省內的所有市的某數據分布情況信息。現有該省的數據分布情況以及該省的行政區划數據。我通過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