python 影像拼接

import gdal
##拼接，不只是简单地把两张影像拼接，更重要的是地理位置配准
#打开、读取影像
filename1='C:\\Users\\admin\\Desktop\\空间信息处理与系统开发资料\\第五周\\Lesson5影像数据\\拼接影像\\test1.tif'
filename2='C:\\Users\\admin\\Desktop\\空间信息处理与系统开发资料\\第五周\\Lesson5影像数据\\拼接影像\\test2.tif'
ds1=gdal.Open(filename1)
ds2=gdal.Open(filename2)
print(ds1.GetProjection())
#读取行列数
#test1.tif
cols1=ds1.RasterXSize
rows1=ds1.RasterYSize
#test2.tif
cols2=ds2.RasterXSize
rows2=ds2.RasterYSize
##读取仿射变换参数
transform1=ds1.GetGeoTransform()
transform2=ds2.GetGeoTransform()
##根据test1空间坐标信息，并根据像素值大小计算右下角坐标
#提取影像左上角坐标
minX1=transform1[0]
maxY1=transform1[3]
#提取水平、垂直分辨率
w_e_pixel1=transform1[1]
n_s_pixel1=transform2[5]
#计算右下角坐标
maxX1=minX1+cols1*w_e_pixel1
minY1=maxY1+rows1*n_s_pixel1
##test2
minX2=transform2[0]
maxY2=transform2[3]
w_e_pixel2=transform2[1]
n_s_pixel2=transform2[5]

maxX2=minX2+cols2*w_e_pixel2
minY2=maxY2+rows2*n_s_pixel2
#设置拼接影像的坐标
minX=min(minX1,minX2)
maxX=max(maxX1,maxX2)
minY=min(minY1,minY2)
maxY=max(maxY1,maxY2)
#计算拼接影像的行与列
cols=int((maxX-minX)/w_e_pixel1)
rows=int((maxY-minY)/abs(n_s_pixel1))
#创建拼接影像的空间信息（仿射变换参数）
geotransform=[minX,w_e_pixel1,0,maxY,0,n_s_pixel1]
#创建新的栅格数据集，并定义放射矩阵和投影信息
filename_out='C:\\Users\\admin\\Desktop\\空间信息处理与系统开发资料\\第五周\\Lesson5影像数据\\拼接影像\\pinjie.tif'
driver=gdal.GetDriverByName('GTiff')
dataset_out=driver.Create(filename_out,cols,rows,5,gdal.GDT_Float32)
dataset_out.SetProjection(ds1.GetProjection())
dataset_out.SetGeoTransform(geotransform)
#计算test1在拼接图像左上角的偏移值
xOffset1=int((minX1-minX)/w_e_pixel1)
yOffset1=int((maxY1-maxY)/n_s_pixel1)
#计算test2在拼接图像左上角的偏移值
xOffset2=int((minX2-minX)/w_e_pixel1)
yOffset2=int((maxY2-maxY)/n_s_pixel1)

for i in range(5):#将图像1中每个波段写入拼接影像相对应的波段
    ds1_band=ds1.GetRasterBand(i+1)
    ds1_data=ds1_band.ReadAsArray(0,0,cols1,rows1)
    # 从拼接影像栅格数据集的xOffset1,yOffset1位置输入
    dataset_out.GetRasterBand(i+1).WriteArray(ds1_data,xOffset1,yOffset1)
    for j in range(5):#将图像2中每个波段写入拼接影像相对应的波段
        ds2_band = ds2.GetRasterBand(j+1)
        ds2_data = ds2_band.ReadAsArray(0, 0, cols2, rows2)
        dataset_out.GetRasterBand(j + 1).WriteArray(ds2_data, xOffset2, yOffset2)
#清空缓存，关闭数据集
ds1.FlushCache()
ds2.FlushCache()
dataset_out.FlushCache()
del ds1
del ds2
del dataset_out

拼接前

           

拼接后