本文整體思路：在Python中使用Geopandas庫，依次讀取shp文件的每一個面狀要素，獲取其空間邊界信息並裁剪對應的柵格影像，計算所裁剪影像Value值的眾數，將其設置為對應面狀要素的NewTYPE值，所有要素屬性值都改好之后保存為新的shp文件。
使用Python處理空間數據確實用的不多，所以一個星期以來一直深受這個程序的折磨，官方文檔、博客、谷歌、百度、論文，能用的方法都給用了，但是進度還是很慢，特別是當看到這篇博客的時候。。。好氣啊。。
不過幸虧頭比較鐵，雖敗不餒，慢慢一步一步調試找問題，一個一個解決，終於啃下了這根硬骨頭。（對寫程序來說，調試真的是第一法寶啊！）
好吧進入正題。。。
Pandas是一款高性能的Python數據分析庫，而Geopandas是由Shapely、Fiona、PyProj、matplotlib以及其他必須的庫共同構建的Pandas地理空間擴展，它既可以處理空間數據、也可以處理屬性數據。看到有些博客說在讀取shp文件的時候用Geopandas庫，而在編輯、導出的時候用pyshp比較好，其實不然，Geopandas也提供了功能完備的導出接口，而且使用特別方便，只不過需要注意一個小的細節問題，否則就會報錯。Rasterio是基於GDAL庫二次封裝的主要用於空間柵格數據處理的Python庫，本程序需要對柵格影像進行裁剪，因此也需要引用這個庫。兩個庫的官方文檔如下，參考的時候要注意版本問題，不同的版本有些接口可能已經改變。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后可以看到效果如下：
至此全部完成！