1 自定义镶嵌函数

遥感图像的镶嵌，主要分为5大步骤：
step1： 1）对于每一幅图像，计算其行与列；
2）获取左上角X,Y
3）获取像素宽和像素高
4）计算max X 和 min Y，切记像素高是负值

maxX1 = minX1 + (cols1 * pixelWidth)
minY1 = maxY1 + (rows1 * pixelHeight）

step2 ：计算输出图像的min X ,max X,min Y,max Y

minX = min(minX1, minX2, …)
maxX = max(maxX1, maxX2, …)

y坐标同理
step3：计算输出图像的行与列

 cols = int((maxX – minX) / pixelWidth)
 rows = int((maxY – minY) / abs(pixelHeight)

step 4：创建一个输出图像
driver.create()
step 5：1）计算每幅图像左上角坐标在新图像的偏移值
2）依次读入每幅图像的数据并利用1)计算的偏移值将其写入新图像中

#mosica 两张图像
import os, sys, gdal
from gdalconst import *
os.chdir(r'F:\algorithm\算法练习\拼接与镶嵌')#改变文件夹路径
# 注册gdal(required)
gdal.AllRegister()
# 读入第一幅图像
ds1 = gdal.Open(r'F:\algorithm\算法练习\拼接与镶嵌\test2_next\4.tif')
band1 = ds1.GetRasterBand(1)
rows1 = ds1.RasterYSize
cols1 = ds1.RasterXSize
bands = ds1.RasterCount
# 获取图像角点坐标
transform1 = ds1.GetGeoTransform()
minX1 = transform1[0]
maxY1 = transform1[3]
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5]#是负值（important）
maxX1 = minX1 + (cols1 * pixelWidth1)
minY1 = maxY1 + (rows1 * pixelHeight1)
# 读入第二幅图像
ds2 = gdal.Open(r'F:\algorithm\算法练习\拼接与镶嵌\test2_next\5.tif')
rows2 = ds2.RasterYSize
cols2 = ds2.RasterXSize
# 获取图像角点坐标
transform2 = ds2.GetGeoTransform()
minX2 = transform2[0]
maxY2 = transform2[3]
pixelWidth2 = transform2[1]
pixelHeight2 = transform2[5]
maxX2 = minX2 + (cols2 * pixelWidth2)
minY2 = maxY2 + (rows2 * pixelHeight2)
# 获取输出图像坐标
minX = min(minX1, minX2)
maxX = max(maxX1, maxX2)
minY = min(minY1, minY2)
maxY = max(maxY1, maxY2)
#获取输出图像的行与列
cols = int((maxX - minX) / pixelWidth1)
rows = int((maxY - minY) / abs(pixelHeight1))
# 计算图1左上角的偏移值（在输出图像中）
xOffset1 = int((minX1 - minX) / pixelWidth1)
yOffset1 = int((maxY1 - maxY) / pixelHeight1)
# 计算图2左上角的偏移值（在输出图像中）
xOffset2 = int((minX2 - minX) / pixelWidth1)
yOffset2 = int((maxY2 - maxY) / pixelHeight1)
# 创建一个输出图像
driver = ds1.GetDriver()
dsOut = driver.Create('mosiac.tif', cols, rows, bands, band1.DataType)#1是bands，默认
# 读图1的数据并将其写到输出图像中
data1 = ds1.ReadAsArray(0, 0, cols1, rows1)
for i in range( bands):
    dsOut.GetRasterBand(i+1).WriteArray(data1[i],xOffset1, yOffset1)
#读图2的数据并将其写到输出图像中
data2 = ds2.ReadAsArray(0, 0, cols2, rows2)
for i in range( bands):
    dsOut.GetRasterBand(i+1).WriteArray(data2[i],xOffset2, yOffset2)
''' 写图像步骤'''
#第二个参数是1的话：整幅图像重度，不需要统计
# 设置输出图像的几何信息和投影信息
geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
dsOut.SetGeoTransform(geotransform)
dsOut.SetProjection(ds1.GetProjection())

2 采用gdal.Warp()提供的接口进行镶嵌

def RasterMosaic():
    print("图像拼接")
    inputrasfile1 = gdal.Open(inputfilePath, gdal.GA_ReadOnly) # 第一幅影像
    inputProj1 = inputrasfile1.GetProjection()
    inputrasfile2 = gdal.Open(referencefilefilePath, gdal.GA_ReadOnly) # 第二幅影像
    inputProj2 = inputrasfile2.GetProjection()
    outputfilePath = 'G:/studyprojects/gdal/GdalStudy/Files/images/RasterMosaic2.tif'
    options=gdal.WarpOptions(srcSRS=inputProj1, dstSRS=inputProj1,format='GTiff',resampleAlg=gdalconst.GRA_Bilinear)
    gdal.Warp(outputfilePath,[inputrasfile1,inputrasfile2],options=options)