Geopandas庫——面狀向量裁剪柵格影像，並依據Value值更改向量屬性

摘要： 本文整體思路：在Python中使用Geopandas庫，依次讀取shp檔案的每一個面狀要素，獲取其空間邊界資訊並裁剪對應的柵格影像，計算所裁剪影像Value值的眾數，將其設定為對應面狀要素的NewTYPE值，所有要素屬性值都改好之後儲存為新的shp檔案。 使用Python處理空間資料確實用...

本文整體思路：在Python中使用Geopandas庫，依次讀取shp檔案的每一個面狀要素，獲取其空間邊界資訊並裁剪對應的柵格影像，計算所裁剪影像Value值的眾數，將其設定為對應面狀要素的NewTYPE值，所有要素屬性值都改好之後儲存為新的shp檔案。

使用Python處理空間資料確實用的不多，所以一個星期以來一直深受這個程式的折磨，官方文件、部落格、谷歌、百度、論文，能用的方法都給用了，但是進度還是很慢，特別是當看到這篇部落格的時候。。。好氣啊。。

不過幸虧頭比較鐵，雖敗不餒，慢慢一步一步除錯找問題，一個一個解決，終於啃下了這根硬骨頭。（對寫程式來說，除錯真的是第一法寶啊！）

好吧進入正題。。。

Pandas是一款高效能的Python資料分析庫，而Geopandas是由Shapely、Fiona、PyProj、matplotlib以及其他必須的庫共同構建的Pandas地理空間擴充套件，它既可以處理空間資料、也可以處理屬性資料。看到有些部落格說在讀取shp檔案的時候用Geopandas庫，而在編輯、匯出的時候用pyshp比較好，其實不然，Geopandas也提供了功能完備的匯出介面，而且使用特別方便，只不過需要注意一個小的細節問題，否則就會報錯。Rasterio是基於GDAL庫二次封裝的主要用於空間柵格資料處理的Python庫，本程式需要對柵格影像進行裁剪，因此也需要引用這個庫。兩個庫的官方文件如下，參考的時候要注意版本問題，不同的版本有些介面可能已經改變。 ofollow,noindex">Geopandas參考文件 ； Rasterio參考文件 。

我是Windows 10系統，在Python中安裝Geopandas庫比較麻煩，不能用pip命令直接安裝，而需要先下載Anaconda再用conda命令安裝，這部分網上有很多參考資料，就不多贅述了。但是安裝完成之後發現它的一些依賴庫不能使用，需要pip命令將其解除安裝之後，再通過此地址： python依賴庫 下載對應的依賴庫並安裝就可以使用了。

本程式完整的程式碼如下：

# -*- coding: utf-8 -*-
from geopandas import *;
import rasterio as rio;
import rasterio.mask;

# 讀入向量和柵格檔案
shpdatafile='D:/PythonConda/Data/ShpData.shp'
rasterfile='D:/PythonConda/Data/Raster.tif'
out_file='D:/PythonConda/Data/OutShpData'
shpdata=GeoDataFrame.from_file(shpdatafile)
rasterdata=rio.open(rasterfile)

out_shpdata = shpdata.copy()
#投影變換，使向量資料與柵格資料投影引數一致
shpdata=shpdata.to_crs(rasterdata.crs)

for i in range(0, len(shpdata)):
# 獲取向量資料的features
geo = shpdata.geometry[i]
feature = [geo.__geo_interface__]
#通過feature裁剪柵格影像
out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, crop=True, nodata=rasterdata.nodata)
#獲取影像Value值，並轉化為list
out_list = out_image.data.tolist()
#除去list中的Nodata值
out_list = out_list[0]
out_data = []
for k in range(len(out_list)):
for j in range(len(out_list[k])):
if out_list[k][j] >=0:
out_data.append(out_list[k][j])
#求資料中的眾數
if len(out_data):
counts = np.bincount(out_data)
new_type = np.argmax(counts)
else:
new_type = None
#依據眾數值更改feature的NewTYPE屬性值
out_shpdata.NewTYPE[i] = new_type

#將屬性更改過的GeodataFrame匯出為shp檔案
out_shpdata.to_file(out_file)

需要注意的問題：

1）兩個檔案需要在同一座標系下，需要將用於裁剪的shp資料進行投影變換，將其投影引數變為柵格資料的投影引數，由以下程式碼實現：

shpdata=shpdata.to_crs(rasterdata.crs)

其中crs表示資料的投影引數，to_crs為投影引數轉換函式。

2）裁剪函式rasterio.mask.mask的引數問題，需要傳入的向量資料為GeoJSON資料，因此讀入的每一個面狀feature都需要用__geo_interface__函式進行格式轉換，這一步可以通過除錯來檢視具體的資料格式是否正確；此函式有兩個返回值，第一個記錄裁剪柵格的資料值，第二個記錄其座標轉換資訊，一般用到第一個返回值比較多；裁剪得到的柵格形狀其實是一個矩形，與feature的外接矩形區域一致，只是位於feature外部的畫素值預設被設定為Nodata，當然這可以通過傳入的引數進行設定。

3）獲取到裁剪的柵格後，通過.data來獲取其Value值，但此時還不能直接用於統計分析，需要將資料轉化為List，函式如下：

out_list = out_image.data.tolist()

此時還需要除錯來進一步確定out_list的資料內容，發現out_list[0]才是我們真正能用到的資料值。

4）獲取眾數的時候需要清除Nodata值的影響，因此用for迴圈把out_list中的非Nodata資料再組成一個新的List，用numpy的自帶函式求其眾數。因為所有的編輯並不能對shpdata源資料進行改變，所以需要構建一個shpdata的copy即out_shpdata，將求得的眾數賦給out_shpdata的NewTYPE，編輯完成之後再將out_shpdata匯出為完整的shp檔案。

5）GeoDataFrame.to_file函式可以將out_shpdata直接匯出為shp檔案，僅需要傳入一個路徑引數即可，但需要注意由於shp檔案包含.shp、.shx、.dbf和.prj，因此路徑只能是一個資料夾，而不能具體到.shp。如下程式碼所示：

out_file='D:/PythonConda/Data/OutShpData'

最後將生成的資料在Arcmap中開啟，設定顯示的Labels後可以看到效果如下：

至此全部完成！