之前分享过一篇分幅裁剪的文章:【Python&RS】基于GDAL遥感影像分幅裁剪,只不过这篇文章当时编写的逻辑是自己输入需要裁剪多少行多少列,由于大家可能并没有直观地希望自己裁剪多少行列,所以非常局限。今天跟大家分享一下使用固定尺寸对遥感影像进行分幅裁剪,即每张裁剪的影像都是一样大的。
原创作者:RS迷途小书童
1. 代码逻辑
逻辑其实很简单,就是利用输入的尺寸计算出需要裁剪多少行多少列。然后按照之前文章的方法进行裁剪即可。当然这里要注意一下边缘处的处理。
2. 代码主函数
由于和之前发的文章原理差不多,所以我这里就不做过多的解释了,代码中都有详细的注释,如果大家有什么问题可以给我留言。
ds = gdal.Open(filepath)
ds_width = ds.RasterXSize
ds_height = ds.RasterYSize
ds_bands = ds.RasterCount
raw = int(ds_height / size) + 1
col = int(ds_width / size) + 1
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_frame
left_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)
ds_geo = ds.GetGeoTransform()
top_left_x = ds_geo[0]
top_left_y = ds_geo[3]
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) # 将无效值设为0
ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中
ds_result = None
del ds_result
3. 总结
标签:RS,Python,frame,裁剪,分幅,geo,ds,GDAL,left From: https://www.cnblogs.com/RSran/p/18117161今天主要分享的是遥感影像的分幅裁剪,大家可以用这段代码减少数据量,也可以用它制作样本集。同时结合之前发的那篇分幅裁剪的文章,基本上概括了所有的分幅裁剪的情况,友友们也可以根据自己的需求修改一下代码。如果大家在学习Python或者RS时有什么问题,可以随时留言交流!如果大家对批量处理有兴趣同样可以留言给博主,博主会分享相关代码以供学习!