你是否，因为嫌弃而重写过ArcGIS的工具？

内容导读 （文末有完整代码）

        在对矢量数据加工处理时，尤其是基于特殊需求和批量操作时，线要素需要根据线交点或临近点，对其进行分割处理。一方面可以使用ArcGIS提供的“在点处分割线”地理处理工具；另一方面，可以按本文提供的思路和方法自己实现“在点处分割线”的工具功能。

有想法，就去试一试，万一实现了呢？ 

（1） mdb（个人地理数据库）转shape file其实并不简单

（2） mdb转gdb实现过程介绍（1）mdb地理数据库结构解析和gdb库的创建

（3）mdb转gdb实现过程介绍（2）三种方式实现GDB数据库的读、写，并将实现方式与ArcGIS环境共存配置

（4）使用Arcpy自定义曲线线段替换为线段的实现

（5）自定义Python工具箱实现mdb转出为shp或gdb格式----终章（工具免费）

一、基于arcpy实现在点处分割线功能

        小编初心是重写ArcGIS Pro中的“在点处分割线”工具，提高工具执行效率；但高版本（3.0+）的arcgis pro已经对此工具进行了优化，优化后大家半斤八两。

        但小编写的脚本工具在使用场景上更加丰富，支持源图层要素直接修改和对分割后的结果另存输出两种方式。并提供代码，与大家一同交流学习。

1.1 ArcGIS Pro“在点处分割线”工具

---

摘 要

根据交叉点或与点要素的邻近性分割线要素。

ArcGIS Pro2.8

---

ArcGIS Pro3.0+

---

        从Pro2.8版本和Pro3.0+版本的帮助文档对比，可以看出，在高版本，“在点处分割线”地理处理工具，增加了部分功能。

        通过实际测试还可以发现，在2.8版本中，140W+条线要素，与10W+条点要素做分割处理，工具运行了3天，工具都没有执行完成；而3.1版本，同样的数据，工具仅运行了34分钟。

        小编猜测主要原因是高版本Pro在执行分割线时，使用了“空间索引”执行搜索，因为小编后面介绍的工具当时就使用了空间索引技术，最后执行效率上不相上下。

1.2 小编写的“在点处分割线”工具

        实际上，当初写这个工具的时候，是无可奈何的，当时使用的Pro版本是2.8.0，我有百万级的线要素，十万级的点要素，需要在点或点附近打断线，工具运行了三天，一直没有结果……

        140W+线要素和10W+点要素，花了40分钟。

01 工具概述

支持的功能：

        （1）它功能上与“在点处分割线”工具相同；

        （2）支持原地修改，即不用输出新的要素类，直接更改到源要素图层，这对于经常数据编辑的小伙伴可能更加贴心。

02 功能流程

        工具打开界面如下图所示：

 03 工具参数介绍如下：

04 工具输出：

        当选中“是否另存输出”参数时，输出与输入线要素同名的分割后的线要素结果；当未选中“是否另存输出”参数，修复结果将直接在输入线要素图层上进行修改。

        所有输出的要素类都存放在输出目录中的scratch.gdb中。若scratch.gdb不存在，则自动创建，若已存在，不会覆盖其中已有的要素。

05 注意事项：

        工具支持是否原地修改，请按实际需要进行选择。

代码实现如下：（后来优化过该部分代码，有需要的私信）

import math
import os

import ConversionUtils
import numpy as np
from rtree import index
import arcpy
import pandas as pd
from shapely.geometry import LineString, MultiPoint, MultiLineString
from arcgis.geometry import Geometry
from shapely.ops import split, snap

arcpy.env.overwriteOutput = True


