Arcpy-重采样记录

===================================================================================
记录缘由：在使用MRT将MCD43A3进行镶嵌裁剪后，发现分辨率有点不满足研究，遂采用arcpy进行重采样；
===================================================================================

执行代码的软件版本：ArcGIS 10.2 自带的 IDLE python (GUI)——python 2. 7 .3
该代码使用的数据前提是tif格式的，没有尝试没有投影的tif，应该也可以重采样成功，应该可能大概没有大的影响，尝试一下就知道了！

一、代码记录——单一文件中所有的tif进行批量重采样

# -*- coding: utf-8 -*-
import arcpy
import os
from arcpy import env
import time
import sys
import codecs# ==============================
# 配置部分
# ==============================
input_folder = r"Z:\mcd43a3_mrt_rep\2003\reproject"         # 输入栅格目录
output_folder = r"Z:\mcd43a3_mrt_rep_res0.1\2003"           # 输出栅格目录
resample_method = "NEAREST"                                 # 重采样方法: NEAREST, BILINEAR, CUBIC
output_cell_size = "0.1 0.1"                                # 重采样的像元大小（如：30米）# ==============================
# 设置控制台输出编码///中文在控制台打印的中文仍然是乱码，使用“gbk”格式也是，不重要，就懒得整了
# ==============================
try:sys.stdout = codecs.getwriter("utf-8")(sys.stdout.detach())
except Exception:pass  # 某些环境可能不支持该方式，可忽略# ==============================
# 函数定义
# ==============================
def resample_raster(in_raster, out_raster, cell_size, method):try:arcpy.Resample_management(in_raster, out_raster, cell_size, method)return True, Noneexcept Exception as e:return False, str(e)def batch_resample():arcpy.env.workspace = input_folderrasters = arcpy.ListRasters()total = len(rasters)print("\n start resampled，total: {} files\n".format(total))if not os.path.exists(output_folder):os.makedirs(output_folder)for idx, raster in enumerate(rasters):input_path = os.path.join(input_folder, raster)output_path = os.path.join(output_folder, raster)  # 使用原始文件名print("[{}/{}] processing : {} ...".format(idx + 1, total, raster))start_time = time.time()success, error = resample_raster(input_path, output_path, output_cell_size, resample_method)if success:print("    -> Finished，run time {:.2f} 秒\n".format(time.time() - start_time))else:print("    -> fail file: {}\n".format(error))if __name__ == '__main__':batch_resample()

====================================================================================

二、如何在ArcMap中运行代码

在下图找到ArcGIS，并在其中找到IDLE（python GUI） ----

接着转到如下：

====================================================================================
====================================================================================

三、数据量大，文件夹多-批量运行代码展示

