剪切tif图像

图像"/>

剪切tif图像

import gdal
import numpy as np
import osdef clip(int_path, out_path, use_band, size):inList = [name for name in os.listdir(int_path)]for file in inList:file = file.lower()if file.endswith('.tif'):print("待裁剪影像：", file)input_path = os.path.join(int_path, file)  # 文件的完整路径output_path = os.path.join(out_path, file.strip(".tif"))  # 文件夹的完整路径if not os.path.exists(output_path):os.makedirs(output_path)in_img = gdal.Open(input_path)  # 读取要切的原图if in_img is None:print("打开失败！")else:print("打开成功！")width = in_img.RasterXSize  # 获取数据宽度height = in_img.RasterYSize  # 获取数据高度bandsize = in_img.RasterCount  # 获取数据波段数im_geotrans = in_img.GetGeoTransform()  # 获取仿射矩阵信息im_proj = in_img.GetProjection()  # 获取投影信息data_array = in_img.ReadAsArray()if 'int8' in data_array.dtype.name:datatype = gdal.GDT_Int16elif 'int16' in data_array.dtype.name:datatype = gdal.GDT_Int16else:datatype = gdal.GDT_Float32print("==================", file, "影像信息======================")print("width:", width, "height:", height, "outbandsize:", bandsize, "datatype",datatype)if bandsize < use_band:print("影像波段数小于所使用波段数！")else:int_data_all = []  # 读取原图中所需的前ues_band个波段，并将前use_band个波段读入“data_all”for inband in range(use_band):int_band = in_img.GetRasterBand(inband + 1)   # 1 2 3 4int_data_all.append(int_band)print('总共打开，并读取到内存的影像波段数目：{0}'.format(len(int_data_all)))# 定义切图的大小（矩形框）if size > width or size > height:print("裁剪尺寸大于原始影像，请重新确定输入！")else:col_num = int(width / size)  # 宽度可以分成几块  # 定义切图的起始点坐标row_num = int(height / size)  # 高度可以分成几块if (width % size != 0):col_num += 1if (height % size != 0):row_num += 1num = 0  # 记录一共有多少块print("row_num:%d   col_num:%d" % (row_num, col_num))for i in range(row_num):  # 0, 1, 2, 3 ... 127offset_y = i * size   # 0, 128, 256, 384, ..., 16256for j in range(col_num):  # 0, 1, 2, 3, ..., 127offset_x = j * size   # 0, 128, 256, 384, ..., 16256# 矩形框的大小b_ysize = min(height - offset_y, size)b_xsize = min(width- offset_x, size)print("width:%d     height:%d    offset_x:%d    offset_y:%d     b_xsize:%d     b_ysize:%d" % (width,        height,     offset_x,      offset_y,       b_xsize,        b_ysize))print("\n")out_data_all = []for band in range(use_band): # 0 1 2 3out_data_band = int_data_all[band].ReadAsArray(offset_x, offset_y, b_xsize, b_ysize)# print(out_data_band)out_data_band[np.where(out_data_band < 0)] = 0out_data_band[np.isnan(out_data_band)] = 0print("min", np.min(out_data_band), "max", np.max(out_data_band))out_data_all.append(out_data_band)num=num+1print("out_data第{0}矩形已成功写入：{1}个波段".format(num, len(out_data_all)))# 获取Tif的驱动，为创建切出来的图文件做准备gtif_driver = gdal.GetDriverByName("GTiff")file = output_path + '\%04d.tif' % numprint("out_file", file)# 创建切出来的要存的文件out_img = gtif_driver.Create(file, size, size, use_band, datatype)print("create new tif file succeed")# 获取原图的原点坐标信息ori_transform = in_img.GetGeoTransform()# 读取原图仿射变换参数值top_left_x = ori_transform[0]  # 左上角x坐标w_e_pixel_resolution = ori_transform[1]  # 东西方向像素分辨率top_left_y = ori_transform[3]  # 左上角y坐标n_s_pixel_resolution = ori_transform[5]  # 南北方向像素分辨率# 根据反射变换参数计算新图的原点坐标top_left_x = top_left_x + offset_x * w_e_pixel_resolutiontop_left_y = top_left_y + offset_y * n_s_pixel_resolutionprint("top_left_x", top_left_x, "top_left_y", top_left_y)# 将计算后的值组装为一个元组，以方便设置dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4],ori_transform[5])# 设置裁剪出来图的原点坐标out_img.SetGeoTransform(dst_transform)# 设置SRS属性（投影信息）out_img.SetProjection(in_img.GetProjection())# 写入目标文件for w_band in range(use_band):  # 0 1 2 3out_img.GetRasterBand(w_band + 1).WriteArray(out_data_all[w_band])# 将缓存写入磁盘out_img.FlushCache()print("=============已成功写入第{}个矩形===============".format(num))del out_imgdel in_imgimg_filepath = "F:\WV3\wv3\MUL"
out_filepath = "F:\WV3\crop\MUL"
band_num = 4
rectangle_size = 128
clip(img_filepath, out_filepath, int(band_num), int(rectangle_size))