# 距离长度转测地线长度
def plane_to_geodesic(p_length, v_angle=0, lat=0):
    # cell 分辨率大小 lat 要素维度
    # 地球半径 r
    r = 6378137

    # 地球周长 R
    R = 2 * math.pi * r

    # 赤道上一度所对应的地表距离    degToMeter = math.pi * r / 180.0
    # 在经线上，变化1纬度是多少距离(对应直角坐标系y轴)
    # 北半球的纬度是0° - 90°，北半球的弧长是 R / 4 ,纬度变化1°，对应的距离 L = R / (90 * 4)
    L = R / (90 * 4)
    y = p_length / L

    # 在纬线上，变化1经度是多少距离(对应直角坐标系x轴)
    # 只要能够得到指定纬度θ所切的圆的半径 r' ，就可以得到在这个纬度θ上所切的圆的周长：R' = 2πr' ，纬度θ上变化 1经度的距离 L = R' / 360
    # 在纬度θ所切圆的周长 R' = 2πr' =  R * cosθ

    R_lat = R * math.cos(math.pi / 180 * lat)
    # 纬度θ的纬线上变化 1经度的距离 L = R' / 360

    degToMeter = R_lat / 360.0
    x = p_length / degToMeter

    degree_angle = (v_angle * np.pi) / 180

    g_length = np.sqrt((x * np.cos(degree_angle)) ** 2 + (y * np.sin(degree_angle)) ** 2)

    return g_length


# 对输入点、线、面创建空间索引，不支持多部件
def create_spatial_index(fc):
    # 创建R-tree
    p = index.Property()
    p.dimension = 2
    idx = index.Index(properties=p)
    desc = arcpy.Describe(fc)
    OID = desc.OIDFieldName

    cursor = arcpy.SearchCursor(fc)
    row = cursor.next()

    if desc.shapeType == 'Polyline':
        ss = "Geometry(row.getValue('SHAPE')).as_shapely.geoms[0]"
    elif desc.shapeType == 'Point':
        ss = "Geometry(row.getValue('SHAPE')).as_shapely"
    elif desc.shapeType == 'Polygon':
        ss = "Geometry(row.getValue('SHAPE')).as_shapely.geoms[0]"
    else:
        raise ValueError("Unsupported shape type: {}".format(desc.shapeType))
    while row:
        if ss is not None:
            sh_geo = eval(ss)
            id = row.getValue(OID)
            feature = (id, sh_geo.bounds, sh_geo)  # 包含obj
            idx.insert(*feature)
        row = cursor.next()
    return idx


# 可以先用点层和线层选一次，把无关的线层线过滤
def select_near_fc(fc1, fc2, dis=2):
    query = "{} Meters".format(dis)
    near_fc = arcpy.management.SelectLayerByLocation(fc1, "WITHIN_A_DISTANCE", fc2, query, "NEW_SELECTION",
                                                     "NOT_INVERT")
    return near_fc


# 在指定距离范围内，创建线与点的临近分析表
def create_near_table(fc_point, fc_polyline, distance=0):
    # 创建点空间索引
    point_index = create_spatial_index(fc_point)

    # 遍历每条线，并查询空间索引以获取在指定距离范围内的点

    desc = arcpy.Describe(fc_polyline)
    OID = desc.OIDFieldName
    spa_type = desc.spatialReference.type
    if spa_type == "Geographic":
        densify_dis = plane_to_geodesic(1)
        buffer_dis = 1e-6
    else:
        densify_dis = 1
        buffer_dis = 1e-3

    cursor = arcpy.SearchCursor(fc_polyline)
    row = cursor.next()
    obj_info = []
    while row:
        geo_row = row.getValue('SHAPE')
        if not geo_row:
            row = cursor.next()
        if geo_row.hasCurves:  # 几何具有曲线段
            densify_row = geo_row.densify('DISTANCE', distance=densify_dis, deviation=densify_dis / 10)
            sh_geo = Geometry(densify_row).as_shapely.geoms[0]
        else:
            sh_geo = Geometry(geo_row).as_shapely.geoms[0]

        line_id = row.getValue(OID)
        if distance > 0:  # 线向两侧缓冲
            buffered_line = sh_geo.buffer(distance, join_style=2, cap_style=2)
            box = buffered_line.bounds
            obj = buffered_line
        else:
            box = sh_geo.bounds
            obj = sh_geo

        line_result = list(point_index.intersection(box, objects=True))

        intersection = [i.object for i in line_result if obj.intersection(i.object.buffer(buffer_dis))]

        # df = df.append({'line_id': line_id, 'point_obj': intersection}, ignore_index=True)
        obj_info.append({'line_id': line_id, 'point_obj': intersection})
        row = cursor.next()
    df = pd.DataFrame(obj_info, columns=['line_id', 'point_obj'])
    return df


