利用Arcpy筛选孤立图斑并批量更改其属性

生态系统分类在自动化分类过程中“针叶”和“阔叶”可能会出现错分的情况，Arcgis自带的栅格计算器和按属性选择具有简单易操作等特点，但是对于广域面积范围，多条件，批量筛选和更改等操作显得力不从心。利用Arcgis自带的python库（Arcpy）能够实现批量筛选和更改。

本文以某地区森林生态系统分类为案例，借助于Arcgis Pro3.0.2提供的Python3.9对分类中出现的孤立图斑进行修改。

---

一、数据源

本文使用的数据源包含2020年自动化分类后的图斑和2000年自动化分类的图斑（分辨率：90m）。

---

二、数据预处理

1. 配准

以2020年分类后的图斑为基准，将2000年分类后的影像配准到2020影像空间。

2. 做差

对2000年影像和2020年影像进行做差处理，得到差值后的影像。

3. 重分类

对做差后的影像进行重分类，分为变化影像和未变化影像。

4. 栅格转矢量

分别对差值后的影像进行栅格转矢量操作，用以数据处理。

---

三、数据处理

1. 整体思路

① 设置面积阈值，筛选出带有“针叶”和“阔叶”的孤立图斑。

② 计算孤立图斑周围一定范围内面积最大并带有“针叶”或者“阔叶”的图斑并筛选出来。

③ 将筛选出来面积最大图斑的分类属性赋值给孤立图斑，即达到了更改孤立图斑分类属性的效果。

整体的流程图如下：

2. isolated_change算法实现

import arcpy

# 设置环境参数
arcpy.env.workspace = "./workspace"
arcpy.env.overwriteOutput = True

# 加载Shapefile文件
input_shp_changed = "./result/changed_result.shp"
input_shp_unchanged = "./result/unchanged_result.shp"

# 创建一个文本文件用于输出过程变量
file = "./isolated_maxarea.txt"

# 获取changed图斑的特定字段
# 若获取所有字段，用["*"]
# "AREA < 10000"为SQL语句，筛选出面积小于10000的图斑
features_changed = arcpy.da.UpdateCursor(input_shp_changed, ["OID@", "SHAPE@", "CLASS", "AREA"], "AREA < 10000")

# 遍历每个图斑
with open(file, 'w', encoding = "utf-8") as f:
    for feature_changed in features_changed:
        oid_changed = feature_changed[0]
        shape_changed = feature_changed[1]
        classify_changed = feature_changed[2]
        
        # 筛选出字段里带有“阔叶”或者“针叶”的图斑
        if "阔叶" in classify_changed or "针叶" in classify_changed:

            # 要用到的arcpy.SelectLayerByLocation_management接口需要输入图层文件，将要查询的文件转换成图层
            arcpy.MakeFeatureLayer_management(input_shp_unchanged, "unchanged_lyr")
            
            # 设置查找的距离为500米
            adjacent_features = arcpy.SelectLayerByLocation_management("unchanged_lyr", "WITHIN_A_DISTANCE", shape_changed, "500 Meters")

            # 找出500米范围内面积最大的图斑
            # 一般这里使用max_area = 0, max_feature = None来初始化变量，但可能会出现"Int, None Type"不可迭代的情况，所以这里使用孤立图斑本身来初始化
            max_area = feature_changed[3]
            max_feature = feature_changed[0]
            max_classify = feature_changed[2]
            adjacent_features_cursor = arcpy.da.SearchCursor(adjacent_features, ["OID@", "AREA", "CLASS"])
            for adjacent_feature in adjacent_features_cursor:
                adjacent_oid = adjacent_feature[0]
                adjacent_area = adjacent_feature[1]
                adjacent_classify = adjacent_feature[2]
                if adjacent_area > max_area:
                    max_area = adjacent_area
                    max_feature = adjacent_oid
                    max_classify = adjacent_classify
            
            # 利用强制换行导出孤立图斑及其对应面积最大的图斑，用以检验代码运行的结果
            f.write(f"{oid_changed, classify_changed, max_feature, max_classify, max_area}\n")

            # 对孤立图斑进行更改，注意这里的缩进要在上面的if "阔叶" ~~~之后，这样才会对带有“阔叶”和“针叶”的图斑进行更改
            if "阔叶" in max_classify or "针叶" in max_classify:
                feature_changed[2] = max_classify
                features_changed.updateRow(feature_changed)

---

四、代码运行结果

共找到70073个孤立图斑及其对应的面积最大的图斑：

 打开更改前和更改后的矢量文件，结果如下：

与生成的txt文件对应，更改成功！

---

五、代码优化

此代码总共检索50多万个图斑，最终对7万多个图斑进行更改，总共耗时3个多小时，笔者对该代码进行更改，但使用了更多的循环和嵌套，运行速度远不如之前的代码，有感兴趣的可以对其进行优化（将找出临近图斑的最大面积封装成一个函数也是个不错的想法），附上代码如下：

import arcpy

arcpy.env.workspace = "./workspace"
arcpy.env.overwriteOutput = True

input_shp_changed = "./result/changed_result.shp"
input_shp_unchanged = "./result/unchanged_result.shp"

file = "./temp.txt"

with open(file, "w", encoding="utf-8") as f:

    with arcpy.da.SearchCursor(input_shp_changed, ["OID@", "SHAPE@", "CLASS", "AREA"], "AREA < 10000") as changed_cursor:
        for changed_row in changed_cursor:
            if "阔叶" in changed_row[2] or "针叶" in changed_row[2]:
                arcpy.MakeFeatureLayer_management(input_shp_unchanged, "unchanged_lyr")
                adjacent_features = arcpy.SelectLayerByLocation_management("unchanged_lyr", "WITHIN_A_DISTANCE", changed_row[1], "500 Meters")
                max_area = changed_row[3]
                max_feature = changed_row[0]
                max_classify = changed_row[2]
                adjacent_features_cursor = arcpy.da.SearchCursor(adjacent_features, ["OID@", "CLASS", "AREA"])
                for adjacent_feature in adjacent_features_cursor:
                    adjacent_oid = adjacent_feature[0]
                    adjacent_area = adjacent_feature[2]
                    adjacent_classify = adjacent_feature[1]
                    if adjacent_feature[2] > max_area:
                        max_area = adjacent_feature[2]
                        max_feature = adjacent_feature[0]
                        max_classify = adjacent_feature[1]
                f.write(f"{changed_row[0], changed_row[2], adjacent_feature[0], adjacent_feature[1], adjacent_feature[2]}\n")

                if "阔叶" in max_classify or "针叶" in max_classify:
                    feature_changed[2] = max_classify
                    features_changed.updateRow(feature_changed)