一、需求背景

        先明确一下“线重叠”的定义。

        ArcGIS拓扑工具集中的拓扑规则：

• 不能自重叠（线） —线要素不得与自身重叠。这些线要素可以交叉或接触但不得有重合的线段。此规则适用于街道等线段可能接触闭合线的要素，但同一街道不得出现两次相同的路线。

• 不能重叠(线) —线不得与同一要素类（或子类型）中的线重叠。例如，当河流要素类中线段不能重复时，使用此规则。线可以交叉或相交，但不能共享线段。

 1.1 拓扑错误修复

        对于“不能自重叠”错误，使用“简化”（从要素中移除导致错误的自重叠或自相交的线段），此修复可以应用到一个或多个错误。

        对于“不能重叠”错误，使用“移除重叠”（移除导致错误的要素重叠部分），此修复只能一条记录一条记录修复。

        需求来了，能不能，一个工具，解决所有的“重叠”错误，并记录修复情况呢？

        毕竟，重叠错误可能是这样的？多如烟海……

1.2 剖析一下

        01 说一说？自重叠

        我们知道，线是由一系列点按给定的顺序依次连接构成的。线自重叠，只需要线节点连接出现往复或压盖即可发生，如下图所示。

 

 

        在理想数据下，实现线自重叠的判断，自重叠部分的查找以及自重叠的修复代码如下：

from shapely.geometry import LineString, Point, MultiLineString
from shapely.ops import linemerge
import numpy as np
from itertools import combinations


# 将折线转换为线段
def create_segments(polyline):
    segments = []

    # 如果输入是列表或ndarray
    if isinstance(polyline, (list, np.ndarray)):
        coords = polyline if isinstance(polyline, list) else polyline.tolist()
        for i in range(len(coords) - 1):
            line = LineString([coords[i], coords[i + 1]])
            segments.append(line)

    # 如果输入是LineString
    elif isinstance(polyline, LineString):
        coords = list(polyline.coords)
        for i in range(len(coords) - 1):
            line = LineString([coords[i], coords[i + 1]])
            segments.append(line)

    # 如果输入是MultiLineString
    elif isinstance(polyline, MultiLineString):
        for line_string in polyline.geoms:
            coords = list(line_string.coords)
            for i in range(len(coords) - 1):
                line = LineString([coords[i], coords[i + 1]])
                segments.append(line)

    else:
        raise TypeError("当前数据类型不支持！")

    return segments


def calculate_slope(line):
    x1, y1 = line.coords[0]
    x2, y2 = line.coords[-1]
    if x2 - x1 != 0:
        slope = (y2 - y1) / (x2 - x1)
    else:
        slope = np.inf  # 垂直线的斜率定义为无穷大
    return slope


def judge_overlapping(multi_line):
    slopes = {}
    for idx, line in enumerate(multi_line):
        slope = calculate_slope(line)
        if slope not in slopes:
            slopes[slope] = []
        slopes[slope].append(idx)

    for slope, line_indices in slopes.items():
        if len(line_indices) <= 1:
            continue
        for i, j in combinations(line_indices, 2):
            processed_points = set()  # 用于记录已处理过的点
            line1 = multi_line[i]
            line2 = multi_line[j]
            if line1.distance(line2) != 0:
                continue
            coincident_points_count = 0
            for p1 in line1.coords:
                if p1 in processed_points:  # 如果点已处理过，跳过计数
                    continue
                if line2.distance(Point(p1)) == 0:
                    coincident_points_count += 1
                    processed_points.add(p1)  # 添加已处理的点
                    if coincident_points_count >= 2:
                        return True  # 发现重叠，直接返回True
            else:
                for p2 in line2.coords:
                    if p2 in processed_points:  # 如果点已处理过，跳过计数
                        continue
                    if line1.distance(Point(p2)) == 0:
                        coincident_points_count += 1
                        processed_points.add(p2)  # 添加已处理的点
                        if coincident_points_count >= 2:
                            return True  # 发现重叠，直接返回True
                else:
                    continue
    return False  # 没有发现重叠，返回False


