使用Python實現子區域數據分類統計

目錄

1. 前言
2. geopandas簡介
3. 子區域數據分類統計
4. 總結

1、前言

       最近碰到一個需求，須要統計某省內的全部市的某數據分佈狀況信息。現有該省的數據分佈狀況以及該省的行政區劃數據。我經過geopandas庫實現了這一需求，在這裏簡單記錄之，供須要的人借鑑。

2、geopandas簡介

想必你們對pandas都不陌生，它是一個開源的強大的Python數據分析工具。pandas確實作到了靈活、快速、高效的進行數據處理，而geopandas是在pandas的基礎上添加了對空間數據的支持，實現了讀取空間數據以及對空間數據進行處理。關於其介紹和安裝等請參考其github主頁，本文再也不贅述。

3、子區域數據分類統計

       直接進入正題，現有某省的分類統計數據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字典當中。

4、總結

       這樣就能夠實現對該省的分類統計數據進行進一步細分，取出每一個市的數據分類信息。固然並必定侷限於省和市，好比全球和國家或者國家和省等。只要存在包含關係便可經過此種方式進行處理。這是雞年的第一篇博客，願全部人今年都能有個好的結果！

---

打個廣告，年初我與朋友成立了一家公司，現長期招聘，以下：

公司名稱：武漢一格空間科技有限公司
官網主頁：http://www.phitrellis.com
公司簡介：咱們是一家初創公司，專一於地理信息系統，目前方向爲研發一套地信大數據管理、可視化、分析平臺
工做地址：甘肅蘭州
招聘崗位：前端、後端都可
職位要求：沒有特殊要求，只要你有能力，有拼搏進取的精神，能吃苦、肯學、肯上進，願意與公司共同成長，咱們都很歡迎
技術框架：前端主要用到React，後端主要用到Python、Scala，包含數據處理（numpy、pandas）、地信基礎分析（geopandas、rasterio等）、Hadoop\Spark分佈式框架。
薪資待遇：薪資面議，只要你足夠優秀，咱們願意提供相應的薪資，公司各項福利優厚，工做環境寬鬆，不強制996，我的愛好和興趣驅動。
聯繫方式：能夠直接經過博客園聯繫，或Email：shoufengwei@phitrellis.com