【Arcpy】与边界成片相连要素的筛选脚本

有的时候我们想要筛选出与某个边界成片相连的所有要素，但在实际操作中，由于内部要素并未与边界直接相接触，而是通过与其他要素相接的方式与边界间接相连。如图中，要素B和C与边界并不相连，但要素A、B、C互相连接成片，而A与边界相连，如果想要简单使用”选择——按位置选择“中的“接触源图层要素边界（BOUNDARY_TOUCHES）”两条要素来选择要素，则只能选中要素A，而要素B和C则会被遗漏。

注意：BOUNDARY_TOUCHES将不会识别相交的情况！

想要同时选中A、B、C应该使用一种“生长”的思路，首先按照位置选择和边界直接相接的要素一，接着按照位置选择与要素一相接的要素二，然后按照位置选择与要素二相接的要素三，不停地循环前面的过程，直到某一个选择步骤不再增加新的选中要素。伪代码如下：

selectConnectToBorder(target_fc, border_fc)

	fc_now <- target_fc.touch(border_fc)
	num_pre <- 0
	num_now <- number(fc_pre)
	
	while num_now ≠ num_pre do
		fc_now <- target_fc.touch(fc_now)
		num_pre <- num_now
		num_now <- number(fc_now)
	end

	save(fc_now)

end

具体实现代码如下：

# -*- coding: utf-8 -*-
import os
import arcpy

arcpy.env.overwriteOutput = True

def selectConnectToBorder(target_fc, border_fc):

    # 结果默认保存在与target相同gdb中
    target_gdb = os.path.split(target_fc)[0]
    # temp要素类用于保存过程临时要素类
    temp_fc_path = os.path.join(target_gdb, "temp")
    # temp要素类复制target要素类
    arcpy.Copy_management(target_fc, temp_fc_path)
    # 删除temp要素类中所有的要素
    arcpy.DeleteFeatures_management(temp_fc_path)
    # 逐行插入temp要素类
    temp_fc_insert = arcpy.InsertCursor(temp_fc_path)

    # 按照与边界接触关系选中target要素类的要素
    arcpy.MakeFeatureLayer_management(target_fc, "target_lyr")
    arcpy.SelectLayerByLocation_management("target_lyr", "BOUNDARY_TOUCHES", border_fc)
    target_shapes = arcpy.da.SearchCursor("target_lyr", ("SHAPE@", "gridcode"))
    
    # 遍历选中的要素，插入到temp要素类中
    for target_shape in target_shapes:
        target_geometry = target_shape[0]
        target_gridcode = target_shape[1]
        temp_row = temp_fc_insert.newRow()
        temp_row.shape = target_geometry
        temp_row.setValue("gridcode", target_gridcode)
        temp_fc_insert.insertRow(temp_row)
    
    iteration_count = 1
    num_pre = 0
    num_now = int(arcpy.GetCount_management("target_lyr").getOutput(0))

    # 如果下次按位置选中的要素个数与上次不一致，说明生长过程未完成继续循环
    while num_pre != num_now:
        arcpy.SelectLayerByLocation_management("target_lyr", "BOUNDARY_TOUCHES", temp_fc_path)
        arcpy.DeleteFeatures_management(temp_fc_path)

        temp_fc_insert = arcpy.InsertCursor(temp_fc_path)
        target_shapes = arcpy.da.SearchCursor('target_lyr', ("SHAPE@", "gridcode"))

        for target_shape in target_shapes:
            target_geometry = target_shape[0]
            target_gridcode = target_shape[1]
            temp_row = temp_fc_insert.newRow()
            temp_row.shape = target_geometry
            temp_row.setValue("gridcode", target_gridcode)
            temp_fc_insert.insertRow(temp_row)

        num_pre = num_now
        num_now = int(arcpy.GetCount_management("target_lyr").getOutput(0))
        iteration_count += 1

    # 搜寻完毕，保存结果
    result_fc_path = os.path.join(target_gdb, "result")
    arcpy.Copy_management(temp_fc_path, result_fc_path)
    arcpy.DeleteFeatures_management(temp_fc_path)


if __name__ == "__main__":
	target_fc = unicode("C:/Users/yenma/Desktop/test.gdb/target", "utf-8")
	border_fc = unicode("C:/Users/yenma/Desktop/test.gdb/border", "utf-8")
	selectConnectToBorder(target_fc, border_fc)