之前分享过一篇分幅裁剪的文章：【Python&RS】基于GDAL遥感影像分幅裁剪，只不过这篇文章当时编写的逻辑是自己输入需要裁剪多少行多少列，由于大家可能并没有直观地希望自己裁剪多少行列，所以非常局限。今天跟大家分享一下使用固定尺寸对遥感影像进行分幅裁剪，即每张裁剪的影像都是一样大的。
原创作者：RS迷途小书童
博客地址：https://blog.csdn.net/m0_56729804?type=blog
1. 代码逻辑
逻辑其实很简单，就是利用输入的尺寸计算出需要裁剪多少行多少列。然后按照之前文章的方法进行裁剪即可。当然这里要注意一下边缘处的处理。
2. 代码主函数
由于和之前发的文章原理差不多，所以我这里就不做过多的解释了，代码中都有详细的注释，如果大家有什么问题可以给我留言。
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
raw = int(ds_height / size) + 1
col = int(ds_width / size) + 1
print("分割后影像行数：", raw)
print("分割后影像列数：", col)
print("分割后影像数为：", raw * col)
for j in range(0, raw):print("正在裁剪第%s行......" % (j + 1))for k in range(0, col):raw_frame = size# 计算每幅图像的高度col_frame = size# 计算每幅图像的宽度left_x = j * raw_frameleft_y = k * col_frame# 计算当前影像的左上角像素坐标raw_frame = min(ds_height-left_x, raw_frame)col_frame = min(ds_width-left_y, col_frame)# 防止幅宽超过整幅影像driver = gdal.GetDriverByName('GTiff') # 载入数据驱动，用于存储内存中的数组ds_result = driver.Create(out_path+"%s_%s.tif" % (j+1, k+1), col_frame, raw_frame,bands=ds_bands, eType=gdal.GDT_Byte) # 创建空tifds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数top_left_x = ds_geo[0] # 原始影像左上角x的投影坐标top_left_y = ds_geo[3] # 原始影像左上角y投影坐标top_left_x = top_left_x + left_y * ds_geo[1]top_left_y = top_left_y + left_x * ds_geo[5]# 计算得到当前影像的左上角投影坐标ds_geo = (top_left_x, ds_geo[1], ds_geo[2], top_left_y, ds_geo[4], ds_geo[5])# 新影像的仿射地理变换参数ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds.GetProjection()) # 导入投影信息array_band = []for i in range(1, ds_bands+1):array_band = ds.GetRasterBand(i).ReadAsArray(left_y, left_x, col_frame, raw_frame).astype(np.float64)# 根据左上角的像素坐标和幅宽读取指定区域内的数据ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中ds_result = Nonedel ds_result3. 总结
今天主要分享的是遥感影像的分幅裁剪，大家可以用这段代码减少数据量，也可以用它制作样本集。同时结合之前发的那篇分幅裁剪的文章，基本上概括了所有的分幅裁剪的情况，友友们也可以根据自己的需求修改一下代码。如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！
)
)