# 根据临近表中的点ID，获取所有ID对应的几何对象
def point_ID_geo(fc_point, df, col):
    # 提取所有列表并去除重复项
    unique_lists = set(df.explode(col)[col])
    # 对列表进行升序排序
    sorted_lists = sorted(unique_lists)
    id_geo = {}
    OID = arcpy.Describe(fc_point).OIDFieldName
    cursor = arcpy.SearchCursor(fc_point)
    row = cursor.next()
    while row:
        id = row.getValue(OID)
        if id in sorted_lists:
            geo = Geometry(row.getValue('SHAPE')).as_shapely
            id_geo[id] = geo
    return id_geo


# 获取临近表中line ID和point OBJ的映射表
def line_ID_to_point_ID(df, col1, col2):
    # 将'col1'列和'col2'列分别作为两个Series对象提取出来
    series_a = df[col1]
    series_b = df[col2]

    # 使用zip()函数将这两个Series对象合并为一个元组列表
    tuple_list = list(zip(series_a, series_b))

    # 将结果转换为字典，以便将第一列的值作为键，第二列的值作为值
    result_dict = dict(tuple_list)
    return result_dict


# 在指定点位置对线进行打断
def split_line_by_point(line, point, tolerance: float = 1.0e-6):
    if line.intersects(point.buffer(tolerance)) and line.is_valid:
        line = line.simplify(0)
        return split(snap(line, point, tolerance), point)
    else:
        return line


# 指定要素，创建空要素类模板
def create_fc_template(fc_path, out_path):
    if os.path.splitext(out_path)[-1] == '.gdb':
        gdb_path = out_path
    else:
        arcpy.env.scratchWorkspace = out_path
        gdb_path = arcpy.env.scratchGDB

    desc = arcpy.Describe(fc_path)

    catalog_path = desc.catalogPath

    name = desc.name.split(".")[0]
    spatial_ref = desc.spatialReference
    geometry_type = desc.shapeType
    inFeatureClass = os.path.join(gdb_path, name)
    if arcpy.Exists(inFeatureClass):
        outFeatureClass = ConversionUtils.GenerateOutputName(inFeatureClass, gdb_path)
        name = os.path.basename(outFeatureClass).split(".")[0]

    has_m = "DISABLED"
    has_z = "DISABLED"

    arcpy.CreateFeatureclass_management(gdb_path, name, geometry_type, catalog_path,
                                        has_m, has_z, spatial_ref)

    return os.path.join(gdb_path, name)


# 根据点打断线
def SplitLineAtPoint(fc_polyline, fc_point, out_path=None, search_distance=0.0):
    fields_name = get_fields_name(fc_polyline)
    desc = arcpy.Describe(fc_polyline)
    spatial_ref = desc.spatialReference
    spa_type = spatial_ref.type
    if spa_type == "Geographic" and search_distance > 0:
        g_radius = plane_to_geodesic(search_distance)
        densify_dis = plane_to_geodesic(1)
        buffer_dis = 1e-6
    else:
        g_radius = search_distance
        densify_dis = 1
        buffer_dis = 1e-3
    # 创建临近表
    df_table = create_near_table(fc_point, fc_polyline, distance=g_radius)

    # 获取分割线ID集合
    dict_df = line_ID_to_point_ID(df_table, "line_id", "point_obj")

    # 在点处分割线，当需要另存分割结果时，在输出目录创建输入要素模板，并插入分割结果；
    # 当原地更新输入要素时，需要构建Row对象，会比较繁琐
    if out_path is not None:
        if arcpy.Exists(out_path):
            new_fc = create_fc_template(fc_polyline, out_path)
            rows = arcpy.InsertCursor(new_fc)
        else:
            raise ValueError("输出目录不存在……")
    else:
        new_fc = fc_polyline

    OID = desc.OIDFieldName
    cursor = arcpy.UpdateCursor(fc_polyline)
    row = cursor.next()
    split_rows = []
    while row:
        points = []
        row_id = row.getValue(OID)

        geo_row = row.getValue('SHAPE')
        if geo_row.hasCurves:  # 几何具有曲线段
            densify_row = geo_row.densify('DISTANCE', distance=densify_dis, deviation=densify_dis / 10)
            sh_geo = Geometry(densify_row).as_shapely.geoms[0]
        else:
            sh_geo = Geometry(geo_row).as_shapely.geoms[0]

        point_obj = dict_df.get(row_id, [])
        if len(point_obj) >= 1:
            for point_geometry in point_obj:  # 当存在搜索距离时，点可能不在线上，需要使用该点找出离线上最近的点，使用最近的点打断线
                if search_distance > 0 and not point_geometry.intersects(sh_geo):  # 判断点是否与线相交
                    point_geometry = sh_geo.interpolate(sh_geo.project(point_geometry))
                points.append(point_geometry)

            multi_point = MultiPoint(points)

            result = split_line_by_point(sh_geo, multi_point, tolerance=buffer_dis)
            if not result:
                row = cursor.next()
                continue

            if isinstance(result, LineString):
                result = MultiLineString([result])
            if len(result.geoms) >= 2:
                # import copy
                for res in result.geoms:
                    geo_coords = Geometry.from_shapely(res)['paths'][0]
                    new_part = arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in geo_coords]),
                                              spatial_reference=spatial_ref)
                    row.setValue("SHAPE", new_part)

                    if out_path is not None:
                        rows.insertRow(row)

                    else:
                        # 使用自定义的复制row对象的类实现row的复制
                        copy_row = MyRow(row, fields_name).copy_row
                        split_rows.append(copy_row)
                        cursor.updateRow(row)
                if out_path is None:
                    split_rows.pop()

            elif out_path is not None:
                rows.insertRow(row)

        elif out_path is not None:
            # copy_row = MyRow(row, fields_name).copy_row
            # split_rows.append(copy_row)

            rows.insertRow(row)

        row = cursor.next()
    del row
    del rows

    # 插入分割出的结果
    if len(split_rows) > 0 and out_path is None:
        rows = arcpy.InsertCursor(new_fc)
        for i in range(len(split_rows)):
            new_row = rows.newRow()
            for f_name, f_value in split_rows[i].items():
                new_row.setValue(f_name, f_value)

            rows.insertRow(new_row)

    return new_fc


