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

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

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

文章目录

• 遥感数据处理系列
• 前言
• 一、栅格数据重采样
• • 1. 原理简介
  • 2. 代码
• 二、多文件夹场景
• 总结
• • • 后记

多源遥感数据的使用如何快速计算？处理为相同维度的栅格数据，那么，一个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牛皮！毕业万岁！中期快乐！
定量遥感中的重采样、按掩膜提取等操作均会对像元值产生一定改变，应注意。尤其是在重采样时，要考虑不同插值方式对结果的影响。

后记

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