python 裁剪影像

# 栅格数据裁剪
# 添加包
import gdal
import numpy as np #导入numpy库
import math

filename = 'C:\\Users\\admin\\Desktop\\空间信息处理与系统开发资料\\第五周\\Lesson5影像数据\\裁剪影像\\test1.tif'
dataset = gdal.Open(filename)
# 一共5个波段，读取每个波段
band = np.ones((5,dataset.RasterYSize,dataset.RasterXSize)) #新建三维数组，且初始值为1 三个参数以此为页、行，列
for i in range(dataset.RasterCount):
    band[i,:,:]=dataset.GetRasterBand(i+1).ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)
   # print(band[i])

# 读取影像的行列数
#print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))

offset_x = 0
offset_y = 0
# 定义切图的大小（矩阵框）
block_xsize = dataset.RasterXSize * 0.5  # 列
block_ysize = dataset.RasterYSize * 0.5  # 行


out_band=np.ones((5,math.ceil(block_ysize),math.ceil(block_xsize))) 
for i in range(dataset.RasterCount):
    out_band[i,:,:] = band[i,offset_y:math.ceil(block_ysize),offset_x:math.ceil(block_xsize)]
    #print(out_band[i])

##创建数据集
gtif_driver = gdal.GetDriverByName('GTiff')
out_ds = gtif_driver.Create('C:\\Users\\admin\\Desktop\\空间信息处理与系统开发资料\\第五周\\Lesson5影像数据\\裁剪影像\\test1_裁剪后.tif',
                            math.ceil(block_xsize),math.ceil(block_ysize),5,gdal.GDT_Float32)  #列、行
# 定义仿射矩阵和投影信息
ori_transform = dataset.GetGeoTransform()
####读取原图仿射变换参数值
######左上角坐标
top_left_x = ori_transform[0]
######东西方向像素分辨率
w_e_pixel = ori_transform[1]
######左上角y坐标
top_left_y = ori_transform[3]
######南北方向像素分辨率
n_s_pixel = ori_transform[5]

top_left_x = top_left_x + offset_x * w_e_pixel
top_left_y = top_left_y + offset_y * n_s_pixel


##存到数组里
dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])
##设置仿射变化参数与投影信息
out_ds.SetGeoTransform(dst_transform)
out_ds.SetProjection(dataset.GetProjection())
#print(out_band.shape)
#print(out_ds.RasterYSize,out_ds.RasterXSize)

#将裁剪过的图像的5个波段像素值写入栅格数据集
for i in range(dataset.RasterCount):
    out_ds.GetRasterBand(i+1).WriteArray(out_band[i])

# 清空缓存，删除数据集
out_ds.FlushCache()
dataset.FlushCache()
del out_ds
del dataset

在ArcGIS中显示