利用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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

效果如下，原始为一个完整的图层，拆分后为单独的多个图层。
原矢量：

拆分后矢量：

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

原栅格：

拆分后栅格：

3、逐栅格统计信息
根据个人需要逐个统计相应信息，这里做的是统计坡谱信息熵。