栅格数据的批量镶嵌（附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：duanquntaoyx@163.com

本文版权归作者和博客园共有，欢迎转载、交流，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接，如果觉得本文对您有益，欢迎点赞、探讨。

巴特西

栅格数据的批量镶嵌（附Python脚本）

1.脚本处理情况说明

2.脚本代码

3.问题总结

最新文章

热门文章