arcpy处理要素合并,要求要素面积小于5平方公里的与相邻要素合并,且在合并后的属性表中写上被合并要素的属性值

以下是该操作的代码:


import arcpy
import math

# 设置输入和输出路径
input_feature = r"C:\data.gdb\input_feature"
output_feature = r"C:\data.gdb\output_feature"
max_area = 5000 # 最大面积(单位:平方米)

# 初始化变量
merged_features = []
processed_features = set()

# 获取要素的相关属性字段
fields = [field.name for field in arcpy.ListFields(input_feature) if not field.required]

# 遍历所有要素
with arcpy.da.SearchCursor(input_feature, ["OID@", "SHAPE@"] + fields) as cursor:
    for row in cursor:
        feature_id = row[0]
        feature_shape = row[1]
        feature_area = feature_shape.getArea("PLANAR")
        if feature_area >= max_area:
            # 如果当前要素的面积超过最大面积限制,则将其直接写入输出图层
            feature_row = [feature_id, feature_shape] + list(row[3:])
            merged_features.append(feature_row)
        else:
            # 如果当前要素的面积小于最大面积限制,则查找所有与之相交的要素,进行合并
            if feature_id not in processed_features:
                intersecting_features = []
                with arcpy.da.SearchCursor(input_feature, ["OID@", "SHAPE@"] + fields,
                                           f"OID <> {feature_id} AND SHAPE.STIntersects({feature_shape})") as sub_cursor:
                    for sub_row in sub_cursor:
                        sub_feature_id = sub_row[0]
                        sub_feature_shape = sub_row[1]
                        sub_feature_area = sub_feature_shape.getArea("PLANAR")
                        if sub_feature_area < max_area:
                            # 如果相邻要素的面积也小于最大面积限制,则将其合并
                            feature_shape = feature_shape.union(sub_feature_shape)
                            processed_features.add(sub_feature_id)
                            intersecting_features.append((sub_feature_id, sub_row))
                    feature_row = [feature_id, feature_shape] + list(row[3:])
                    merged_features.append(feature_row)
                    for (intersecting_id, intersecting_row) in intersecting_features:
                        merged_row = [feature_id, feature_shape] + list(intersecting_row[3:])
                        merged_features.append(merged_row)

# 将所有合并后的要素写入输出图层
arcpy.CreateFeatureclass_management("C:/data.gdb", "output_feature", "POLYGON", spatial_reference=input_feature)

with arcpy.da.InsertCursor(output_feature, ["SHAPE@"] + fields) as cursor:
    for row in merged_features:
        feature_shape = row[1]
        feature_fields = row[2:]
        if feature_shape is not None:
            cursor.insertRow([feature_shape] + feature_fields)

print(f"已将{len(merged_features)}个要素合并为{arcpy.GetCount_management(output_feature)}个要素:{output_feature}")

上述代码会遍历指定图层(`input_feature`)的每一个要素,并判断其面积是否小于5平方公里。如果小于等于5平方公里,则查找所有与之相邻且面积也小于5平方公里的要素,并将其合并成一个新的要素,同时在新要素的属性表中写上被合并要素的属性值。如果当前要素的面积超过5平方公里,则直接写入输出图层(`output_feature`)。最后,所有合并后的要素都会写入输出图层中。

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

认真学GIS

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值