def find_overlapping_sections(multi_line):
    slopes = {}
    for idx, line in enumerate(multi_line):
        slope = calculate_slope(line)
        slope = np.round(slope, decimals=8)
        if slope not in slopes:
            slopes[slope] = []
        slopes[slope].append(idx)

    overlapping_sections = []
    for slope, line_indices in slopes.items():
        if len(line_indices) <= 1:
            continue
        for i, j in combinations(line_indices, 2):
            processed_points = set()  # 用于记录已处理过的点
            line1 = multi_line[i]
            line2 = multi_line[j]
            if not np.isclose(line1.distance(line2), 0):
                continue
            coincident_points_count = 0
            for p1 in line1.coords:
                if p1 in processed_points:  # 如果点已处理过，跳过计数
                    continue
                if np.isclose(line2.distance(Point(p1)), 0):
                    coincident_points_count += 1
                    processed_points.add(p1)  # 添加已处理的点
            for p2 in line2.coords:
                if p2 in processed_points:  # 如果点已处理过，跳过计数
                    continue
                if np.isclose(line1.distance(Point(p2)), 0):
                    coincident_points_count += 1
                    processed_points.add(p2)  # 添加已处理的点
            if coincident_points_count >= 2:
                overlapping_sections.append((i, j))

    return overlapping_sections


def export_overlapping(multi_line, overlapping_sections):

    # 输出重叠部分
    intersecting_sections = []
    for i, j in overlapping_sections:
        over_part = multi_line[i].intersection(multi_line[j], grid_size=0)
        intersecting_sections.append(over_part)
    union_lines = linemerge(intersecting_sections)
    if isinstance(union_lines, MultiLineString):
        # 如果结果是一个 MultiLineString，则将其拆解为单独的几何体
        overlapping = list(union_lines.geoms)
    elif isinstance(union_lines, LineString):
        # 如果结果是一个单独的 LineString，则将其包装成列表
        overlapping = [union_lines]
    else:
        # 如果结果是其他类型，这里可以根据具体情况进行处理
        overlapping = []

    return overlapping


def repair_overlapping(segments, overlapping_sections):
    multi_line = segments[::]

    # 修复多段线
    for i, j in overlapping_sections:
        if not multi_line[i] or not multi_line[j]:
            continue
        if multi_line[i].equals(multi_line[j]):
            multi_line[i] = []

        elif multi_line[i].contains(multi_line[j]):
            multi_line[j] = []

        elif multi_line[j].contains(multi_line[i]):
            multi_line[i] = []

        else:
            multi_line[i] = multi_line[i].difference(multi_line[j])

    repaired_line = linemerge([i for i in multi_line if i])

    if isinstance(repaired_line, MultiLineString):
        # 如果结果是一个 MultiLineString，则将其拆解为单独的几何体
        overlapping = list(repaired_line.geoms)
    elif isinstance(repaired_line, LineString):
        # 如果结果是一个单独的 LineString，则将其包装成列表
        overlapping = [repaired_line]
    else:
        # 如果结果是其他类型，这里可以根据具体情况进行处理
        overlapping = []

    return overlapping


if __name__ == "__main__":

    # 示例：检查多段线是否存在重叠部分
    test_line = LineString([(0, 0), (1, 1), (2, 2), (3, 3), (0.5, 0.5)])
    example_multiline = create_segments(test_line)
    result = find_overlapping_sections(example_multiline)
    print("重叠部分的线段编号：", result)

    export = export_overlapping(example_multiline, result)
    print("重叠部分：", export)

    repaired = repair_overlapping(example_multiline, result)
    print("修复后线段：", repaired)

        线自重叠一般很少发生，但是通过对自重叠的研究分析后，发现一些很有意思的事。如下图：

        按我的估计，80%以上的人会认为162、164和163三个节点在一条直线上，即线162-163的斜率和线162-164的斜率一定相等。为什么使用斜率来检测是否共线呢？

        斜率，这是一个初略过滤两条线段是否重叠的基础，是一种非常高效简洁的检测思路。

        事实上，从坐标值来看，这3个点还真不一定共线（实际上就碰到了）。但在ArcGIS中分析和查看，他们却一定是共线的。这里涉及到一个知识点，XY容差。

