ArcGIS批量拼接栅格数据（ArcPy方法）

遥感数据处理系列

一些项目及科研中遇到的小需求，一方面记录自己的学习历程，另一方面帮助大家学习。本系列文章的开发环境为：ArcGIS 10.2.2 + Python 2.7、ENVI 5.3 + IDL 8.5

ArcGIS批量计算栅格数据平均值（ArcPy方法）
GLDAS数据下载及处理（NC转TIF）
ArcGIS批量裁剪栅格数据
ArcGIS批量栅格重采样（ArcPy方法）
ArcGIS批量裁剪栅格数据（ArcPy方法）

IDL多进程批处理遥感数据
ArcGIS批量拼接栅格数据（ArcPy方法）

---

---

前言

如何拼接两幅相邻的遥感影像？一大堆的遥感数据如何批量拼接？又如何把一个文件夹里的所有栅格数据拼接为一景影像？两个文件夹场景呢？那么，一个NB的批处理脚本派上了用场！祭出ArcPy

---

一、栅格数据拼接

1. 原理简介

栅格数据拼接主要使用ArcPy的MosaicToNewRaster函数。

函数使用：

MosaicToNewRaster(input_rasters, output_location, raster_filename, {coordinate_system}, {pixel_type}, {cellsize}, number_of_bands, {mosaic_method}, {mosaic_colormap_mode})

常用参数简介：

	input_rasters：待拼接的栅格数据集。须为相同的波段数和相同的位深度
	output_location：文件输出路径的文件夹
	raster_filename：输出文件名
	coordinate_system：输出栅格数据集的坐标系
	pixel_type：镶嵌数据集的位深度
		-> 8_BIT_UNSIGNED — 8位无符号数据类型。支持的值为 0 到 255
		-> 8_BIT_SIGNED   — 8位有符号数据类型。支持的值为 -128 到 127
		-> 16_BIT_UNSIGNED — 16位无符号数据类型。取值范围为 0 到 65,535
		-> 16_BIT_SIGNED  — 16位有符号数据类型。取值范围为 -32,768 到 32,767
		-> 32_BIT_UNSIGNED — 32位无符号数据类型。取值范围为 0 到 4,294,967,295
		-> 32_BIT_SIGNED — 32位有符号数据类型。取值范围为 -2,147,483,648 到 2,147,483,647
		-> 32_BIT_FLOAT — 支持小数的32位数据类型
		-> 64_BIT — 支持小数的64位数据类型
	cellsize：像元大小
	number_of_bands：输出栅格的波段数
	mosaic_method：镶嵌重叠的方法，设置重叠区域像元值
		->FIRST  —  第一个栅格数据集中的值
		->LAST   —  最后一个栅格数据集中的值（默认）
		->BLEND  —  各像元值的水平加权计算结果
		->MEAN   —  叠置像元的平均值
		->MINIMUM — 叠置像元的最小值
		->MAXIMUM — 叠置像元的最大值
		->SUM    —  叠置像元的总和
	mosaic_colormap_mode：输出的色彩映射表进行选择

注意：栅格拼接的参数较多，选择时须核实下数据类型。

2. 代码

文件组织架构：

输入：

1. 一个含有若干栅格数据的文件夹 inws_1 （本例为“.tif”格式）
2. 另一个文件夹 inws_2 中的文件数量与顺序和 inws_1 中相同

输出： 在输出路径文件夹下重构文件夹组织形式，并生成若干合并空间范围的tif格式数据

代码实例：

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

'''
适用于：两个文件夹场景
思路：
    两个输入路径 -> 两个范围的栅格数据，两个文件夹中的文件数量及顺序相同，仅空间范围不同
'''

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

# 输入路径  应该注意，中文路径，会导致读不出文件
inws_1 = r"C:\Mosaic\BottomRight"
inws_2 = r"C:\Mosaic\TopLeft"

# 输出路径
outws = r"C:\Mosaic\Out"

path_1_list = os.listdir(inws_1)  # 存了三年的数据，里面有三个文件夹，每个文件夹直接放了每年的数据
path_2_list = os.listdir(inws_2)  # 存了三年的数据，里面有三个文件夹，每个文件夹直接放了每年的数据

