ArcGIS栅格图层批量归一化：原理剖析与Python自动化实践

1次阅读

共计 2829 个字符，预计需要花费 8 分钟才能阅读完成。

在遥感分析和环境监测中，栅格数据归一化是常见预处理步骤。传统手动操作在 ArcGIS 中需要重复执行以下流程：

打开每个栅格文件
手动计算最大值最小值
应用栅格计算器公式
保存结果

这种操作方式存在三个明显缺陷：

时间成本高 ：处理 100 个文件可能需要整天时间
人为错误风险 ：公式输入错误可能导致整批数据作废
可复现性差 ：难以保证每次处理的参数完全一致

特性	GDAL	ArcPy
处理速度	较快	中等
功能完整性	基础功能	完整 ArcGIS 功能集
学习曲线	陡峭	平缓
许可证要求	开源	需 ArcGIS 许可

选择 ArcPy 的核心优势在于：

直接集成 ArcGIS 空间分析功能
兼容 Esri 格式的工作空间
可利用 arcpy.sa 模块的栅格代数运算

import arcpy
from arcpy.sa import *
import os
import time

# 检查空间分析扩展许可
arcpy.CheckOutExtension("Spatial")

def normalize_raster(input_raster, output_folder):
    """
    执行归一化：(值 - 最小值)/(最大值 - 最小值)
    :param input_raster: 输入栅格路径
    :param output_folder: 输出目录
    """
    try:
        # 获取输入文件名（无扩展名）basename = os.path.basename(input_raster).split('.')[0]

        # 计算统计值（忽略 NoData）min_val = arcpy.GetRasterProperties_management(input_raster, "MINIMUM")
        max_val = arcpy.GetRasterProperties_management(input_raster, "MAXIMUM")

        # 执行栅格代数运算
        out_raster = (Raster(input_raster) - float(min_val.getOutput(0))) / \
                     (float(max_val.getOutput(0)) - float(min_val.getOutput(0)))

        # 保存结果
        output_path = os.path.join(output_folder, f"{basename}_normalized.tif")
        out_raster.save(output_path)

        return True
    except Exception as e:
        print(f"Error processing {input_raster}: {str(e)}")
        return False

def batch_normalize(input_workspace, output_folder):
    """
    批量处理工作空间内的所有 TIFF 文件
    :param input_workspace: 输入工作空间（文件夹或地理数据库）:param output_folder: 输出目录
    """
    start_time = time.time()
    processed_count = 0

    # 创建输出目录（如果不存在）if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    # 遍历工作空间
    raster_list = arcpy.ListRasters("*.tif")
    print(f"Found {len(raster_list)} TIFF files to process")

    for raster in raster_list:
        if normalize_raster(raster, output_folder):
            processed_count += 1

    # 释放空间分析许可
    arcpy.CheckInExtension("Spatial")

    print(f"Processed {processed_count}/{len(raster_list)} files in"
          f"{time.time()-start_time:.2f} seconds")

分块处理 ：对大文件使用 arcpy.sa.TileIterator

tiles = arcpy.sa.TileIterator(input_raster, tile_size="512 512")
for tile in tiles:
    # 对每个分块执行计算

使用 NumPy 加速 ：

arr = arcpy.RasterToNumPyArray(input_raster)
norm_arr = (arr - np.min(arr)) / (np.max(arr) - np.min(arr))
result = arcpy.NumPyArrayToRaster(norm_arr)

由于 ArcGIS 的许可证限制，建议采用以下方案：

将文件列表拆分为多个子集
为每个子集创建独立的 Python 进程
每个进程单独获取 ArcGIS 许可

解决方法：在处理前统一投影

# 获取目标坐标系（以第一个文件为准）ref_raster = raster_list[0]
spatial_ref = arcpy.Describe(ref_raster).spatialReference

# 对输出结果应用统一坐标系
arcpy.DefineProjection_management(output_path, spatial_ref)

关键代码修改：

# 在归一化计算前排除 NoData
out_raster = Con(IsNull(input_raster), 0, 
                (Raster(input_raster) - min_val) / (max_val - min_val))

记录已处理文件：

with open("processed.log", "a") as f:
    f.write(f"{input_raster}\n")

中断后恢复：

# 读取日志文件跳过已处理
if os.path.exists("processed.log"):
    with open("processed.log") as f:
        processed = set(line.strip() for line in f)
    raster_list = [r for r in raster_list if r not in processed]

本方案可轻松扩展至其他栅格运算场景：