GDAL+Python实现栅格影像处理之拼接镶嵌Mosaic

镶嵌概念

镶嵌是指将有重叠区域的多个图像根据其地理坐标将其拼接生成一个图像的过程。
实现原理：

• 1、计算各个待镶嵌图像的四至范围，然后对计算的各个待镶嵌的图像的四至范围求并得
  到整个镶嵌结果的四至范围
• 2、通过指定输出图像的分辨率以及计算的四至范围计算输出图像的大小，并创建输出图像
• 3、循环处理每个待镶嵌图像，将待镶嵌图像的像元读取出来并写入镶嵌结果图像中的对应位置
• 4、处理完所有的待镶嵌图像，关闭结果图像，清理资源等

使用方法

在讲述方法之前，列举一下numpy中的几个例子，因为下面我们会使用。

    a=numpy.array([(2,23,4),(23,3,4),(2345,56,8)])
    # 矩阵切片替换 将a矩阵2行 数据替换
    a[1,:]=(56,3,5) 
    b=numpy.array([(54,34,23)])
    # 矩阵增加维度（行数和列数） 在a矩阵下增加其行数
    print(numpy.r_[a,b]

• 方法1 ：通过镶嵌原理自定义方法实现。由于我的两幅图像位置固定，未根据其他情况进行拼接。读者可以参考这种方式，因为理解思路最重要。

代码实现

def RasterMosaic():
    ds=gdal.Open(inputfilePath,gdal.GA_ReadOnly)
    db=gdal.Open(referencefilefilePath,gdal.GA_ReadOnly)
    # todo 判断两影像相对位置（计算一副图像顶点对应在另一幅影像中的行列号）
    Originx = int((ds.GetGeoTransform()[0]-db.GetGeoTransform()[0])/db.GetGeoTransform()[1])
    Originy = int((ds.GetGeoTransform()[3]-db.GetGeoTransform()[3])/db.GetGeoTransform()[5])
    # 第一步判断顶点是否在其中
    if(Originx<db.RasterXSize and Originy<db.RasterYSize):
        print("顶点在其内")
        # 第二步判断右下角是否在其内
        OriginRightx=Originx+ds.RasterXSize
        OriginRighty=Originy+ds.RasterYSize
        # 根据两顶点判断其上下左右的情况
        if(OriginRightx<db.RasterXSize and OriginRighty<db.RasterYSize):
            print("右下角也在其内")
        elif(OriginRightx<db.RasterXSize and OriginRighty <db.RasterySize):
            print("右下角不在其内且右边界未超出最大边界")
            # todo 将两幅影像公共部分分别都表示出来
            data1 = ds.GetRasterBand(1).ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)[0:db.RasterYSize - Originy,
                    0:ds.RasterXSize]
            data2 = db.GetRasterBand(1).ReadAsArray(0, 0, db.RasterXSize, db.RasterYSize)[
                    Originy:db.RasterYSize - Originy, Originx:Originx + ds.RasterXSize]
            # 表示出在外部的那一部分数据
            dataexcept = ds.GetRasterBand(1).ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)[
                         db.RasterYSize - Originy:ds.RasterYSize, 0:ds.RasterXSize]
            print(ds.RasterYSize - (db.RasterYSize - Originy))
            datacommon = numpy.zeros(shape=(ds.RasterYSize - (db.RasterYSize - Originy), db.RasterYSize), dtype=int)
            datacommon[:, Originx:ds.RasterXSize + Originx] = dataexcept
            data = numpy.r_[db.GetRasterBand(1).ReadAsArray(0, 0, db.RasterXSize, db.RasterYSize), datacommon]
        elif(OriginRightx>db.RasterXSize and OriginRighty<db.RasterYSize):
            print("右下角不在其内且右边界超出最大边界并且下边界未超出最大边界")
        elif(OriginRightx>db.RasterXSize and OriginRighty>db.RasterYSize):
            print("右下角不在其内且右边界超出最大边界并且下边界未超出最大边界")
            width=db.RasterYSize
            height=OriginRighty
        outputfilePath = 'G:/studyprojects/gdal/GdalStudy/Files/images/RasterMosaic.tif'
        ds = gdal.Open(inputfilePath, gdal.GA_ReadOnly)
        driver = gdal.GetDriverByName('GTiff')
        output = driver.Create(outputfilePath, Originx+ds.RasterXSize, db.RasterYSize, db.RasterCount, db.GetRasterBand(1).DataType)
        output.SetGeoTransform(db.GetGeoTransform())
        output.SetProjection(db.GetProjection())
        for i in range(db.RasterCount):
            # todo 将两幅影像公共部分分别都表示出来
            data1 = ds.GetRasterBand(i+1).ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)[0:ds.RasterYSize,
                    0:db.RasterXSize - Originx]
            data2 = db.GetRasterBand(i+1).ReadAsArray(0, 0, db.RasterXSize, db.RasterYSize)[
                    Originy:ds.RasterYSize + Originy, Originx:db.RasterXSize]
            # 表示出在外部的那一部分数据  这里是我已知我的图像位置，只写了一种情况。
            dataexcept = ds.GetRasterBand(i+1).ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)[0:ds.RasterYSize,
                         db.RasterXSize - Originx:ds.RasterXSize]
            datacommon = numpy.zeros(shape=(db.RasterYSize, ds.RasterXSize - (db.RasterXSize - Originx)), dtype=int)
            datacommon[Originy:ds.RasterYSize + Originy, :] = dataexcept
            data = numpy.c_[db.GetRasterBand(i+1).ReadAsArray(0, 0, db.RasterXSize, db.RasterYSize), datacommon]
            output.GetRasterBand(i + 1).WriteArray(data, 0, 0)
            output.GetRasterBand(i + 1).SetNoDataValue(0)  # todo 给各波段设置nodata值
            output.GetRasterBand(i + 1).FlushCache()  # todo 波段统计量
    else:
        print("顶点不在其内")

• 方法2：方法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)

效果展示

原始两幅图像位置：

新的影像：

如动图所示：