Arcpy案例:用地用海图斑行政区分割

188 阅读3分钟

需求说明

在城乡规划项目中,需要统计行政区下一级的用地用海现状数据,例如:对于镇级的项目,需要知道这些项目地理位置分布在哪些村中,为此,可以以行政区表中的村级行政区为分割字段,对用地用海项目进行分割提取。

行政区表里面的村级行政区字段

解决思路

  1. 在行政区表中对村级行政区字段进行分割提取,获得村级行政区图斑,函数:Select_analysis。参考:Select_analysis 官方文档
  2. 使用生成的村级行政区图斑对用地用海项目图斑进行裁剪操作,得到村级行政区内的用地用海项目,函数:Clip_analysis。参考:Clip_analysis 官方文档

代码实现

主体实现代码

def XZQ_divide_feature(in_feature, clip_feature, clip_code_field, clip_name_field, square_field, res_path, xy_tolerance):
    """
    :param in_feature: 输入的待分割的源要素
    :param clip_feature: 用于进行分割要素,含有行政区字段信息
    :param clip_code_field: 分割的行政区代码字段
    :param clip_name_field: 分割的行政区名称字段
    :param square_field: 面积字段
    :param res_path: 结果保存路径
    :param xy_tolerance: XY容差
    :return: 
    """
    # 用于保存分割的行政区信息,包括行政区代码和行政区名称
    xzq_table = []
    # 访问分割要素的属性表数据
    data_table = arcpy.SearchCursor(clip_feature)
    for row in data_table:
        xzq_table.append([row.getValue(clip_code_field), row.getValue(clip_name_field)])
    # 保存分割结果的文件路径
    res_paths = []
    # 记录任务总量,用于设置进度条信息
    total = len(xzq_table) + 1
    # 设置进度条
    arcpy.SetProgressor("step", "Clipping YDYH", 0, total, 1)
    for row in xzq_table:
        # 行政区代码
        area_code = row[0]
        # 行政区名称
        area_name = row[1]
        # 选择条件
        where_express = ' "{0}"=\'{1}\' '.format(clip_name_field, area_name)
        clip_layer_path = os.path.join(res_path, area_name)
        # 如果行政区名称存在,则删除
        if arcpy.Exists(clip_layer_path):
            arcpy.Delete_management(clip_layer_path)
        # 根据行政区名称字段,将下一级行政区从分割要素图层提取出来,作为分割用的中间文件
        arcpy.Select_analysis(clip_feature, clip_layer_path, where_express)
        # 分割结果路径
        res_layer_path = os.path.join(res_path, area_name + 'div')
        res_paths.append(res_layer_path)
        # 如果分割结果存在,则删除
        if arcpy.Exists(res_layer_path):
            arcpy.Delete_management(res_layer_path)
        # 使用行政区图层对源要素进行分割
        arcpy.Clip_analysis(in_feature, clip_layer_path, res_layer_path, xy_tolerance)
        # 使用字段计算器更新结果中的行政区代码、行政区名称、面积数据
        arcpy.CalculateField_management(res_layer_path, clip_code_field, area_code, "PYTHON_9.3")
        arcpy.CalculateField_management(res_layer_path, clip_name_field, u"'{0}'".format(area_name), "PYTHON_9.3")
        arcpy.CalculateField_management(res_layer_path, square_field, "!Shape.geodesicArea!", "PYTHON_9.3")
        # 删除用于分割的行政区中间文件
        arcpy.Delete_management(clip_layer_path)
        # 更新进度条状态
        arcpy.SetProgressorLabel("{0} Finished".format(area_name))
        arcpy.SetProgressorPosition()
    # 图层合并路径
    merged_layers=os.path.join(res_path, "MergedLayers")
    # 如果合并路径存在,则删除
    if arcpy.Exists(merged_layers):
        arcpy.Delete_management(merged_layers)
    # 将所有分割完成的图层合并为一个新的图层
    arcpy.Merge_management(res_paths, merged_layers)
    arcpy.SetProgressorLabel("Merge all layers")
    arcpy.SetProgressorPosition()
    arcpy.ResetProgressor()

参数获取代码

# 待分割的源要素
source_feature = arcpy.GetParameterAsText(0)
# 带有下一级行政区信息的分割要素
divide_feature = arcpy.GetParameterAsText(1)
# 分割的行政区代码
divide_code = arcpy.GetParameterAsText(2)
# 分割的行政区名称
divide_name = arcpy.GetParameterAsText(3)
# 面积字段
square = arcpy.GetParameterAsText(4)
# 输出路径
out_gdb = arcpy.GetParameterAsText(5)
# XY容差
tolerance = float(arcpy.GetParameterAsText(6))

ArcGIS工具导入

参数设置:

插件运行界面效果:

数据分割效果: