以最常见的3通道影像为例,假设要合并2幅3通道遥感影像为一幅6通道影像,并且不损失地理信息,就必须要请出我们的主角"GDAL"。
要使用GDAL'的功能,首先要安装GDAL,具体的安装教程在GDAL安装教程。
安装好了之后就可以简答使用了,比如gdal的读写功能,裁剪拼接功能,,,可以自己小试牛刀一下!
python GDAL 影像处理系列之裁剪与合并影像
python GDAL 影像处理系列之给遥感影像加地理信息
二话不说,直接上代码:
import gdal
import os
import numpy as np
#读取tif文件函数
def readTif(fileName):dataset = gdal.Open(fileName)if dataset == None:print(fileName+"文件无法打开")returnim_width = dataset.RasterXSize #栅格矩阵的列数im_height = dataset.RasterYSize #栅格矩阵的行数im_bands = dataset.RasterCount #波段数im_data = dataset.ReadAsArray(0,0,im_width,im_height)#获取数据im_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息im_proj = dataset.GetProjection()#获取投影信息im_blueBand = im_data[0,0:im_height,0:im_width]#获取蓝波段im_greenBand = im_data[1,0:im_height,0:im_width]#获取绿波段im_redBand = im_data[2,0:im_height,0:im_width]#获取红波段#im_nirBand = im_data[3,0:im_height,0:im_width]#获取近红外波段im_dtype = im_data.dtype.namereturn im_data, im_width,im_height, im_bands, im_geotrans ,im_proj,im_dtype, im_blueBand, im_greenBand, im_redBand#保存tif文件函数def writeTiff(im_data,im_width,im_height,im_bands,im_geotrans,im_proj,path):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif len(im_data.shape) == 2:im_data = np.array([im_data])else:im_bands, (im_height, im_width) = 1,im_data.shape#创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(path, im_width, im_height, im_bands, datatype)if(dataset!= None):dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数dataset.SetProjection(im_proj) #写入投影for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del datasetdef get_file_names(data_dir, file_type = ['tif','tiff']):result_dir = [] result_name = []for maindir, subdir, file_name_list in os.walk(data_dir):for filename in file_name_list:apath = maindir+'/'+filenameext = apath.split('.')[-1] if ext in file_type:result_dir.append(apath)result_name.append(filename)else:passreturn result_dir, result_namein_dir1 = 'F:/project/cut256_1'
in_dir2 = 'F:/project/cut256_2'
out_dir = 'F:/project/combine256_6bands'file_type = 'tif'
data_dir_list1,_ = get_file_names(in_dir1, file_type)
data_dir_list2,_ = get_file_names(in_dir2, file_type)
#data_dir_list = data_dir_list1 + data_dir_list2for each_index, each_dir in enumerate(data_dir_list1):img1, width1, height1, bands1, geotrans1, proj1,dtype1, blueband1, greenband1, redband1 = readTif(each_dir)img2, width2, height2, bands2, geotrans2, proj2,dtype2, blueband2, greenband2, redband2 = readTif(data_dir_list2[each_index])#print(each_dir)#print(dtype1)if 'int8' in dtype1:datatype = gdal.GDT_Byte elif 'int16' in dtype1:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32driver = gdal.GetDriverByName("GTiff")# print(type(driver))each_out_dir = out_dir + '/cut256_201903_201906_' + each_dir.split('_')[-4] + '_' + each_dir.split('_')[-3] + '_' + each_dir.split('_')[-2] + '_' + each_dir.split('_')[-1]#each_out_dir = 'C:/Users/Dell/Desktop/guigang/trans_6bands.tif'#datatype = 'uint8'print('each_out_dir: ', each_out_dir)new_dataset = driver.Create(each_out_dir, width1, height1, bands1+bands2, gdal.GDT_Byte)print(type(new_dataset))#print(each_out_dir)new_dataset.SetGeoTransform(geotrans1)new_dataset.SetProjection(proj1)new_dataset.GetRasterBand(1).WriteArray(blueband1)new_dataset.GetRasterBand(2).WriteArray(greenband1)new_dataset.GetRasterBand(3).WriteArray(redband1)new_dataset.GetRasterBand(4).WriteArray(blueband2)new_dataset.GetRasterBand(5).WriteArray(greenband2)new_dataset.GetRasterBand(6).WriteArray(redband2)new_dataset = Noneprint('combine over')
假如你的数据不是3通道的,直接在读取完图像之后获取每个波段的数据,然后在合并的时候,就可以直接在“new_dataset.GetRasterBand(6).WriteArray(redband2)”后面添加你自己的数据格式就ok了。也可以不像我这样rgbr'g'b'排列,你也可以选择rr'gg'bb'等格式,这个看你需求了!
在代码调试过程中问题,在下方留言就ok