栅格資料的批量鑲嵌（附Python腳本）

部落格小序：在資料處理的過程中，會遇到需要大量鑲嵌的情況，當資料較多時手動鑲嵌較為麻煩，自己最近對分省的DEM資料進行鑲嵌，由于利用python進行鑲嵌較為友善，特撰此博文以記之。

參考部落格：

https://blog.csdn.net/qq_15642411/article/details/79187787

https://blog.csdn.net/XBR_2014/article/details/85255412

1.腳本處理情況說明

本執行個體中，處理的資料是分省的DEM資料，每個省由若幹數量不同的tif檔案，本腳本主要的功能隻有一個即： 對分省的DEM資料進行分省鑲嵌

2.腳本代碼

#添加一個計時器
import time
start = time.time()

import os, shutil, glob
from osgeo import gdal

# 定義一個鑲嵌的函數，傳入的參數是需要鑲嵌的資料的清單，以及輸出路徑
def mosaic(data_list, out_path):

    # 讀取其中一個栅格資料來确定鑲嵌圖像的一些屬性
    o_ds = gdal.Open(data_list[0])
    # 投影
    Projection = o_ds.GetProjection()
    # 波段資料類型
    o_ds_array = o_ds.ReadAsArray()

    if 'int8' in o_ds_array.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in o_ds_array.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    # 像元大小
    transform = o_ds.GetGeoTransform()
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    del o_ds

    minx_list = []
    maxX_list = []
    minY_list = []
    maxY_list = []

    # 對于每一個需要鑲嵌的資料，讀取它的角點坐标
    for data in data_list:

        # 讀取資料
        ds = gdal.Open(data)
        rows = ds.RasterYSize
        cols = ds.RasterXSize

        # 擷取角點坐标
        transform = ds.GetGeoTransform()
        minX = transform[0]
        maxY = transform[3]
        pixelWidth = transform[1]
        pixelHeight = transform[5]  # 注意pixelHeight是負值
        maxX = minX + (cols * pixelWidth)
        minY = maxY + (rows * pixelHeight)

        minx_list.append(minX)
        maxX_list.append(maxX)
        minY_list.append(minY)
        maxY_list.append(maxY)

        del ds

    # 擷取輸出圖像坐标
    minX = min(minx_list)
    maxX = max(maxX_list)
    minY = min(minY_list)
    maxY = max(maxY_list)

    # 擷取輸出圖像的行與列
    cols = int((maxX - minX) / pixelWidth)
    rows = int((maxY - minY) / abs(pixelHeight))# 注意pixelHeight是負值

    # 計算每個圖像的偏移值
    xOffset_list = []
    yOffset_list = []
    i = 0

    for data in data_list:
        xOffset = int((minx_list[i] - minX) / pixelWidth)
        yOffset = int((maxY_list[i] - maxY) / pixelHeight)
        xOffset_list.append(xOffset)
        yOffset_list.append(yOffset)
        i += 1

    # 建立一個輸出圖像
    driver = gdal.GetDriverByName("GTiff")
    dsOut = driver.Create(out_path + ".tif", cols, rows, 1, datatype) 
    bandOut = dsOut.GetRasterBand(1)

    i = 0
    #将原始圖像寫入新建立的圖像
    for data in data_list:
        # 讀取資料
        ds = gdal.Open(data)
        data_band = ds.GetRasterBand(1)
        data_rows = ds.RasterYSize
        data_cols = ds.RasterXSize

        data = data_band.ReadAsArray(0, 0, data_cols, data_rows)
        bandOut.WriteArray(data, xOffset_list[i], yOffset_list[i])

        del ds
        i += 1

    # 設定輸出圖像的幾何資訊和投影資訊
    geotransform = [minX, pixelWidth, 0, maxY, 0, pixelHeight]
    dsOut.SetGeoTransform(geotransform)
    dsOut.SetProjection(Projection)

    del dsOut

def main():

    input_folder = "D:\\cnblogs\\data\\china"
    file_list = glob.glob(input_folder + "\\*")

    out_file = "D:\\cnblogs\\data\\china_moasic"
    if os.path.exists(out_file):
        shutil.rmtree(out_file)
        os.mkdir(out_file)
    else:
        os.mkdir(out_file)

    for file in file_list:
        basename = os.path.basename(file)
        out_path = out_file + "\\" + basename

        data_list = glob.glob(file + "\\" + "*.tif")
        print(data_list)

        try:
            mosaic(data_list, out_path)
            print(file + "鑲嵌結束")
        except:
            bad_list.append(file)
            print(file + "資料超過4G或其他原因導緻無法鑲嵌")

bad_list = []
main()
print("無法鑲嵌的檔案包括如下")
print (bad_list)

end = time.time()
print ("程式運作時間{:.2f}分鐘".format((end-start)/60.0))
           

3.問題總結

1.多幅栅格資料鑲嵌時，由于鑲嵌後的栅格資料大小超過4G，python的會報錯。

2.此腳本隻考慮了栅格資料隻有一個波段的情形，多波段栅格資料的鑲嵌未來有機會遇到再補充。

3.由于新生成的鑲嵌資料是規則的矩形，但是用于鑲嵌的資料不一定能剛好完美的覆寫新的圖層，導緻沒有資料用來鑲嵌的部分的值0，此問題較為重要，需要注意。

本文作者：DQTDQT

限于作者水準有限，如文中存在任何錯誤，歡迎不吝指正、交流。

聯系方式：

QQ:1426097423

e-mail：[email protected]

本文版權歸作者和部落格園共有，歡迎轉載、交流，但未經作者同意必須保留此段聲明，且在文章頁面明顯位置給出原文連結，如果覺得本文對您有益，歡迎點贊、探讨。

轉載于:https://www.cnblogs.com/xsman/p/11439894.html