ArcGIS批量栅格重采样（ArcPy方法）

遥感数据处理系列

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

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

---

---

前言

多源遥感数据的使用如何快速计算？处理为相同维度的栅格数据，那么，一个NB的批处理脚本派上了用场！祭出ArcPy

---

一、栅格数据重采样

1. 原理简介

栅格数据重采样主要使用ArcPy的Resample_management函数。

函数使用：

Resample_management (in_raster, out_raster, {cell_size}, {resampling_type})

常用参数简介：

	in_raster：想要更改空间分辨率的栅格数据集
	out_raster：要创建的数据集的名称、位置和格式
	cell_size：使用现有栅格数据集的新栅格的像元大小或指定其宽度 (x) 和高度 (y)
	resampling_type：根据待处理数据类型选择相应的重采样方式
	        NEAREST — 适用于离散数据，例如土地覆被   |    处理速度最快，对像素值的更改最小（常用）
	        BILINEAR — 适用于连续数据               |    平均化（距离权重）周围 4 个像素的值
	        CUBIC — 适用于连续数据                  |    拟合周围16个像素的平滑曲线
	        MAJORITY — 适用离散数据                 |    赋值像元 3 x 3 窗口中的最高频值

注意：重采样时应根据不同的需求进行插值方式选取，尽量避免对像素值的改变，从而避免对地物光谱信息或分类属性的错解。

2. 代码

输入： 一个含有若干栅格数据的文件夹（本例为“.tif”格式）

输出： 在输出路径文件夹下生成若干tif格式数据

参考像元： 本例为MYD09A1的500m分辨率大小

代码实例：

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

import arcpy

# 输入路径  应该注意，中文路径，会导致读不出文件
inws = r"F:\Mekong\TIF-8day\2017"

# 参考文件路径    使用栅格数据集（从其导入方形像元大小）的路径
inf = r"H:\MYD09A1\b1\MYD09A1_A2008001.MekongRiver_basin.sur_refl_b01.tif"

cellsize = "{0} {1}".format(arcpy.Describe(inf).meanCellWidth,arcpy.Describe(inf).meanCellHeight)

# 输出路径
outws = r"F:\TIF-8day-re\2017"

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

# 循环rasters中的所有影像，进行按掩模提取操作
for raster in rasters:
    print str(raster)

    nameT = os.path.basename(raster).split(".")[0] + "_re.tif"	# 自定义文件名
    outname = os.path.join(outws, nameT)  # 合并输出文件名+输出路径

    arcpy.Resample_management(raster, outname, cellsize, "NEAREST")

print("OK!")

上例可实现对输入路径文件夹下的所有栅格数据的重采样，但多年的数据都需要重采样怎么办？
修改逻辑为进入二级目录、三级目录…然后调用核心函数即可。

二、多文件夹场景

栅格重采样费时、费事、费心思？想挂机处理？Perfect，继续来！

代码如下：

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

# 输入路径  应该注意，中文路径，会导致读不出文件
inws = r"E:\TIF-8day-glt"

# 参考文件路径    使用栅格数据集（从其导入方形像元大小）的路径
inf = r"J:\MYD09A1\b1\MYD09A1_A2008001.MekongRiver_basin.sur_refl_b01.tif"
cellsize = "{0} {1}".format(arcpy.Describe(inf).meanCellWidth,arcpy.Describe(inf).meanCellHeight)

# 输出路径	需做到输入文件夹组织形式 与 输出文件夹组织形式相同（可进一步自动化，懒）
outws = r"E:\TIF-8day-glt-re"

# 以下逻辑为：
#   获得子目录全路径，并进入子目录
#   执行文件夹批量重采样操作
path_list = os.listdir(inws)  # 如果这里是存了三年的数据，那里面有三个文件夹，每个文件夹直接放了每年的数据
target_path_list = os.listdir(outws)  # 输出路径，也是里面有三个年份的文件夹

for i in range(len(path_list)):
    path = inws + "\\" + path_list[i]  # 文件夹列表   # 那这里就得到了年份的路径
    target_path = outws + "\\" + target_path_list[i]  # 这里，应该是进入了年份的路径，是处理后的数据存放的地方

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

    # 循环rasters中的所有影像，进行重采样操作
    for raster in rasters:
        print str(raster)

        nameT = os.path.basename(raster) + "_re.tif"       # 构造 输出文件名
        outname = os.path.join(target_path, nameT)                       # 合并输出文件名+输出路径 = 输出文件全路径
        arcpy.Resample_management(raster, outname, cellsize, "NEAREST")  # 执行操作

    print(path_list[i] + " --- OK!")

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

总结

ArcPy牛皮！毕业万岁！中期快乐！
定量遥感中的重采样、按掩膜提取等操作均会对像元值产生一定改变，应注意。尤其是在重采样时，要考虑不同插值方式对结果的影响。

后记

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