如何在ArcGIS中统计DEM栅格中高程大于150的像元数量？

Question

在ArcGIS中统计DEM栅格中高程大于150的像元数量时，常见问题为：**执行“Con”或“Raster Calculator”后得到的结果栅格无属性表，无法直接查看像元计数，或输出为浮点型栅格导致“Get Raster Properties”工具返回空值/错误**。根本原因在于——ArcGIS默认不为浮点型栅格（如多数DEM）自动生成属性表（Attribute Table），而像元计数需依赖属性表中的COUNT字段或通过“Zonal Statistics as Table”“Get Raster Properties”等工具获取；若未先将条件结果转为整型（如用Int(Con("dem" > 150, 1))），后续统计将失败。此外，忽略像元大小、投影坐标系与地理坐标系差异，可能导致面积换算偏差；未检查NoData值处理逻辑（如Con函数中未显式定义else分支），亦会遗漏或误计无效像元。该问题高频出现在ArcGIS Pro 3.x与ArcMap 10.8的批量处理脚本及模型构建器中。

Answer 1

```html

一、现象层：典型报错与界面异常表现

• 执行 Raster Calculator 输入 "dem" > 150 后，输出栅格右键 → “属性” → “源”选项卡中显示 “无属性表”；
• 调用 Get Raster Properties 工具并设置 Property = COUNT 时，返回空字符串或 ERROR 000864；
• 在 ArcGIS Pro 3.1 的 Python 窗口中运行 arcpy.GetRasterProperties_management(raster, "COUNT") 抛出 RuntimeError: ERROR 000864；
• 模型构建器中连接 Con 输出至 Zonal Statistics as Table 时，提示“输入栅格不包含有效区域”；
• 导出为 TIFF 后用 QGIS 打开，发现像元值为浮点型（如 1.0、0.0），且无统计直方图支持。

二、机制层：ArcGIS 栅格属性表生成的底层逻辑

ArcGIS 仅对 整型（Integer）栅格 自动构建属性表（Attribute Table），其本质是基于唯一像元值建立 VALUE 与 COUNT 的关系表。而 DEM 多为 Float32，Con("dem" > 150, 1) 默认继承输入栅格数据类型，输出仍为浮点型——此时 COUNT 字段无法生成，所有依赖该字段的工具均失效。

栅格类型	是否自动生成属性表	COUNT 可读性	GetRasterProperties 支持
Float32（如原始 DEM）	❌ 否	❌ 不可用	❌ 报错
Int16（如 Int(Con(...))）	✅ 是（需手动构建或首次访问触发）	✅ VALUE/COUNT 可见	✅ 返回整数计数

三、工程层：健壮性统计流程设计（含坐标系与 NoData 控制）

1. 强制类型转换：使用 Int(Con("dem" > 150, 1, 0)) —— 显式定义 else 分支为 0，避免 NoData 扩散；
2. 投影校验：通过 arcpy.Describe(raster).spatialReference.type 判断是否为 "Projected"，若为地理坐标系（GCS），须先投影再计算面积；
3. 像元面积补偿：调用 arcpy.GetRasterProperties_management(out_int_raster, "CELLSIZEX") 获取分辨率，结合空间参考计算真实地表面积；
4. 批量容错封装：Python 脚本中嵌套 try/except 捕获 ERROR 000864 并自动触发 arcpy.BuildRasterAttributeTable_management。

四、代码层：生产级 Python 实现（ArcPy + 错误恢复）

import arcpy
from arcpy import env

env.overwriteOutput = True
dem_path = r"C:\data\dem.tif"
threshold = 150

# 步骤1：条件转整型（显式NoData处理）
int_cond = arcpy.sa.Int(arcpy.sa.Con(arcpy.Raster(dem_path) > threshold, 1, 0))
int_cond.save(r"C:\data\dem_gt150_int.tif")

# 步骤2：强制构建属性表（关键！）
arcpy.BuildRasterAttributeTable_management(int_cond, overwrite=True)

# 步骤3：安全获取COUNT
try:
    count = int(arcpy.GetRasterProperties_management(int_cond, "COUNT").getOutput(0))
    cellsize_x = float(arcpy.GetRasterProperties_management(int_cond, "CELLSIZEX").getOutput(0))
    spatial_ref = arcpy.Describe(int_cond).spatialReference
    area_km2 = count * (cellsize_x ** 2) / 1e6 if spatial_ref.type == "Projected" else "需重投影后计算"
    print(f"像元数量：{count} | 面积估算：{area_km2} km²")
except RuntimeError as e:
    print(f"属性表读取失败：{e}")

五、架构层：模型构建器与脚本协同的推荐范式

graph TD
  A[输入DEM] --> B{Con条件表达式}
  B --> C[Int转换]
  C --> D[Build Raster Attribute Table]
  D --> E[Get Raster Properties COUNT]
  D --> F[Zonal Statistics as Table]
  E --> G[Calculate Value Area]
  F --> H[Join to Zone Features]
  style B fill:#ffe4b5,stroke:#ff8c00
  style D fill:#98fb98,stroke:#2e8b57

六、陷阱层：被忽视的三大隐性风险

• NoData 传播陷阱：未写 Con(..., 1, 0) 而仅写 Con(..., 1)，导致原 DEM 中 NoData 区域在输出中仍为 NoData，被 GetRasterProperties 完全忽略，计数严重偏低；
• 地理坐标系面积失真：在 WGS84 下直接用 CELLSIZEX × CELLSIZEY × COUNT 计算面积，高纬度地区误差可达 300% 以上；
• 位深度截断风险：对 >65535 个满足条件像元使用 Byte 类型（如 Int(Con(...)) 未指定输出数据类型），引发整数溢出，部分像元值变为 0 或 255。

```