def replace_NaT(values) -> list:
    out_value = []
    for value in values:
        if value is pd.NaT:
            out_value.append(None)
        else:
            out_value.append(value)
    return out_value


def get_oid_name(fc_path):
    desc = arcpy.Describe(fc_path)
    oid_name = arcpy.AddFieldDelimiters(fc_path, desc.OIDFieldName).strip('"')
    return oid_name


def get_fields_name(fc_path):
    oid_name = get_oid_name(fc_path)
    fds = [i.name for i in arcpy.ListFields(fc_path)]
    return list(set(fds) - set([oid_name, "Shape_Length", 'Shape_Area']))


class MyRow(object):
    def __init__(self, orignal_row, fields):
        self.orignal_row = orignal_row
        self.fields = fields
        self.geo = orignal_row.getValue('SHAPE')
        self.spa = self.geo.spatialReference
        self.feature_info = Geometry(self.geo).coordinates()

    def _create_polygon(self):
        part_array = arcpy.Array()
        for feature in self.feature_info:
            array = arcpy.Array([arcpy.Point(*coords) for coords in feature])

            array.append(array[0])
            part_array.append(array)

        return arcpy.Polygon(part_array, spatial_reference=self.spa)

    def _create_polyline(self):
        part_array = arcpy.Array()
        for feature in self.feature_info:
            part_array.append(
                arcpy.Array([arcpy.Point(*coords) for coords in feature]))

        return arcpy.Polyline(part_array, spatial_reference=self.spa)

    def _create_point(self):
        pt = arcpy.Point(self.feature_info[0])
        return arcpy.PointGeometry(pt, spatial_reference=self.spa)

    @property
    def copy_row(self):
        if self.geo.type == 'polyline':
            f_info = self._create_polyline()
        elif self.geo.type == 'point':
            f_info = self._create_point()
        elif self.geo.type == 'polygon':
            f_info = self._create_polygon()
        else:
            raise ValueError("Unsupported shape type: {}".format(self.geo.type))
        row = {}
        for fd in self.fields:
            value = self.orignal_row.getValue(f"{fd}")
            row[fd] = value if value is not pd.NaT else None

        row['Shape'] = f_info
        return row


if __name__ == "__main__":
    fc_line = r"D:\5.python项目\data.gdb\dt\test"
    fc_point = r"D:\5.python项目\data.gdb\dt\split_points"
    out_path = r"D:\5.python项目"
    # SplitLineAtPoint(fc_line, fc_point, out_path=None, search_distance=0)
    SplitLineAtPoint(fc_line, fc_point, out_path=out_path, search_distance=0)