---

X,Y 容差指的是坐标之间的最小距离，小于该距离的坐标将合并到一起。

默认 X,Y 容差设置为 0.001 米，或以数据集的坐标系单位表示的等效值。例如，如果坐标系以英尺为单位，则此默认值是 0.003281 英尺（0.03937 英寸）。默认值是默认 X,Y 分辨率的 10 倍，且在大多数情况下均推荐此设置。如果坐标以经纬度表示，则默认 X,Y 容差为 0.0000000556 度。

---

        而点位坐标的存储精度，是通过X,Y分辨率来表达的。X,Y 分辨率（表示非常小的距离）是指用于存储 x,y 坐标值的有效数字的位数。

        接着说……

        接着之前的疑问，当162、164和163三个点坐标值不共线，但在ArcGIS拓扑分析时，为什么这条折线会出现自重叠呢？

        我们知道，在创建完成拓扑数据集后，需要添加拓扑规则，最后验证拓扑，就能按拓扑规则的内容找出错误的数据记录。

        也就是这个过程，发现了一些奥妙……

        在对“重叠”验证的过程，实际上对线要素的节点进行了分析和规整，它在出现重叠的点位（或附近）进行了插入节点处理，并对线节点的顺序进行了重新调整，以确保线的is_vaild属性和is_simple属性能正确识别线。

        162-163-164节点构成的线，在使用拓扑验证前，该条线有3个节点，is_simple属性是TRUE，在使用拓扑验证后，is_simple属性为FALSE，要素的节点变成4个，其中在164节点的位置增加了一个节点。

---

验证拓扑包含以下四个过程：

1. 对要素折点进行裂化和聚类以查找共享同一位置（具有通用坐标）的重合要素。

2. 将共有坐标折点插入到共享几何的重合要素中。

3. 运行一系列完整性检查以确定是否违反了为拓扑定义的规则。

4. 在要素数据集中创建潜在拓扑错误的错误日志。

---

        为了实现该“验证”处理的内容，方便在以后写工具时解决线要素的重整，小编附上实现代码：

from shapely.geometry import LineString, Point
from rtree import index
import numpy as np


def create_point_rtree(points):
    p = index.Property()
    p.dimension = 2
    idx = index.Index(properties=p)

    for i, pt in enumerate(points):
        point = Point(pt)
        idx.insert(i, point.bounds, point)

    return idx


def verification_line(line):
    line_coords = np.array(line.coords)

    idx = create_point_rtree(line_coords)

    new_coords = []

    for i in range(len(line_coords) - 1):
        segment_coords = [line_coords[i], line_coords[i + 1]]
        segment = LineString(segment_coords)

        intersecting_nodes = idx.intersection(segment.bounds, objects=True)
        ins_seg_nodes = [i.object for i in intersecting_nodes if segment.intersects(i.object.buffer(1e-8))]

        ls_nodes = [node for node in ins_seg_nodes if not np.isclose(np.array(segment.coords) - np.array(node.coords), 0).any()]

        new_coords.append(segment_coords[0])
        if len(ls_nodes) == 0:
            continue

        elif len(ls_nodes) > 1:

            # 排序后的节点插入到新坐标串
            sorted_nodes = sorted(ls_nodes,
                                  key=lambda node: np.linalg.norm(np.array(node.coords) - line_coords[0]))

            # 逐一添加排序后的节点坐标到 new_coords
            for node in sorted_nodes:
                new_coords.append(np.array(node.coords))
        else:
            new_coords.append(np.array(ls_nodes[0].coords))

    new_coords.append(line_coords[-1])

    new_line_coords = np.vstack(new_coords)

    return LineString(new_line_coords)


