一、曲线线段修复
在GIS中，曲线一般有两种表示方式：一类是由多个直线段（多段线）组成的曲线，另一类是由贝塞尔曲线、圆弧和椭圆弧等数学曲线表示的曲线。
-
多段线（LineString）： 这种曲线是由一系列相邻的直线段连接而成。在GIS中，多段线通常用于近似曲线或折线特征，如道路、河流等。多段线是由一系列节点（顶点）连接而成，可以通过这些节点来近似实际的曲线。
-
贝塞尔曲线、圆弧和椭圆弧： 这些曲线是由数学公式生成的光滑曲线，其路径由控制点、半径等参数决定。在GIS中，这些曲线通常用于绘制更为光滑和精确的曲线，如图形设计、CAD制图、地图标注等。它们可以更准确地模拟自然曲线，提供更高的图形质量。
在GIS中，多段线和贝塞尔曲线、圆弧、椭圆弧等曲线都有各自的应用场景和优势。
1.1 多段线（LineString）的应用
-
折线地物表示： 多段线通常用于近似表示折线状的地物，如河流、道路等。由一系列直线段组成，适用于模拟具有直线特征的地物。
-
轨迹记录： 在移动物体轨迹的记录中，多段线可以用于表示物体在不同时间点之间的运动路径。
-
地图标注： 多段线可以用于绘制地图上的标注线，例如连接地图元素之间的注释线。
1.2 贝塞尔曲线、圆弧和椭圆弧的应用
-
曲线地物表示： 贝塞尔曲线、圆弧和椭圆弧等数学曲线可以用于更准确地表示自然界中的光滑曲线地物，如环岛、山脊的轮廓等。
-
图形设计： 在GIS中，用于图形设计的曲线类型可以提供更高的图形质量，用于绘制符号、标签等，使地图更具艺术性。
-
CAD绘图： 在CAD（计算机辅助设计）软件中，贝塞尔曲线、圆弧和椭圆弧常用于精确绘制建筑、工程图纸等。
1.3 将曲线增密为多段线的原因
-
数据存储和传输： 曲线的数学表示可能需要更多的数据存储和传输资源，因为曲线的定义通常需要更多的参数以准确描述其形状。这可能包括曲线的方程、曲线上的控制点等。与之相比，直线或多段线的数学表示可能更简单，只需存储端点坐标即可。
-
分析和处理效率： 曲线的分析和处理可能相对复杂，而多段线的分析通常更为高效。在某些GIS分析任务中，将曲线增密为多段线可能更容易实现和处理。
-
符号化和显示： 在某些情况下，曲线的图形表示可能无法很好地适应显示设备或符号化需求。增密为多段线可以更好地适应不同的显示要求。
二、如何将曲线线段（贝塞尔、圆弧和椭圆弧）替换为线段？
2.1 使用ArcGIS解决
使用ArcGIS提供的地理处理工具“增密”，通过指定“增密方法”为“偏移”，实现对曲线线段替换为线段。工具直接修改源要素图层，没有其他输出。
增密工具参数示意如下：
实际上，因为没有输出，这种“做好事不留名”的做法，在一些场景下是不友好的。比如，我不希望直接在源数据上进行更改。
从数据质检需求和需要保持源数据不变的原则下，当我们需要对曲线线段（贝塞尔、圆弧和椭圆弧）进行查找和修复时，需要单独输出曲线线段的记录和最后修复后的结果。
2.2 使用Arcpy自定义曲线线段替换为线段的实现
“曲线线段修复”工具将曲线线段（贝塞尔、圆弧和椭圆弧）替换为线段。
工具概述
（1）将输入要素中存在曲线线段（贝塞尔、圆弧和椭圆弧）的要素输出；
（2）修复输入要素中的曲线线段（贝塞尔、圆弧和椭圆弧），并输出到指定目录。
（3）对线要素（多部件）的曲线部分进行增密处理，不破坏原有要素类的其他属性，如字段属性，属性域关系等。
（4）默认“最大偏移偏差（输出线段与原始线段之间的最大距离）”为0.1米；
功能流程
工具执行界面如下图所示：
工具运行结果展示如下图：
工具参数介绍如下：
工具输出：
在“输出路径”下，输出与输入同名的要素类，并将被转换为折线的曲线部分输出，图层名为“curve_repair”。
所有输出的要素类都存放在输出目录中的scratch.gdb中。若scratch.gdb不存在，则自动创建，若已存在，不会覆盖其中已有的要素。
注意事项：
无。
实现代码如下：
import math
import os
from collections import defaultdict
from itertools import groupby
from operator import itemgetter
import arcpy
from arcgis.geometry import Geometry
from pyxmy.feature.database import Row, create_feature
from pyxmy.feature.field import TextField, FieldTable
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 * math.pi) / 180
g_length = math.sqrt((x * math.cos(degree_angle)) ** 2 + (y * math.sin(degree_angle)) ** 2)
return g_length
# 指定要素，创建空要素类模板
def create_fc_template(fc_path, out_path):
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
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 parse_curve(curve_coords):
rs = defaultdict(list)
curvePaths = curve_coords['curvePaths']
for n, part in enumerate(curvePaths):
for i, p in enumerate(part):
if isinstance(p, dict):
rs[n].append(i)
for key in rs:
values = rs[key]
groups = []
for k, g in groupby(enumerate(values), lambda i_x: i_x[0] - i_x[1]):
groups.append(list(map(itemgetter(1), g)))
rs[key] = groups
return rs
# 计算两个坐标值之间的距离
def distance(p1, p2):
x1, y1 = p1
x2, y2 = p2
return (x2 - x1) ** 2 + (y2 - y1) ** 2
# 找出p在ls中，与ls中最接近的点位
def min_index(p, ls):
distances = [distance(p, point) for point in ls]
min_index = distances.index(min(distances))
return min_index
def repair_curve(fc_polyline, out_path):
# 遍历每条线，并查询空间索引以获取在指定距离范围内的点
new_fc = create_fc_template(fc_polyline, out_path)
rows = arcpy.InsertCursor(new_fc)
desc = arcpy.Describe(fc_polyline)
OID = desc.OIDFieldName
spatial_ref = desc.spatialReference
cursor = arcpy.UpdateCursor(fc_polyline)
row = cursor.next()
spa_type = spatial_ref.type
if spa_type == "Geographic":
g_radius = plane_to_geodesic(1)
else:
g_radius = 1
records = []
arcpy.env.scratchWorkspace = out_path
gdb_path = arcpy.env.scratchGDB
while row:
geo_row = row.getValue('SHAPE')
line_id = row.getValue(OID)
if geo_row.hasCurves: # 几何具有曲线段
curve_coords = Geometry(geo_row)
rs = parse_curve(curve_coords)
# 可以沿线要素添加折点，将曲线线段（贝塞尔、圆弧和椭圆弧）替换为线段
densify_row = geo_row.densify('DISTANCE', distance=g_radius, deviation=g_radius/10)
densify_curve_coords = Geometry(densify_row)
# 只对曲线线段进行增密处理，即替换为线段
for part, values in rs.items(): # part为部件索引，values为部件内曲线线段点连续索引集合
part_coords = densify_curve_coords['paths'][part]
curve_part_coords = curve_coords['curvePaths'][part]
current_index = 0
for value in values:
start_i = value[0] + current_index
end_i = value[-1] + current_index
try:
if start_i - 1 >= 0:
start_p = curve_part_coords[start_i - 1]
else:
start_p = curve_part_coords[start_i]
except IndexError:
start_p = curve_part_coords[start_i]
try:
if end_i + 1 < len(Geometry(geo_row).get_part(part)):
end_p = curve_part_coords[end_i + 1]
else:
end_p = curve_part_coords[end_i]
except IndexError:
end_p = curve_part_coords[end_i]
# 将曲线线段（贝塞尔、圆弧和椭圆弧）替换为线段
rp_start_i = curve_part_coords.index(start_p)
rp_end_i = curve_part_coords.index(end_p)
if isinstance(start_p, dict):
first_key = next(iter(start_p))
start_p = start_p[first_key][0]
if isinstance(end_p, dict):
first_key = next(iter(end_p))
end_p = end_p[first_key][0]
new_start_i = min_index(start_p, part_coords)
new_end_i = min_index(end_p, part_coords)
# 曲线增密部分的坐标串
replace_curve_coords = part_coords[new_start_i:new_end_i + 1]
curve_coords['curvePaths'][part][rp_start_i:rp_end_i+1] = replace_curve_coords
# 将曲线替换为圆弧的部分输出
f_info = arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in replace_curve_coords]), spatial_reference=spatial_ref)
f_row = Row([line_id], 'POLYLINE', geometry=f_info, attr_name=['IN_FID'])
records.append(f_row)
current_index += (new_end_i - new_start_i) - (rp_end_i - rp_start_i)
# 构建新的线段
for feature in curve_coords['curvePaths']:
new_part = arcpy.Polyline(
arcpy.Array([arcpy.Point(*coords) for coords in feature]), spatial_reference=spatial_ref)
row.setValue("SHAPE", new_part)
rows.insertRow(row)
else:
rows.insertRow(row)
row = cursor.next()
del row
del rows
tf = TextField(name='IN_FID', alias='输入ID', length=10, null_able=False)
out_fc = os.path.join(gdb_path, "{}".format("curve_repair"))
cursor = create_feature(feature_name=out_fc, field=FieldTable(tf), geometry_type='POLYLINE', srs=spatial_ref)
cursor.insert(records)
if __name__ == "__main__":
fc_line = r"E:\6.测试代码\Default.gdb\aa"
out_path = r'E:\6.测试代码'
repair_curve(fc_line, out_path)
2403
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
