利用Python實現矢量逐個圖斑裁剪柵格，形成圖斑對應的柵格文件

平時工作中存在，利用矢量裁剪柵格的要求，但多數情況下基於完整的單個矢量裁剪柵格，非利用矢量中某個圖斑裁剪柵格，因此做以下工作。
1、將矢量按照單個圖斑要素拆分成shp
這裡用了縣區的矢量。代碼如下：

from osgeo import gdal
import osgeo.ogr as ogr
input_shape = r"C:/分類/縣區投影.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(input_shape, 1)
layer = dataSource.GetLayer()
print('the length of layer:', len(layer))
for i, feature in enumerate(layer):
# 新建DataSource，Layer
fid = feature.GetField('DISTNAME')#讀取當前Feature某一字段的屬性值
out_ds = driver.CreateDataSource(fid+".shp")
out_lyr = out_ds.CreateLayer(fid+".shp", layer.GetSpatialRef(), ogr.wkbPolygon)
def_feature = out_lyr.GetLayerDefn()
# 生成Shapefile文件
# current_union = layer[0].Clone()
geometry = feature.GetGeometryRef()
current_union = geometry.Clone()
current_union = current_union.Union(geometry).Clone()
out_feature = ogr.Feature(def_feature)
out_feature.SetGeometry(current_union)
out_lyr.ResetReading()
out_lyr.CreateFeature(out_feature)

效果如下，原始為一個完整的圖層，拆分後為單獨的多個圖層。
原矢量：

拆分後矢量：

2、利用拆分後的矢量循環裁剪柵格
這裡使用的掩膜提取的方法進行裁剪，裁剪不規則的范圍，非外接矩形。
代碼如下：

from osgeo import gdal
import osgeo.ogr as ogr
# tif輸入路徑，打開文件
input_raster = r"C:/DEM/坡度.tif"
# 柵格文件路徑，打開柵格文件
input_raster=gdal.Open(input_raster)
#匹配文件名字，可以編寫讀取文件夾文件來替換
name =['開縣',..........., '石柱土家族自治縣']
for n in name:
#開始裁剪，一行代碼
ds = gdal.Warp(n+".tif",#生成的柵格
input_raster,
format = 'GTiff',
#矢量文件
cutlineDSName = n+".shp",
#cutlineWhere="FIELD = 'whatever'",
dstNodata = 0)

原柵格：

拆分後柵格：

3、逐柵格統計信息
根據個人需要逐個統計相應信息，這裡做的是統計坡譜信息熵。