if __name__ == "__main__":
    # Example usage:
    multi_line = LineString([(0, 0), (3, 3), (5, 5), (0.4999999999, 0.5)])
    updated_multi_line = verification_line(multi_line)
    print(updated_multi_line)

        拓扑容差，在ArcGIS中普遍使用。 拓扑容差是一条术语，通常用于表示两种容差：X,Y 容差和 Z 容差。拓扑容差的默认值是坐标分辨率的 10 倍。

---

实用提示

以下是一些有关拓扑容差的实用提示：

1. 设置容差大小

• 建议使用 10 倍于 X 和 Y 分辨率的容差值，这通常会得到良好的结果。

• 典型的 X 和 Y 容差要比数据采集的实际精度小。

• 需要平衡 X、Y 容差，太小可能导致无法正确处理重叠边界的线作业，而太大可能导致要素坐标重叠，影响地图表达的精度。

2. X、Y 容差和数据采集精度

• X、Y 容差不应接近数据采集的精度，这可能导致坐标移动过大或过小。

• 在地图比例为 1:12,000 时，即使 1/50 英寸也对应 20 英尺，这种精度难以在数字化和扫描转换中实现。

3. 坐标等级设置

• 在拓扑中，可以为每个要素类设置坐标精度等级，使更精确的要素具有更高的等级。

• 更高等级的要素可以调整较低等级但更不精确的要素，以确保坐标的一致性。

4. 控制移动要素类

• 在聚类过程中，可以设置拓扑规则以控制可能移动的要素类。

• 通过为要素类分配不同的等级，在拓扑容差范围内，低等级要素的折点会被捕捉到高等级要素的折点，以及同级要素折点位置的几何平均。

---

        02 线重叠

        线重叠示意图如下：

        线重叠，就比较好理解，就是两条或多条线要素，在空间位置上存在重叠和覆盖。如上图所示，折线1、9、2和折线10发生重叠，折线3与折线4发生重叠。若只是为了找出这种重叠，可使用“成对相交”工具或对该要素类创建拓扑数据集，使用拓扑规则“不能重叠”解决。

        但要实现“重叠”部分的修复，就不是一件简单的事。

        比如下图中，修复要素8和要素9重叠部分时，“移除重叠”功能需要你选择要保留的要素。

        从几何的角度分析，选择保留的要素，需要考虑线的走向和连接；从属性的角度分析，选择保留的要素，某些属性值应该与线的上下连接要素保持一致。

        数据实际处理过程中，往往是优先考虑几何的走向，属性值采取属性传递或单独填写。一是属性值一般不靠谱，都需要分析判断后选择，二是属性值靠谱，无论选谁都可以。基于这两点情况，我们在做批量“线重叠”修复时，先考虑几何上的合理性。

        具体算法如下图所示描述：

二、线重叠检查修复工具

        线重叠检查修复工具，对线要素类图层中，线要素或线要素间存在重叠的部分进行查找和自动修复。

2.1 工具概述

        “线重叠检查修复修复”工具，支持的功能如下：

1. 对线自重叠部分进行查找输出，支持自动修复功能；

2. 对线重叠部分进行查找输出，支持批量、合理性自动修复。

2.2 功能流程

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

        （2）工具参数介绍如下：

        （3）工具输出：

        当要素类中存在自重叠错误时：

        输出“{要素类名}_self_overlay”，记录自重叠修复的位置。

        当要素类中存在线重叠错误时：

        输出{要素类名}_overlay”，记录线重叠修复的位置。

        若“是否另存输出”未选中，将在scratch.gdb中备份输入要素，并对其进行重叠修复。默认将直接对输入要素进行重叠修复。

        （4）注意事项：

        当工具参数“是否另存输出”未选中时，将对输入的要素进行原地修改，运行工具前，应注意数据备份。