ArcPy批量计算栅格面积总和关联日期存储报错求助

ArcPy 解决方案 #

核心问题排查与修复 #

1. 日期提取逻辑优化
   针对两种文件名格式，用正则表达式精准提取8位日期字符串，无需区分格式位置：

   import re
   def extract_date(filename):
       match = re.search(r'\d{8}', filename)
       return match.group() if match else None
   

2. 空表与SystemError修复要点

   • 确认栅格属性表的面积字段名（常见为AREA或Shape_Area，需与实际字段名一致）
   • 处理大数量文件时添加异常捕获，避免单个文件报错中断全流程
   • 确保输出表字段类型正确（日期字段设为文本/日期型，面积字段设为浮点型）
   • 检查文件路径权限，避免因访问限制导致栅格无法读取

完整ArcPy代码示例 #

import arcpy
import os
import re

# 配置参数
raster_folder = r"C:\your_raster_folder"
output_table = r"C:\output\area_summary.dbf"
area_field = "AREA"  # 替换为实际面积字段名

# 创建输出表
if arcpy.Exists(output_table):
    arcpy.Delete_management(output_table)
arcpy.CreateTable_management(os.path.dirname(output_table), os.path.basename(output_table))
arcpy.AddField_management(output_table, "Date", "TEXT", field_length=8)
arcpy.AddField_management(output_table, "Total_Area", "DOUBLE")

# 获取所有栅格文件
raster_list = arcpy.ListRasters("*", "TIF")  # 替换为你的栅格格式

# 循环处理每个栅格
with arcpy.da.InsertCursor(output_table, ["Date", "Total_Area"]) as cursor:
    for raster in raster_list:
        raster_path = os.path.join(raster_folder, raster)
        try:
            # 提取日期
            date_str = extract_date(raster)
            if not date_str:
                print(f"跳过无日期的文件: {raster}")
                continue
            
            # 计算面积总和
            total_area = 0
            with arcpy.da.SearchCursor(raster_path, [area_field]) as scursor:
                for row in scursor:
                    total_area += row[0]
            
            # 插入数据
            cursor.insertRow([date_str, total_area])
            print(f"已处理: {raster}")
        except Exception as e:
            print(f"处理失败 {raster}: {str(e)}")
            continue

print("处理完成")

QGIS 解决方案 #

PyQGIS 批量处理代码 #

针对大数量栅格，使用PyQGIS可高效完成任务，无需依赖图形界面：

from qgis.core import (
    QgsApplication, QgsRasterLayer, QgsVectorLayer,
    QgsField, QgsFeature, QgsProject
)
from PyQt5.QtCore import QVariant
import os
import re

# 初始化QGIS应用
QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.28", True)  # 替换为你的QGIS安装路径
qgs = QgsApplication([], False)
qgs.initQgis()

# 配置参数
raster_folder = "C:/your_raster_folder"
output_table = "C:/output/area_summary.gpkg"
area_field = "AREA"  # 替换为实际面积字段名

# 创建输出矢量表（GeoPackage格式）
if os.path.exists(output_table):
    os.remove(output_table)
fields = [
    QgsField("Date", QVariant.String, "", 8),
    QgsField("Total_Area", QVariant.Double)
]
writer = QgsVectorFileWriter(
    output_table, "UTF-8", fields, QgsWkbTypes.NoGeometry,
    QgsCoordinateReferenceSystem(), "GPKG"
)

# 定义日期提取函数
def extract_date(filename):
    match = re.search(r'\d{8}', filename)
    return match.group() if match else None

# 遍历栅格文件
for root, dirs, files in os.walk(raster_folder):
    for file in files:
        if file.endswith(".tif"):  # 替换为你的栅格格式
            raster_path = os.path.join(root, file)
            try:
                # 加载栅格
                raster_layer = QgsRasterLayer(raster_path, "temp_raster")
                if not raster_layer.isValid():
                    print(f"无法加载栅格: {file}")
                    continue
                
                # 提取日期
                date_str = extract_date(file)
                if not date_str:
                    print(f"跳过无日期的文件: {file}")
                    continue
                
                # 读取属性表并计算面积总和
                total_area = 0
                provider = raster_layer.dataProvider()
                attribute_table = provider.attributeTable()
                for row in range(attribute_table.rowCount()):
                    area_val = attribute_table.value(row, area_field)
                    total_area += float(area_val) if area_val else 0
                
                # 写入结果
                feat = QgsFeature()
                feat.setAttributes([date_str, total_area])
                writer.addFeature(feat)
                print(f"已处理: {file}")
            except Exception as e:
                print(f"处理失败 {file}: {str(e)}")
                continue

# 清理资源
del writer
qgs.exitQgis()
print("处理完成")

图形界面快速处理（小批量推荐） #

1. 打开QGIS，加载所有栅格文件
2. 打开Processing Toolbox，找到Vector analysis > Summary statistics工具
3. 开启Batch Processing，将所有栅格的属性表作为输入，设置统计字段为面积字段，统计类型选Sum
4. 处理完成后，使用Field Calculator在每个统计结果表中提取日期字段
5. 使用Merge Vector Layers工具将所有结果表合并为一张总表

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