# 进入二级目录
for i in range(len(path_1_list)):
    path_1 = inws_1 + "\\" + path_1_list[i]  # 年份文件夹的路径
    path_2 = inws_2 + "\\" + path_2_list[i]

    # 重构输出文件夹
    outFolder = outws + "\\" + path_1_list[i]
    isExists = os.path.exists(outFolder)
    if not isExists:
        os.makedirs(outFolder)

    # 利用glob包，将三级目录下的所有tif文件读存放到rasters中
    rasters_1 = glob.glob(os.path.join(path_1, "*.tif"))
    rasters_2 = glob.glob(os.path.join(path_2, "*.tif"))

    base = rasters_1[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

    rasters = []  # 提取待拼接影像的文件名，且中间以;隔开，例如：a.tif;b.tif;c.tif

    # 循环rasters中的所有影像，进行累加操作
    for ii in range(len(rasters_1)):    # 这里是，年份文件夹里的栅格数据数量，两个文件夹的第一个文件进行拼接，第二个进行拼接，...
        rasters.append(rasters_1[ii])
        rasters.append(rasters_2[ii])

        ras_list = ";".join(rasters)  # 字符串拼接

        nameT = os.path.basename(base)+ '.tif'

        # 执行拼接操作
        arcpy.MosaicToNewRaster_management(ras_list, outFolder, nameT, out_coor_system,
                                           "32_BIT_FLOAT",
                                           cellwidth, bandcount, "LAST", "FIRST")
        del rasters[:]   # 清空List   注意： del a[:]会删除列表中的所有元素，但是不会删除列表a;       del a 不仅会删除列表中的所有元素，同时会删除列表a
        print os.path.basename( nameT ) + "         ---- OK！ ----         "

    print path_1_list[i] + "         ---- OK！ ----         \n"

print("   ---   All project is OK!   ---   ")

上例可实现对两个输入路径文件夹下的所有栅格数据的按顺序拼接

二、单文件夹场景

在处理MYD09时，遇到了没有坐标的情况，而直接在NASA官网选择预处理后的输出结果存在拼接错误。故而仿MRT处理逻辑搞了个单文件夹拼接，首先把同一天的所有数据放到一个文件夹下，然后，获取该文件夹内的所有栅格数据进行拼接操作

输入路径文件组织架构：

inws  -> 2008、2009、2010...（若干年份文件夹里存放了若干tif文件）

代码如下：

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

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

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

# 输入路径  应该注意，中文路径，会导致读不出文件
inws = r"C:\Mosaic\BottomRight"

# 输出路径
outws = r"C:\Mosaic\Out"

# 以下逻辑为：
#   获得子目录全路径，并进入子目录
path_list = os.listdir(inws)  # 存了三年的数据，里面有三个文件夹，每个文件夹直接放了每年的数据

# 进入二级目录
for i in range(len(path_list)):
    path = inws + "\\" + path_list[i]  # 年份文件夹的路径

    # 重构输出文件夹
    outFolder = outws + "\\" + path_list[i]
    isExists = os.path.exists(outFolder)
    if not isExists:
        os.makedirs(outFolder)

    # 利用glob包，将三级目录下的所有tif文件读存放到rasters中
    rasters = glob.glob(os.path.join(path, "*.tif"))

    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 = os.path.basename(base) + '.tif'

    # 执行拼接操作
    arcpy.MosaicToNewRaster_management(mosaic_ras, outFolder, nameT, out_coor_system,
                                       "32_BIT_FLOAT",
                                       cellwidth, bandcount, "LAST", "FIRST")
    print os.path.basename(nameT) + "         ---- OK！ ----         "

print("   ---   All project is OK!   ---   ")

注意：这里只是展示在合并单文件夹内所有文件的业务逻辑，并不是处理MYD09没坐标的那个数据的代码。不过，基本思路是相同的

总结

ArcPy牛皮！毕业万岁！中期快乐！
批量合并栅格数据，ArcPy太棒了吧，python真香喽~

后记

写博客的初衷是分享我的一些经验，同时也方便自己在其他电脑上进行数据处理。帮了很多人，但评论区小伙伴也有遇到问题的，那么：知识付费，我的时间和经验正好可以解决你的问题。