遙感圖像的拼接和鑲嵌


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)

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM