快捷導(dǎo)航

Python基于GDAL鑲嵌拼接遙感影像

更新時(shí)間：2023年10月23日 14:42:21 作者：RS迷途小書(shū)童

這篇文章主要介紹了Python基于GDAL鑲嵌拼接遙感影像, 這里有一點(diǎn)需要注意的就是,用這個(gè)方法進(jìn)行鑲嵌拼接操作時(shí),影像有一條明顯的拼接線(xiàn),不知道是不是我數(shù)據(jù)的問(wèn)題,你們可以自己嘗試一下,只要修改主函數(shù)中的路徑即可,需要的朋友可以參考下

沒(méi)啥好說(shuō)的，處理高分辨率影像時(shí)，數(shù)據(jù)高達(dá)幾十G。用ENVI或者ArcGIS進(jìn)行影像的拼接時(shí)，往往會(huì)出現(xiàn)未響應(yīng)的情況。出現(xiàn)未響應(yīng)的話(huà)，運(yùn)氣好等個(gè)一晚上可能會(huì)動(dòng)一動(dòng)，運(yùn)氣不好就等著強(qiáng)制關(guān)閉重做吧。

所以搞了一個(gè)Python進(jìn)行拼接操作的代碼，雖然速度不算快，但至少不會(huì)未響應(yīng)。同時(shí)如果對(duì)代碼進(jìn)行一些改進(jìn)，還可以進(jìn)行批量拼接的操作，百利而無(wú)一害。

一、導(dǎo)入GDAL庫(kù)

from osgeo import gdal

二、查看影像信息

為了湊字?jǐn)?shù)的，可以查看影像的投影、寬度、高度、波段數(shù)等信息。不過(guò)需要注意的是在ENVI中沒(méi)有投影坐標(biāo)系，只有地理坐標(biāo)系是做不了鑲嵌拼接的。

這個(gè)代碼我還不太清楚能不能不要投影坐標(biāo)系進(jìn)行拼接，你們可以自己試試。但最好還是用包含投影坐標(biāo)系的影像進(jìn)行拼接。所以在拼接之前就可以用這段代碼先看一看。

def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打開(kāi)數(shù)據(jù)集dataset
    ds_width = ds.RasterXSize  # 獲取數(shù)據(jù)寬度
    ds_height = ds.RasterYSize  # 獲取數(shù)據(jù)高度
    ds_bands = ds.RasterCount  # 獲取波段數(shù)
    ds_geo = ds.GetGeoTransform()  # 獲取仿射地理變換參數(shù)
    ds_prj = ds.GetProjection()  # 獲取投影信息
    print("影像的寬度為：" + str(ds_width))
    print("影像的高度為：" + str(ds_height))
    print("仿射地理變換參數(shù)為：" + str(ds_geo))
    print("投影坐標(biāo)系為：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以數(shù)組的形式讀取整個(gè)數(shù)據(jù)集

三、鑲嵌模塊

這里用到了GDAL中的Warp函數(shù)，是不是有點(diǎn)熟悉。之前裁剪也是用的這個(gè)函數(shù)，不得不說(shuō)這個(gè)函數(shù)是真的強(qiáng)大，之前有文章介紹了其中的函數(shù)，大家有興趣可以去看下，同時(shí)給個(gè)贊吧！【Python&RS】GDAL批量裁剪遙感影像/柵格數(shù)據(jù)

代碼中的srcSRS,dstSRS分別是輸入投影和輸出投影，這里用一樣的就行了。因?yàn)槲覀冏龅氖氰偳恫僮?，肯定是不用?dòng)原始坐標(biāo)系的。其他的參數(shù)都在代碼中表明了，這里就不介紹了，如果大家有什么問(wèn)題，可以留言交流。

def Mosaic_GDAL(path_image1, path_image2, path_out):
    """
    :param path_image1: 需要鑲嵌的影像
    :param path_image2: 需要鑲嵌的影像
    :param path_out: 鑲嵌后輸出的影像路徑
    :return: None
    """
    image1 = gdal.Open(path_image1, gdal.GA_ReadOnly)  # 第一幅影像
    input_proj = image1.GetProjection()
    image2 = gdal.Open(path_image2, gdal.GA_ReadOnly)  # 第二幅影像
    options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                               resampleAlg=gdal.GRA_NearestNeighbour,callback=Show_Progress)
    # 輸入投影，輸出投影，輸出格式，重采樣方法
    gdal.Warp(path_out, [image1, image2], options=options)
    # 輸出路徑，需要鑲嵌的數(shù)據(jù)，參數(shù)配置
    ds = gdal.Open(path_out, gdal.GA_ReadOnly)
    ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
    # 創(chuàng)建金字塔
    del image1, image2, ds

四、回調(diào)函數(shù)（沒(méi)啥用）

上一篇博文已經(jīng)介紹過(guò)了，就兩點(diǎn)：1.水字?jǐn)?shù)，2.記錄進(jìn)度。省的看著代碼發(fā)呆。

def Show_Progress(percent, msg, tag):
    """
    :param percent: 進(jìn)度，0~1
    :param msg:
    :param tag:
    :return:
    """
    if 25 <= percent*100 <= 26:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")
    if 50 <= percent*100 <= 51:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")
    if 75 <= percent*100 <= 76:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")

五、完整代碼

# -*- coding: utf-8 -*-
"""
@Time ： 2023/6/25 16:28
@Auth ： RS迷途小書(shū)童
@File ：Raster Data Mosaic.py
@IDE ：PyCharm
@Purpose：遙感影像、柵格數(shù)據(jù)的鑲嵌拼接（有明顯拼接線(xiàn)）
"""
from osgeo import gdal
def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打開(kāi)數(shù)據(jù)集dataset
    ds_width = ds.RasterXSize  # 獲取數(shù)據(jù)寬度
    ds_height = ds.RasterYSize  # 獲取數(shù)據(jù)高度
    ds_bands = ds.RasterCount  # 獲取波段數(shù)
    ds_geo = ds.GetGeoTransform()  # 獲取仿射地理變換參數(shù)
    ds_prj = ds.GetProjection()  # 獲取投影信息
    print("影像的寬度為：" + str(ds_width))
    print("影像的高度為：" + str(ds_height))
    print("仿射地理變換參數(shù)為：" + str(ds_geo))
    print("投影坐標(biāo)系為：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以數(shù)組的形式讀取整個(gè)數(shù)據(jù)集
def Show_Progress(percent, msg, tag):
    """
    :param percent: 進(jìn)度，0~1
    :param msg:
    :param tag:
    :return:
    """
    if 25 <= percent*100 <= 26:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")
    if 50 <= percent*100 <= 51:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")
    if 75 <= percent*100 <= 76:
        print("進(jìn)度：" + "%.2f" % (percent*100) + "%")
def Mosaic_GDAL(path_image1, path_image2, path_out):
    """
    :param path_image1: 需要鑲嵌的影像
    :param path_image2: 需要鑲嵌的影像
    :param path_out: 鑲嵌后輸出的影像路徑
    :return: None
    """
    image1 = gdal.Open(path_image1, gdal.GA_ReadOnly)  # 第一幅影像
    input_proj = image1.GetProjection()
    image2 = gdal.Open(path_image2, gdal.GA_ReadOnly)  # 第二幅影像
    options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                               resampleAlg=gdal.GRA_NearestNeighbour, callback=Show_Progress)
    # 輸入投影，輸出投影，輸出格式，重采樣方法
    gdal.Warp(path_out, [image1, image2], options=options)
    # 輸出路徑，需要鑲嵌的數(shù)據(jù)，參數(shù)配置
    ds = gdal.Open(path_out, gdal.GA_ReadOnly)
    ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
    # 創(chuàng)建金字塔
    del image1, image2, ds
if __name__ == "__main__":
    path1 = "B:/proj_clip"
    path2 = "B:/2_proj_clip"
    path_output = "warp1.tif"
    print("開(kāi)始拼接......")
    Mosaic_GDAL(path1, path2, path_output)
    print("拼接完成......")