gdal：矢量裁剪栅格数据

#!usr/bin/python
# -*- coding: utf-8 -*-
from osgeo import gdal
from osgeo import ogr
from osgeo.gdalconst import GDT_Byte, GDT_Float32, GDT_UInt16
import os
from tqdm import tqdm

shp_path = '/media/dell/CH3/fishnet/fishnet.shp'
img_path ='/media/dell/CH3/imges.tif'
save_path ='/media/dell/CH3/savepath/'
#field_paths = os.chdir(field_path)

os.chdir("/media/dell/CH3/savepath/")
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open(shp_path, 0)
layer = ds.GetLayer()
numFeatures = layer.GetFeatureCount()

for i in tqdm(range(numFeatures)):
    feature = layer.GetNextFeature()
    id = feature.GetField('Name')
    save_name=save_path+id+'.tif'  
    select= "Name=" + "'"+ str(id)+ "'"
    result = gdal.Warp(
                save_name,
                img_path,
                format = 'GTiff',
                cutlineDSName = shp_path,
                cutlineWhere=select, # 用于裁剪的矢量
                cropToCutline = True, # 是否使用cutlineDSName的extent作为输出的界线
                dstNodata = 255 # 输出数据的nodata值
            )
    reclass=result.GetRasterBand(1).ReadAsArray().astype(int) #重分类过程
    '''
    '''   
    band=result.GetRasterBand(1)
    gtiff_driver =gdal.GetDriverByName('GTiff')
    result1 = gtiff_driver.Create(save_name,band.XSize, band.YSize, 1, GDT_Byte)
    result1.SetProjection(result.GetProjection())  
    result1.SetGeoTransform(result.GetGeoTransform())  
    result1.WriteArray(reclass) #也可直接result
    result1.FlushCache()

    del result, reclass
print("finish!!!!")