# -*- coding: utf-8 -*-
import arcpy
import os
from arcpy import env
import time
import sys
import codecs# ==============================
# 配置部分——以及简要信息说明
# =============================
#
# tif文件夹的路径信息如下，代码可以直接遍历子目录，只要找到了tif就会进行重采样，
# 但是输出路径有一个问题，设置的输出路径只有一个，所有年份\*的文件夹最后保证是只有一个子文件夹里面有tif,
# 不然输出结果就会混在一起了，我就需要一个，所以就不单独设置了，输出路径就直接按照输入路径里面的年份进行创建了
# Z:\mcd43a3_mrt_rep\2000\reproject
# Z:\mcd43a3_mrt_rep\2001\reproject
# Z:\mcd43a3_mrt_rep\2002\reproject
#
# 输入父目录，包含多个年份的子文件夹，每个子文件夹下有'reproject'文件夹存放待处理栅格
input_parent_dir = r"Z:\\mcd43a3_mrt_rep"# 输出根目录，重采样结果会保存到对应年份的子文件夹内，保持文件夹结构一致
output_base_root = r"Z:\\mcd43a3_mrt_rep_res0.1"# 重采样方法，支持："NEAREST"（最邻近插值）、"BILINEAR"（双线性插值）、"CUBIC"（三次卷积插值）
resample_method = "NEAREST"# 输出像元大小，格式为 "X_size Y_size" ，这里示例是0.1度*0.1度 
# 注意，本身要有经纬度，要有参考坐标，不然可能会报错，未进行尝试，所以在有空间参考下，可执行
output_cell_size = "0.1 0.1"# 日志文件名，程序运行时会将日志内容写入此文件（UTF-8编码）
log_file = "Z:\\mcd43a3_mrt_rep\\resample_log.txt"# ==============================
# 控制台和日志双输出函数
# ==============================
def log_print(text):"""负责打印日志信息，双重输出：1. 控制台打印（尝试用GBK编码避免Windows cmd乱码）2. 写入UTF-8编码的日志文件"""try:# 打印到控制台，使用GBK编码忽略无法编码字符，防止乱码报错print(text.encode('gbk', 'ignore').decode('gbk'))except Exception:# 出错时直接打印原文本print(text)try:# 写入日志文件，追加模式，保证多次调用都写入文件末尾with open(log_file, "a", encoding="utf-8") as f:f.write(text + "\n")except:# 忽略写文件时的异常，保证程序不中断pass# ==============================
# 核心功能：重采样栅格函数
# ==============================
def resample_raster(in_raster, out_raster, cell_size, method):"""使用 arcpy 的 Resample_management 工具对单个栅格文件进行重采样。参数:in_raster: 输入栅格文件路径（字符串）out_raster: 输出栅格文件路径（字符串）cell_size: 输出栅格的像元大小，如 "0.1 0.1"（字符串）method: 重采样方法，如 "NEAREST"（字符串）返回:tuple(bool, str or None): 成功返回 (True, None)，失败返回 (False, 错误信息)"""try:arcpy.Resample_management(in_raster, out_raster, cell_size, method)return True, Noneexcept Exception as e:return False, str(e)# ==============================
# 批量处理主函数
# ==============================
def batch_resample():"""遍历 input_parent_dir 目录下所有年份子目录的 'reproject' 文件夹，对其中所有 .tif 文件执行重采样操作，并将结果按年份和目录结构输出到 output_base_root 下，同时记录日志信息。"""# 如果日志文件存在，先删除，保证每次运行日志清晰if os.path.exists(log_file):os.remove(log_file)# 1. 自动扫描所有年份子文件夹的'reproject'路径，加入待处理列表input_root_years = []for year_folder in os.listdir(input_parent_dir):reproj_path = os.path.join(input_parent_dir, year_folder, "reproject")if os.path.isdir(reproj_path):input_root_years.append(reproj_path)# 2. 遍历每个年份的'reproject'目录，递归处理所有子目录for root_input_folder in input_root_years:# 根据路径解析出年份，例如 Z:\mcd43a3_mrt_rep\2000\reproject -> year = "2000"year = os.path.basename(os.path.dirname(root_input_folder))# 构造对应的输出根目录，保证输出目录中包含年份文件夹root_output_folder = os.path.join(output_base_root, year)# 递归遍历所有子文件夹，保持结构一致for foldername, subfolders, filenames in os.walk(root_input_folder):# 相对路径，相对于当前年份的reproject文件夹rel_path = os.path.relpath(foldername, root_input_folder)# 结合输出根目录，构造输出子目录路径output_folder = os.path.join(root_output_folder, rel_path)# 设置arcpy工作空间为当前文件夹，方便调用ListRastersarcpy.env.workspace = foldername# 获取当前目录所有.tif栅格文件rasters = arcpy.ListRasters("*.tif")if not rasters:continue  # 当前目录无tif文件则跳过total = len(rasters)log_print("\nFolder: {}\nStart resampling, total: {} files\n".format(foldername, total))# 如果输出目录不存在，则创建if not os.path.exists(output_folder):os.makedirs(output_folder)# 对每个栅格文件执行重采样for idx, raster in enumerate(rasters):input_path = os.path.join(foldername, raster)output_path = os.path.join(output_folder, raster)  # 输出保持原文件名log_print("[{}/{}] Processing: {} ...".format(idx + 1, total, input_path))start_time = time.time()# 调用重采样函数success, error = resample_raster(input_path, output_path, output_cell_size, resample_method)# 根据结果输出日志if success:log_print("    -> Finished, elapsed time {:.2f} sec\n".format(time.time() - start_time))else:log_print("    -> Failed: {}\n".format(error))# ==============================
# 主程序入口
# ==============================
if __name__ == '__main__':batch_resample()

这里的结果看不出来，只是贴个图，看结果可以直接使用ArcMap查看属性即可。

在运行程序之前确保：重采样工具需要 ArcGIS Spatial Analyst 扩展许可；