如何使用Python/arcpy计算多幅大尺寸栅格的像元平均值

• 你当前写的逐像元双层循环写法效率极低，哪怕内存足够，跑全球尺度栅格也要花数倍于正常方法的时间，实现逻辑完全没必要这么复杂
• 32位Python内存上限只有2G左右，处理全球MODIS栅格必然触发内存溢出，硬凑64位后台处理因为权限问题跑不通属于ArcGIS桌面版的常见坑，不用死磕这个环境

方法1：直接调用ArcPy空间分析工具（最简便，零内存压力） #

ArcGIS的CellStatistics工具原生支持NoData值跳过，不需要自己写循环转数组处理，工具会自动按块读取栅格计算，不会把所有数据一次性加载进内存，对大栅格兼容性极好。
直接用下面的代码替换原有逻辑即可：

import arcpy, os, glob
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True

# 路径配置
inws = "G:/data0610/MODIS_VI/EVI/EVI_pro/"
outws = "G:/data0610/MODIS_VI/EVI/EVI_pro/"
ref_raster = "G:/data0610/MODIS_VI/sm_mean.tif" # 参考栅格，用于对齐坐标系、分辨率、范围

# 读取所有输入tif
rasters = glob.glob(os.path.join(inws, "*.tif"))

# 环境预配置：提前设置坐标系、范围、像元大小、捕捉栅格，输出自动和参考栅格对齐，无需后续单独重采样、重投影
arcpy.env.outputCoordinateSystem = Raster(ref_raster)
arcpy.env.extent = Raster(ref_raster)
arcpy.env.cellSize = Raster(ref_raster)
arcpy.env.snapRaster = Raster(ref_raster)

# 预处理：把MODIS小于0的填充值设为NoData，避免无效值参与计算
valid_rasters = []
for ras_path in rasters:
    ras = Raster(ras_path)
    valid_ras = SetNull(ras < 0, ras)
    valid_rasters.append(valid_ras)

# 计算像元均值，参数"DATA"表示计算时自动跳过NoData值
mean_ras = CellStatistics(valid_rasters, "MEAN", "DATA")

# 直接输出结果
out_path = os.path.join(outws, "evi_mean.tif")
mean_ras.save(out_path)
print("均值栅格计算完成，输出路径：", out_path)

注意：如果你的MOD13C2产品原始填充值为-3000这类特殊值，只需要调整SetNull里的判断条件即可，和原有逻辑里判断有效值≥0的规则保持一致。

方法2：numpy分块计算方案（无空间分析许可时使用） #

如果需要用numpy处理，不要把整幅全球栅格一次性转成数组加载，也不要写逐像元的Python双层循环（Python循环比numpy向量化运算慢几十上百倍），用分块读取+向量化计算的方式编写，哪怕32位环境也能稳定运行：

import arcpy, os, glob, numpy
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True

inws = "G:/data0610/MODIS_VI/EVI/EVI_pro/"
outws = "G:/data0610/MODIS_VI/EVI/EVI_pro/"
ref_raster = "G:/data0610/MODIS_VI/sm_mean.tif"
rasters = glob.glob(os.path.join(inws, "*.tif"))

# 读取参考栅格基础信息
r_ref = Raster(ref_raster)
lowerLeft = arcpy.Point(r_ref.extent.XMin, r_ref.extent.YMin)
cellW = r_ref.meanCellWidth
cellH = r_ref.meanCellHeight
no_data_val = r_ref.noDataValue
rows_total = r_ref.height
cols_total = r_ref.width

# 分块大小设置，每次读取1024*1024像素块，内存占用可控
block_size = 1024
out_path = os.path.join(outws, "evi_mean_block.tif")

# 逐块计算
for i in range(0, rows_total, block_size):
    for j in range(0, cols_total, block_size):
        # 计算当前块的实际行列数（边缘块可能不足设定的block_size）
        cur_rows = min(block_size, rows_total - i)
        cur_cols = min(block_size, cols_total - j)
        # 计算当前块的左下角坐标
        block_xmin = lowerLeft.X + j * cellW
        block_ymin = lowerLeft.Y + i * cellH
        block_corner = arcpy.Point(block_xmin, block_ymin)

        sum_arr = numpy.zeros((cur_rows, cur_cols), dtype=numpy.float64)
        count_arr = numpy.zeros((cur_rows, cur_cols), dtype=numpy.int32)

        # 读取所有栅格在当前块的数值
        for ras_path in rasters:
            ras = Raster(ras_path)
            block_arr = arcpy.RasterToNumPyArray(
                ras, block_corner, cur_cols, cur_rows, no_data_val
            )
            # 向量化判断有效值，无需逐像素循环
            valid_mask = block_arr >= 0
            sum_arr[valid_mask] += block_arr[valid_mask]
            count_arr[valid_mask] += 1
        
        # 计算当前块均值，无效位置设为NoData
        mean_arr = numpy.full((cur_rows, cur_cols), no_data_val, dtype=numpy.float32)
        valid_pixel = count_arr > 0
        mean_arr[valid_pixel] = sum_arr[valid_pixel] / count_arr[valid_pixel]

        # 将当前块写入输出栅格
        block_raster = arcpy.NumPyArrayToRaster(
            mean_arr, block_corner, cellW, cellH, no_data_val
        )
        if i ==0 and j ==0:
            arcpy.DefineProjection_management(block_raster, r_ref.spatialReference)
            block_raster.save(out_path)
        else:
            arcpy.Mosaic_management(block_raster, out_path, "LAST", "MATCH")
        print(f"完成块({i},{j})计算")

print("分块计算完成，输出路径：", out_path)

原有代码的核心问题说明 #

• 逐像元编写Python双层for循环效率极低，全球MODIS栅格单幅就有上千万像素，251幅的循环次数达天量，跑数天都未必能跑完
• 一次性把整幅栅格加载进numpy数组，32位Python内存上限不足必然报错
• 无需手动编写重投影、重采样步骤，提前在arcpy环境里设置好参考栅格的坐标系、像元大小、捕捉栅格，工具输出会自动对齐，省掉冗余处理步骤

内容的提问来源于stack exchange，提问作者Lorraine