博客小序：在数据处理的过程中，会遇到需要大量镶嵌的情况，当数据较多时手动镶嵌较为麻烦，自己最近对分省的DEM数据进行镶嵌，由于利用python进行镶嵌较为方便，特撰此博文以记之。
参考博客：1.脚本处理情况说明
本实例中，处理的数据是分省的DEM数据，每个省由若干数量不同的tif文件，本脚本主要的功能只有一个即： 对分省的DEM数据进行分省镶嵌
2.脚本代码
#添加一个计时器import timestart = time.time()import os, shutil, globfrom 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 dsOutdef 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本文版权归作者和博客园共有，欢迎转载、交流，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接，如果觉得本文对您有益，欢迎点赞、探讨。