Skip to content

档案馆借的地形数据是分幅的,使用python批量镶嵌

使用 arcpy.MosaicToNewRaster_management方法。

python

# -*- coding: UTF-8 -*-
import arcpy
import glob
import os


'''
数据目录如下:
dem
    |-- 分幅1
        tif数据
    |-- 分幅2
        tif数据
'''

'''
思路:
    两个输入路径 -> 两个范围的栅格数据
'''

arcpy.CheckOutExtension("Spatial")  # 检查是否注册该模块

# 输入路径  应该注意,中文路径,会导致读不出文件
inws = r"D:\dem\G48"

# 输出路径
outws = r"Z:\lm\dem"

# 以下逻辑为:
#   获得子目录全路径,并进入子目录
path_list = os.listdir(inws)  

rasters = []

# 进入二级目录
for i in range(len(path_list)):
    path = inws + "\\" + path_list[i]
    ras= glob.glob(os.path.join(path, "*.bil"));
    rasters.append(ras[0])
    # 因为不太熟悉python,在for循环里没法修改全局变量
    # 4090是指该文件夹下的全部数据
    if(i == 4090):

        base = rasters[0]  # 指定该工作空间下的一副影像为基础影像,为后面的参数提取做准备
        out_coor_system = arcpy.Describe(base).spatialReference  # 坐标系统
        dataType = arcpy.Describe(base).DataType  # 数据类型
        piexl_type = arcpy.Describe(base).pixelType  # 像素类型
        cellwidth = arcpy.Describe(base).meanCellWidth  # 像元宽度
        bandcount = arcpy.Describe(base).bandCount  # 获取bandCount
        
        mosaic_ras = []  # 待拼接影像的文件名,且中间以;隔开,例如:a.tif;b.tif;c.tif
        for ras in rasters:
            mosaic_ras.append(ras)
        
        mosaic_ras = ";".join(rasters)  # 字符串拼接
        
        nameT = 'hebing' + '.bil'
        
        # 执行拼接操作
        arcpy.MosaicToNewRaster_management(mosaic_ras, outws, nameT, out_coor_system,
                                           "32_BIT_FLOAT",
                                           cellwidth, bandcount, "LAST", "FIRST")
        print os.path.basename(nameT) + "         ---- OK! ----         "
        
        print("   ---   All project is OK!   ---   ")