数据融合工具（10）线重叠检查修复

craybb

于 2024-07-17 09:08:30 发布

阅读量534

点赞数 21

分类专栏：高效数据融合工具集文章标签：线重叠拓扑数据融合

本文链接：https://blog.csdn.net/craybb/article/details/140484134

版权

高效数据融合工具集专栏收录该内容

10 篇文章 0 订阅

订阅专栏

一、需求背景

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

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节点的位置增加了一个节点。

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

对要素折点进行裂化和聚类以查找共享同一位置（具有通用坐标）的重合要素。
将共有坐标折点插入到共享几何的重合要素中。
运行一系列完整性检查以确定是否违反了为拓扑定义的规则。
在要素数据集中创建潜在拓扑错误的错误日志。

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


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. 设置容差大小