遥感图像的拼接和镶嵌

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)计算的偏移值将其写入新图像中

import os
import glob
import math

import gdal

def get_extent(fn):
    ds = gdal.Open(fn)
    gt = ds.GetGeoTransform()

    r = ds.RasterYSize
    c = ds.RasterXSize
    return (gt[0], gt[3], gt[0] + gt[1] * c, gt[3] + gt[5] * r)


def mosiac(in_files):
    min_X, max_Y, max_X, min_Y = get_extent(in_files[0])
    for fn in in_files[1:]:
        minx, maxy, maxx, miny = get_extent(fn)
        min_X = min(min_X, minx)
        max_Y = max(max_Y, maxy)
        max_X = max(max_X, maxx)
        min_Y = min(min_Y, miny)


    ds1 = gdal.Open(in_files[0])
    band1 = ds1.GetRasterBand(1)
    transform1 = ds1.GetGeoTransform()
    pixelWidth1 = transform1[1]
    pixelHeight1 = transform1[5]
    bands = ds1.RasterCount

    # 获取输出图像的行与列
    cols = int((max_X - min_X) / pixelWidth1)
    rows = int((max_Y - min_Y) / abs(pixelHeight1))

    driver = ds1.GetDriver()
    dsOut = driver.Create(r'F:\algorithm\算法练习\拼接与镶嵌\mosiac1.tif', cols, rows, bands, band1.DataType)#1是bands，默认

    for file in in_files:
        ds2 = gdal.Open(file)
        rows2 = ds2.RasterYSize
        cols2 = ds2.RasterXSize
        data2 = ds2.ReadAsArray(0, 0, cols2, rows2)
        transform2 = ds2.GetGeoTransform()
        minX2 = transform2[0]
        maxY2 = transform2[3]
        # 计算每张图左上角的偏移值（在输出图像中）
        xOffset2 = int((minX2 - min_X) / pixelWidth1)
        yOffset2 = int((maxY2 - max_Y) / pixelHeight1)
        for i in range( bands):
            dsOut.GetRasterBand(i+1).WriteArray(data2[i],xOffset2, yOffset2)

    geotransform = [min_X, pixelWidth1, 0, max_Y, 0, pixelHeight1]
    dsOut.SetGeoTransform(geotransform)
    dsOut.SetProjection(ds1.GetProjection())


if __name__ == '__main__':
    os.chdir(r'F:\algorithm\算法练习\拼接与镶嵌\test2_next')
    in_files = glob.glob('*.tif')
    print(in_files)
    mosiac(in_files )

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)