Commit c7cf5b29 authored by xuebingbing's avatar xuebingbing

熊超的积累更新

parent e0c51e5d
Pipeline #501 failed with stages
......@@ -7,12 +7,12 @@ import csv
from connect_solution.coordinate_transformation import *
from util.common import *
from math_tools.data_interpolation import interpolation_dataset_3d
from connect_solution.geographic_tool import *
from util.geographic_tool import *
from math_tools.polynomial_fitting import *
from math_tools.savitsky_golay_filter import *
from math_tools.douglas_peucker_2d import *
from util import const
from smoothing_elevation import SmoothingElevation
class ConnectMeshByMeshList:
......@@ -425,7 +425,7 @@ class ConnectTwoMeshByTopo:
line_blh_points = line_feature.GetGeometryRef().GetPoints()
line_cpt_points = blh_to_cpt(line_blh_points, lon, lat, up, pitch, roll, azimuth)
# 去除很近点,防止后面报错
new_line_cpt_points = self.remove_near_points(interpolation_dataset_3d(line_cpt_points, const.INTERPOLATION_COUNT1))
new_line_cpt_points = self.remove_near_points(interpolation_dataset_2d(line_cpt_points, const.INTERPOLATION_COUNT1))
point_index1, point_index2 = nearest_m_point_of_dataset_2d([0, 0, 0], new_line_cpt_points, 2)
cpt_points = [new_line_cpt_points[point_index1], new_line_cpt_points[point_index2]]
vector = [cpt_points[1][0] - cpt_points[0][0], cpt_points[1][1] - cpt_points[0][1]]
......@@ -436,7 +436,7 @@ class ConnectTwoMeshByTopo:
utm_topo_point = transform_to_utm(topo_blh_point, source_srs)[0][0]
utm_line_points = transform_to_utm(line_blh_points, source_srs)[0]
new_line_utm_points = interpolation_dataset_3d(utm_line_points, const.INTERPOLATION_COUNT1)
new_line_utm_points = interpolation_dataset_2d(utm_line_points, const.INTERPOLATION_COUNT1)
point_index = nearest_m_point_of_dataset_2d(utm_topo_point, new_line_utm_points, 1)[0]
return two_point_distance_2d(utm_topo_point, new_line_utm_points[point_index])
......@@ -529,8 +529,8 @@ class ConnectTwoMeshByTopo:
return True
blh_link_points1 = link_feature1.GetGeometryRef().GetPoints()
blh_link_points2 = link_feature2.GetGeometryRef().GetPoints()
utm_link_points1 = interpolation_dataset_3d(transform_to_utm(blh_link_points1, source_srs)[0], 5)
utm_link_points2 = interpolation_dataset_3d(transform_to_utm(blh_link_points2, source_srs)[0], 5)
utm_link_points1 = interpolation_dataset_2d(transform_to_utm(blh_link_points1, source_srs)[0], 5)
utm_link_points2 = interpolation_dataset_2d(transform_to_utm(blh_link_points2, source_srs)[0], 5)
elevation_bias = self.cal_elevation_bias(utm_topo_point, utm_link_points1, utm_link_points2)
if elevation_bias > const.ELEVATION_THRESHOLD:
output_info = "gps_id为" + str(topo_feature[const.ID]) + \
......@@ -712,8 +712,8 @@ class ConnectTwoMeshByTopo:
def match_line(self, lines1, lines2, line_layer_name, id):
# 对前后得线段集合进行匹配,id为参考线的主键
# lines1和lines2都是cpt坐标系下的点
interpol_lines1 = [interpolation_dataset_3d(line, 2) for line in lines1]
interpol_lines2 = [interpolation_dataset_3d(line, 2) for line in lines2]
interpol_lines1 = [interpolation_dataset_2d(line, 2) for line in lines1]
interpol_lines2 = [interpolation_dataset_2d(line, 2) for line in lines2]
cal_x_mean = lambda line: np.mean([point[0] for point in line[-const.COMPARISON_NUM2:]])
line_info1 = [[cal_x_mean(lines1[i]), i] for i in range(len(interpol_lines1))]
line_info2 = [[cal_x_mean(lines2[j]), j] for j in range(len(interpol_lines2))]
......@@ -814,8 +814,8 @@ class ConnectTwoMeshByTopo:
@staticmethod
def cal_smoothing_count(cpt_dataset1, cpt_dataset2):
# 根据x值得间距计算平滑次数,两数据都指向接边边界
new_cpt_dataset1 = interpolation_dataset_3d(cpt_dataset1, 3)
new_cpt_dataset2 = interpolation_dataset_3d(cpt_dataset2, 3)
new_cpt_dataset1 = interpolation_dataset_2d(cpt_dataset1, 3)
new_cpt_dataset2 = interpolation_dataset_2d(cpt_dataset2, 3)
mean_x1 = np.mean([point[0] for point in new_cpt_dataset1[-const.COMPARISON_NUM2:]])
mean_x2 = np.mean([point[0] for point in new_cpt_dataset2[-const.COMPARISON_NUM2:]])
return int(15 + 15 * np.log(1 + abs(mean_x1 - mean_x2)) / np.log(2))
......@@ -1095,8 +1095,8 @@ class ConnectTwoMeshByTopo:
get_middle_point = lambda point1, point2: [(point1[0] + point2[0]) / 2, (point1[1] + point2[1]) / 2,
(point1[2] + point2[2]) / 2]
connect_boundary.append(get_middle_point(new_cpt_link_points1[-1], new_cpt_link_points2[-1]))
interpolation_cpt_link_points1 = interpolation_dataset_3d(new_cpt_link_points1, 4)
interpolation_cpt_link_points2 = interpolation_dataset_3d(new_cpt_link_points2, 4)
interpolation_cpt_link_points1 = interpolation_dataset_2d(new_cpt_link_points1, 4)
interpolation_cpt_link_points2 = interpolation_dataset_2d(new_cpt_link_points2, 4)
compare_x_list1 = [point[0] for point in interpolation_cpt_link_points1[-const.COMPARISON_NUM1:]]
compare_x_list2 = [point[0] for point in interpolation_cpt_link_points2[-const.COMPARISON_NUM1:]]
......@@ -1122,9 +1122,12 @@ class ConnectTwoMeshByTopo:
:param dataset: utm坐标系下的数据仅计算x, y坐标
:return: 返回的是utm坐标系下的点集
"""
if len(dataset) == 1:
return dataset
dataset1 = dataset[::-1] + []
target_dataset = []
point1 = dataset1.pop()
point2 = dataset1[0]
target_dataset.append(point1)
while dataset1:
point2 = dataset1.pop()
......@@ -1134,6 +1137,10 @@ class ConnectTwoMeshByTopo:
continue
target_dataset.append(point2)
point1 = point2
if not is_same_point_3d(target_dataset[-1], point2):
if two_point_distance_2d(target_dataset[-1], point2) < 0.03:
target_dataset.pop()
target_dataset.append(point2)
return target_dataset
@staticmethod
......@@ -1187,8 +1194,8 @@ class ConnectTwoMeshByTopo:
self.output_info1.append(output_info)
return False, []
utm_dataset1, utm_dataset2 = self.remove_near_points(utm_dataset1), self.remove_near_points(utm_dataset2)
interpol_utm_dataset1 = interpolation_dataset_3d(utm_dataset1, 3) # 按每米3个点进行插值
interpol_utm_dataset2 = interpolation_dataset_3d(utm_dataset2, 3)
interpol_utm_dataset1 = interpolation_dataset_2d(utm_dataset1, 3) # 按每米3个点进行插值
interpol_utm_dataset2 = interpolation_dataset_2d(utm_dataset2, 3)
compare_dataset1 = interpol_utm_dataset1[-const.COMPARISON_NUM2:]
intersection_point1 = intersection_of_two_line(compare_dataset1[0], compare_dataset1[-1],
utm_boundary[0], utm_boundary[-1])
......@@ -1224,7 +1231,7 @@ class ConnectTwoMeshByTopo:
new_utm_dataset2 = []
else:
new_utm_dataset2 = new_utm_dataset2[:len(new_utm_dataset2) - operate_point_num]
compare_dataset = interpolation_dataset_3d(new_compare_dataset1 + new_compare_dataset2[::-1], 2)
compare_dataset = interpolation_dataset_2d(new_compare_dataset1 + new_compare_dataset2[::-1], 2)
flag, info = self.smoothing_dataset(compare_dataset, smoothing_count, utm_boundary, gps_id)
if not flag:
......@@ -1388,12 +1395,11 @@ class ConnectTwoMeshByTopo:
flag, node_feature = self.get_node_by_id(node_layer, node_fields[2], node_id,
mesh_id, ori_node_point, layer_name, source_srs)
if not flag:
return False
return False, None
geom = create_point_geometry(blh_node_point)
node_feature.SetGeometry(geom)
node_layer.SetFeature(node_feature)
node_feature.Destroy()
return True
return True, node_feature
def is_intersect_by_multi_lines(self, node_layer, blh_point, mesh_id, gps_id, layer_name, ids_field, source_srs):
# 判断一个点处是否是多条线的端点,如在路口或者匝道处
......@@ -1447,7 +1453,7 @@ class ConnectTwoMeshByTopo:
return True
return False
def connect_line(self, source_srs, utm_boundary, gps_id, attribute_info, connect_info):
def connect_line(self, source_srs, utm_boundary, gps_id, attribute_info, connect_info, elevation_info):
"""
连接线图层以及对应的点图层,先连接线图层,再连接点图层
:param source_srs: 图层的坐标系
......@@ -1455,6 +1461,7 @@ class ConnectTwoMeshByTopo:
:param gps_id: 接边附近的一个topo点
:param attribute_info: 接边的属性信息
:param connect_info: 接边的连接信息,根据计算获取
:param elevation_info: 平滑高程的时候使用的信息
:return: boolean
"""
if not attribute_info or not connect_info:
......@@ -1512,14 +1519,28 @@ class ConnectTwoMeshByTopo:
line_layer.SetFeature(feature1)
line_layer.SetFeature(feature2)
# 连接线图层对应的节点图层
if not self.connect_point(feature1, blh_points1[-1], direction1, blh_intersection_point,
node_layer, node_fields, layer_name, source_srs) or \
not self.connect_point(feature2, blh_points2[-1], direction2, blh_intersection_point,
node_layer, node_fields, layer_name, source_srs):
flag1, node_feature1 = self.connect_point(feature1, blh_points1[-1], direction1, blh_intersection_point,
node_layer, node_fields, layer_name, source_srs)
if not flag1:
output_info = "gps_id为:" + str(gps_id) + "附近的接边处," + \
str(layer_name) + "的节点没有连接成功,请手动接边"
self.output_info1.append(output_info)
return False
flag2, node_feature2 = self.connect_point(feature2, blh_points2[-1], direction2, blh_intersection_point,
node_layer, node_fields, layer_name, source_srs)
if not flag2:
output_info = "gps_id为:" + str(gps_id) + "附近的接边处," + \
str(layer_name) + "的节点没有连接成功,请手动接边"
self.output_info1.append(output_info)
return False
smoothing_elevation_info = [node_layer, line_layer, [node_feature1, node_feature2], [feature1, feature2],
[direction1, direction2], source_srs, elevation_info[0], elevation_info[1]]
solution = SmoothingElevation()
solution.models(smoothing_elevation_info)
if solution.info:
message = "gps_id为:" + str(gps_id) + "附近的接边处,接边成功,但是高程平滑失败,请手动修改高程"
self.output_info1.append(message)
self.output_info1.extend(solution.info)
return True
@staticmethod
......@@ -1534,13 +1555,6 @@ class ConnectTwoMeshByTopo:
source_srs = self.link_layer.GetSpatialRef() # 获取图层的坐标系
for topo_feature in topo_features:
try:
# if int(topo_feature[const.ID]) == 6602679627452:
# print("dflksjhd")
# else:
# continue
link_ids = [] # 存储已经接过边的道路参考线id,主要应用于多次同路轨迹
mesh_id1, mesh_id2 = self.reset_mesh_id(topo_feature, mesh_id1, mesh_id2)
if self.is_gap_or_big_curve(topo_feature, source_srs):
......@@ -1603,11 +1617,15 @@ class ConnectTwoMeshByTopo:
const.HD_LINK, const.LINK_IDS]
link_connect_info = [[link_feature1], [link_feature2], [[link_direction1, link_direction2]], [[0, 0]],
[link_smoothing_count]]
if self.connect_line(source_srs, utm_connect_boundary, gps_id, link_attribute_info, link_connect_info):
node_fields = [const.LINK_IDS, const.NODE_ID]
link_fields = [const.LINK_ID, const.START_NODE_ID, const.END_NODE_ID]
link_elevation_info = [node_fields, link_fields]
if self.connect_line(source_srs, utm_connect_boundary, gps_id,
link_attribute_info, link_connect_info, link_elevation_info):
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,参考线接边成功"
self.output_info2.append(output_info)
else:
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,参考线接边失败,请手动接边"
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,参考线接边失败或者高程平滑失败,请手动接边"
self.output_info1.append(output_info)
continue
......@@ -1615,12 +1633,15 @@ class ConnectTwoMeshByTopo:
lane_edge_attribute_info = [self.lane_edge_layer, self.lane_edge_node_layer,
[const.START_EDGE_NODE_ID, const.END_EDGE_NODE_ID, const.LANE_EDGE_NODE_ID],
const.HD_LANE_EDGE, const.EDGE_IDS]
lane_edge_node_fields = [const.EDGE_IDS, const.LANE_EDGE_NODE_ID]
lane_edge_fields = [const.LANE_EDGE_ID, const.START_EDGE_NODE_ID, const.END_EDGE_NODE_ID]
lane_edge_elevation_info = [lane_edge_node_fields, lane_edge_fields]
if self.connect_line(source_srs, utm_connect_boundary, gps_id,
lane_edge_attribute_info, lane_edge_connect_info):
lane_edge_attribute_info, lane_edge_connect_info, lane_edge_elevation_info):
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,边界线接边成功"
self.output_info2.append(output_info)
else:
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,边界线接边失败,请手动接边"
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,边界线接边失败或者高程平滑失败,请手动接边"
self.output_info1.append(output_info)
continue
......@@ -1628,12 +1649,15 @@ class ConnectTwoMeshByTopo:
lane_link_attribute_info = [self.lane_link_layer, self.lane_node_layer,
[const.START_LANE_NODE_ID, const.END_LANE_NODE_ID, const.LANE_NODE_ID],
const.HD_LANE_LINK, const.LANE_IDS]
lane_node_fields = [const.LANE_IDS, const.LANE_NODE_ID]
lane_link_fields = [const.LANE_LINK_ID, const.START_LANE_NODE_ID, const.END_LANE_NODE_ID]
lane_link_elevation_info = [lane_node_fields, lane_link_fields]
if self.connect_line(source_srs, utm_connect_boundary, gps_id,
lane_link_attribute_info, lane_link_connect_info):
lane_link_attribute_info, lane_link_connect_info, lane_link_elevation_info):
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,中心线接边成功"
self.output_info2.append(output_info)
else:
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,中心线接边失败,请手动接边"
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,中心线接边失败或者高程平滑失败,请手动接边"
self.output_info1.append(output_info)
continue
......@@ -1641,13 +1665,17 @@ class ConnectTwoMeshByTopo:
link_edge_attribute_info = [self.link_edge_layer, self.link_edge_node_layer,
[const.START_EDGE_NODE_ID, const.END_EDGE_NODE_ID, const.LINK_EDGE_NODE_ID],
const.HD_LINK_EDGE, const.EDGE_IDS]
if self.connect_line(source_srs, utm_connect_boundary,
gps_id, link_edge_attribute_info, link_edge_connect_info):
link_edge_node_fields = [const.EDGE_IDS, const.LINK_EDGE_NODE_ID]
link_edge_fields = [const.LINK_EDGE_ID, const.START_EDGE_NODE_ID, const.END_EDGE_NODE_ID]
link_edge_elevation_info = [link_edge_node_fields, link_edge_fields]
if self.connect_line(source_srs, utm_connect_boundary, gps_id, link_edge_attribute_info,
link_edge_connect_info, link_edge_elevation_info):
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,道路边缘接边成功"
self.output_info2.append(output_info)
else:
if self.link_edge_layer:
output_info = "gps_id为:" + str(gps_id) + "的topo的接边处,道路边缘接边失败,请手动接边"
output_info = "gps_id为:" + str(gps_id) + \
"的topo的接边处,道路边缘接边失败或者高程平滑失败,请手动接边"
self.output_info1.append(output_info)
continue
# 运行到此处,说明所有该接边的图层都接边了
......
......@@ -7,6 +7,7 @@ import sys
import os
from osgeo import osr, ogr, gdal
import csv
from PyQt5.QtCore import *
from util import const
from connect_solution.connect_mesh import ConnectMeshByMeshList
......@@ -77,7 +78,7 @@ def get_adjacent_mesh_box(mesh_id, box_layer, output_info1):
mesh_ids = []
for _ in range(mesh_box_count):
mesh_box_feature = box_layer.GetNextFeature()
if mesh_box_feature is None or mesh_box_feature[const.MESH_ID] == mesh_id:
if mesh_box_feature is None:
continue
mesh_ids.append(mesh_box_feature[const.MESH_ID])
box_layer.SetSpatialFilter(None) # 消除空间过滤
......@@ -85,8 +86,15 @@ def get_adjacent_mesh_box(mesh_id, box_layer, output_info1):
return mesh_ids
def get_current_time():
data = QDateTime.currentDateTime()
currTime = data.toString("yyyy-MM-dd hh:mm:ss")
return currTime
def run(connect_postgresql, sheet_designation_path, mesh_box_path, output_document_path):
# 参数依次为pg库连接字符串、图幅号文件路径,图幅框文件路径、输出文档路径
print(str(get_current_time()) + ": 程序运行开始")
output_info1, output_info2, output_info3 = [], [], []
output_info = "*****" + '-' * 13 + "有关topo图层/图幅框图层没有字段的异常信息,大概率是上游数据的命名不符合规范的问题" + \
'-' * 14 + "*****"
......@@ -99,47 +107,47 @@ def run(connect_postgresql, sheet_designation_path, mesh_box_path, output_docume
output_info = "传入的数据库处理字符串不是PG库,请确认后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
hd_data_source = ogr.Open(connect_postgresql, 1)
if hd_data_source is None:
output_info = "根据连接字符串获取不到数据库,请确认后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
if os.path.exists(sheet_designation_path) is False:
output_info = "根据输入的图幅号列表文件路径,找不到指定的文件,请确认后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
if sheet_designation_path[-3:].lower() != 'txt':
output_info = "传入的图幅号列表文件的格式不是txt格式,运行结束"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
if os.path.exists(mesh_box_path) is False:
output_info = "根据输入的图幅框文件路径没有获取到图幅框文件,请确认后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
if mesh_box_path[-4:] != 'gpkg':
output_info = "传入的图幅框文件格式不是gpkg格式,请确认后重试"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
mesh_box_data_source = ogr.Open(mesh_box_path, 1)
if mesh_box_data_source is None:
output_info = "图幅框文件中没有任何数据,请确认后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
mesh_box_layer = mesh_box_data_source.GetLayerByName(const.MESH_BOX_LAYER_NAME)
if mesh_box_layer is None:
output_info = "从输入的图幅框文件中无法获取图幅框图层,请查验后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
if mesh_box_layer.FindFieldIndex(const.MESH_ID, False) == -1:
output_info = "图幅框图层没有图幅号字段或者图幅号字段的字段名不符合规格文件的命名,程序自动改名"
output_info2.append(output_info)
......@@ -149,7 +157,7 @@ def run(connect_postgresql, sheet_designation_path, mesh_box_path, output_docume
output_info = "图幅框的图幅号字段更改失败,请反馈给研发"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
else:
field_name = ogr.FieldDefn(const.MESH_ID, ogr.OFTString)
if mesh_box_layer.AlterFieldDefn(field_index, field_name, ogr.ALTER_NAME_FLAG) == ogr.OGRERR_NONE:
......@@ -163,28 +171,29 @@ def run(connect_postgresql, sheet_designation_path, mesh_box_path, output_docume
output_info = "输入接边信息文档的路径不是有效的路径,请确认后重试"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
mesh_ids = get_mesh_list(sheet_designation_path, output_info1)
if len(mesh_ids) < 1:
output_info = "传入的图幅号列表文件中没有图幅号,请查验后重试!"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
try:
for mesh_id in mesh_ids:
new_mesh_ids = get_adjacent_mesh_box(mesh_id, mesh_box_layer, output_info1)
if not new_mesh_ids:
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
for new_mesh_id in new_mesh_ids:
if not RefreshTopoNext().run(hd_data_source, new_mesh_id):
output_info = "调用上游程序,未能成功刷新图幅边界上topo点的前驱后继关系"
output_info1.append(output_info)
write_info(output_document_path, output_info1, output_info2, output_info3)
return
return 1
print("调用上游程序修改topo点成功")
except Exception as e:
output_info = "调用上游刷新topo程序时,出现异常,异常为:" + str(e)
output_info1.append(output_info)
......@@ -194,39 +203,49 @@ def run(connect_postgresql, sheet_designation_path, mesh_box_path, output_docume
hd_data_source.Destroy()
if mesh_box_data_source:
mesh_box_data_source.Destroy()
return
return 1
try:
write_info(output_document_path, output_info1, output_info2, output_info3)
ConnectMeshByMeshList(hd_data_source, mesh_box_layer, mesh_ids, output_document_path).main()
print("接边已全部运行完毕,请确认!")
print(str(get_current_time()) + ": 接边已全部运行完毕,请确认!")
if hd_data_source:
hd_data_source.Destroy()
if mesh_box_data_source:
mesh_box_data_source.Destroy()
return 0
except Exception as e:
print("接边程序出发异常,异常为:" + str(e))
finally:
if hd_data_source:
hd_data_source.Destroy()
if mesh_box_data_source:
mesh_box_data_source.Destroy()
return 1
if __name__ == '__main__':
# arg1 = r"PG:dbname='test12' host='localhost' port='5432' user='postgres' password='19931018'"
# # arg1 = r"PG:dbname='qgis_hb_sync' host='win-1.corp.roadlinks.cn' port='5432' user='postgres' password='juefx0123'"
# arg2 = r"C:\Users\熊超\Desktop\edge_match_file\mesh_ids.txt"
# # arg3 = r"C:\Users\熊超\Desktop\edge_match_file\mesh_grid.gpkg"
# arg3 = r"C:\Users\熊超\Desktop\new1\gps_grids.gpkg"
# arg4 = r"C:\Users\熊超\Desktop\edge_match_file"
# run(arg1, arg2, arg3, arg4)
arg1, arg2, arg3, arg4 = None, None, None, None
try:
arg1 = sys.argv[1] # 数据库连接字符串
arg2 = sys.argv[2] # 图幅号的txt文件
arg3 = sys.argv[3] # 图幅框文件
arg4 = sys.argv[4] # 接边报表路径
current_path = os.path.dirname(os.path.abspath(__file__))
gdal_data = current_path + '\\proj'
osr.SetPROJSearchPath(gdal_data)
except Exception as e:
print("传入的参数有问题,程序发生异常,异常为:" + str(e))
#arg1 = r"PG:dbname='qgismap' host='hb.roadlinks.cn' port='5432' user='qgis' password='qgis@123"
# arg2 = r"C:\Users\熊超\Desktop\edge_match_file\mesh_ids.txt"
#arg4 = r"C:\Users\熊超\Desktop\edge_match_file"
#arg3 = r"C:\Users\熊超\Desktop\edge_match_file\mesh_grid.gpkg"
arg1 = r"PG:dbname='qgis_hb_sync' host='win-1.corp.roadlinks.cn' port='5432' user='postgres' password='juefx0123"
arg2 = r"D:\data\meshids.txt"
arg3 = r"D:\data\mesh_grid.gpkg"
arg4 = r"D:\data\edge_match_file"
run(arg1, arg2, arg3, arg4)
# arg1, arg2, arg3, arg4 = None, None, None, None
# try:
# arg1 = sys.argv[1] # 数据库连接字符串
# arg2 = sys.argv[2] # 图幅号的txt文件
# arg3 = sys.argv[3] # 图幅框文件
# arg4 = sys.argv[4] # 接边报表路径
# current_path = os.path.dirname(os.path.abspath(__file__))
# gdal_data = current_path + '\\proj'
# osr.SetPROJSearchPath(gdal_data)
# except Exception as e:
# print("传入的参数有问题,程序发生异常,异常为:" + str(e))
# sys.exit(1)
# sys.exit(run(arg1, arg2, arg3, arg4))
......@@ -34,34 +34,6 @@ def interpolation_dataset_2d(dataset, interpolation_mount):
m = len(dataset)
new_dataset = []
for i in range(m - 1):
new_dataset += interpolation_between_two_point_3d(dataset[i], dataset[i + 1], interpolation_mount)
new_dataset.append(dataset[m - 1])
return new_dataset
def two_point_distance_3d(point1, point2):
# 计算二元坐标的距离
return math.sqrt((point2[0] - point1[0]) ** 2 + (point2[1] - point1[1]) ** 2 + (point2[2] - point1[2]) ** 2)
def interpolation_between_two_point_3d(start_point, end_point, ratio):
# 在起点与终点之间进行插值,ratio是插值的数量
distance = two_point_distance_3d(start_point, end_point)
n = int(distance * ratio) # 插值的个数
dataset = []
result = np.linspace(0, 1, n + 2)[:n + 1] # 保留首,不保留尾
for k in result:
x = k * (end_point[0] - start_point[0]) + start_point[0]
y = k * (end_point[1] - start_point[1]) + start_point[1]
z = k * (end_point[2] - start_point[2]) + start_point[2]
dataset.append([x, y, z])
return dataset
def interpolation_dataset_3d(dataset, interpolation_mount):
m = len(dataset)
new_dataset = []
for i in range(m - 1):
new_dataset += interpolation_between_two_point_3d(dataset[i], dataset[i + 1], interpolation_mount)
new_dataset += interpolation_between_two_point_2d(dataset[i], dataset[i + 1], interpolation_mount)
new_dataset.append(dataset[m - 1])
return new_dataset
......@@ -13,14 +13,9 @@ def two_point_distance_2d(point1, point2):
return math.sqrt((point2[0] - point1[0]) ** 2 + (point2[1] - point1[1]) ** 2)
def two_point_distance_3d(point1, point2):
# 计算三元坐标点的距离
return math.sqrt((point2[0] - point1[0]) ** 2 + (point2[1] - point1[1]) ** 2 + (point2[2] - point1[2]) ** 2)
def interpolation_between_two_point(start_point, end_point, ratio):
def interpolation_between_two_point_2d(start_point, end_point, ratio):
# 在起点与终点之间进行插值,ratio是插值的数量
distance = two_point_distance_3d(start_point, end_point)
distance = two_point_distance_2d(start_point, end_point)
n = int(distance * ratio) # 插值的个数
dataset = []
result = np.linspace(0, 1, n + 2)[:n + 1] # 保留首,不保留尾
......@@ -36,10 +31,19 @@ def interpolation_between_two_point(start_point, end_point, ratio):
return dataset
def interpolation_dataset_2d(dataset, interpolation_mount):
# 按照二维坐标插值
m = len(dataset)
new_dataset = []
for i in range(m - 1):
new_dataset += interpolation_between_two_point_2d(dataset[i], dataset[i + 1], interpolation_mount)
new_dataset.append(dataset[m - 1])
return new_dataset
class CalElevation:
"""
使用局部加权计算新的插值点的高程值
通过高斯函数将距离映射为权重,带宽取距离之和除以一个常数
使用局部加权计算高程值,并且通过高斯函数将距离映射为权重
"""
def __init__(self, m, n):
......@@ -68,7 +72,7 @@ class CalElevation:
def get_elevation(self, ori_points, new_points):
# new_points是需要计算高程的点集,ori_points是原始点集
new_dataset = interpolation_dataset(ori_points, self.n)
new_dataset = interpolation_dataset_2d(ori_points, self.n)
ori_xy_list = [[point[0], point[1]] for point in new_dataset]
tree = KDTree(np.array(ori_xy_list), leaf_size=20) # 二维kd数
for i in range(1, len(new_points) - 1):
......@@ -86,15 +90,6 @@ class CalElevation:
return new_points
def interpolation_dataset(dataset, interpolation_mount):
m = len(dataset)
new_dataset = []
for i in range(m - 1):
new_dataset += interpolation_between_two_point(dataset[i], dataset[i + 1], interpolation_mount)
new_dataset.append(dataset[m - 1])
return new_dataset
class ImprovedDouglasPeucker(object):
"""
1. 抽稀后得结果不要两个点两个点得形式
......@@ -141,6 +136,7 @@ class ImprovedDouglasPeucker(object):
self.returned_list.append(points[-1])
self.returned_list.append(points[0])
else:
# 在这里分段会产生一个重点,目的是不让最后产生的点集是两个点两个点的效果,最后需要去重才能得到正确的效果
sequence1 = points[:index + 1]
sequence2 = points[index:]
sequences = [sequence1, sequence2]
......@@ -180,5 +176,5 @@ class ImprovedDouglasPeucker(object):
self.dilute(self.not_dilute_list.pop())
self.returned_list = self.returned_list[::-1] # 在抽稀的过程中反向了
tmp_points = self.remove_duplicated_point(self.returned_list)
new_points = interpolation_dataset(tmp_points, self.interpolation_count)
new_points = interpolation_dataset_2d(tmp_points, self.interpolation_count)
return CalElevation(self.m, self.n).get_elevation(points, new_points)
......@@ -257,7 +257,7 @@ class SavitzkyGolayFilter:
3. 对滤波后的数据进行抽稀和插值
"""
def __init__(self, interpolation_count=0.125):
def __init__(self):
self.trans = PosNegTrans() # 坐标系转移类
self.poly_fit = PolyFit() # 多项式拟合类
......
# _*_ coding: utf-8 _*_ #
# @Time: 2020/6/28 16:35
# @Author: XiongChao
# _*_ coding: utf-8 _*_ #
# @Time: 2020/6/28 10:03
# @Author: XiongChao
# 修改指定图层的高程以达到平滑的要求
from osgeo import gdal
from util.common import *
from util.geographic_tool import *
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") # 支持中文路径
gdal.SetConfigOption("SHAPE_ENCODING", "") # 使属性表支持中文
ogr.RegisterAll() # 注册所有驱动
DEF_Z_RANGE_THRESHOLD = 4 # 判断是否进行高程操作的阈值
DEF_SLOPE_THRESHOLD = np.tan(30 / 180 * np.pi) # detla_z / dis的阈值
DEF_ALPHA = 0.0000001
def two_point_distance_by_blh(blh_point1, blh_point2, source_srs):
# 根据经纬度计算两点的距离,使用的是utm坐标系
utm_point1 = transform_to_utm([blh_point1], source_srs)[0][0]
utm_point2 = transform_to_utm([blh_point2], source_srs)[0][0]
return two_point_distance_2d(utm_point1, utm_point2)
class SmoothingElevation:
"""
对输入的点线的高程进行平滑
"""
def __init__(self):
self.info = []
def is_reset_elevation(self, line_features, line_directions, line_layer_name, primary_key):
# 判断是否调整高程,以及返回节点的高程
# 如果几条线的高程值的极差超过2米,就认为异常,不修复高程
z_list = []
for i in range(len(line_features)):
blh_points = line_features[i].GetGeometryRef().GetPoints()
if line_directions[i] == 1:
z_list.extend([point[2] for point in blh_points[-2:]])
else:
z_list.extend([point[2] for point in blh_points[:2]])
z_range = max(z_list) - min(z_list)
if z_range > DEF_Z_RANGE_THRESHOLD:
info = str(line_layer_name) + "上id: " + str(line_features[0][primary_key]) + \
"为的节点,附近关联的线的z值的极差过大,不进行处理,请修改数据"
self.info.append(info)
return False, 0
return True, np.mean(z_list)
@staticmethod
def reset_node_elevation(node_layer, node_feature, elevation):
# 重置节点的高程
node_point = node_feature.GetGeometryRef().GetPoints()[0]
new_node_point = [node_point[0], node_point[1], elevation]
node_geom = create_point_geometry(new_node_point)
node_feature.SetGeometry(node_geom)
node_layer.SetFeature(node_feature)
@staticmethod
def cal_slope(blh_point1, blh_point2, source_srs):
distance = two_point_distance_by_blh(blh_point1, blh_point2, source_srs)
detla_z = abs(blh_point1[2] - blh_point2[2])
return detla_z / distance
def reset_line_elevation(self, line_layer, line_feature, line_direction,
elevation, source_srs, line_layer_name, primary_key):
blh_points = [list(point) for point in line_feature.GetGeometryRef().GetPoints()]
if line_direction == -1:
blh_points = blh_points[::-1]
blh_points[-1] = [blh_points[-1][0], blh_points[-1][1], elevation]
j = len(blh_points) - 1
if j <= 2:
slope = self.cal_slope(blh_points[len(blh_points) - 1], blh_points[len(blh_points) - 2], source_srs)
# 如果满足条件,直接将点集写入即可
if slope > DEF_SLOPE_THRESHOLD:
# 此处没办法平滑
info = "id: " + str(line_feature[primary_key]) + "的线,其长度很短并且无法满足指定角度的要求"
self.info.append(info)
return False
else:
flag = False
for j in range(len(blh_points) - 2, int((len(blh_points) + 1) / 2) - 1, -1):
slope = self.cal_slope(blh_points[len(blh_points) - 1], blh_points[j], source_srs)
if slope < DEF_SLOPE_THRESHOLD + 0.001:
# 0.001是为了防止浮点数计算带来的误差
flag = True
break
if not flag:
# 表示没有找到可以达到平滑要求的点
info = str(line_layer_name) + "上id: " + str(line_feature[primary_key]) + \
"的线,其要么长度太短,要么高程值异常"
self.info.append(info)
return False
start_point = blh_points[len(blh_points) - 1]
start_elevation, end_elevation = start_point[2], blh_points[j][2]
distance = two_point_distance_by_blh(blh_points[len(blh_points) - 1], blh_points[j], source_srs)
cal_elevation = lambda point: two_point_distance_by_blh(start_point, point, source_srs) / \
distance * (end_elevation - start_elevation) + start_elevation
for k in range(len(blh_points) - 2, j, -1):
blh_points[k][2] = cal_elevation(blh_points[k])
if line_direction == -1:
blh_points = blh_points[::-1]
geom = create_line_geometry(blh_points)
line_feature.SetGeometry(geom)
line_layer.SetFeature(line_feature)
return True
def reset_elevation(self, node_layer, line_layer, line_layer_name, node_features,
line_features, line_directions, line_fields, source_srs):
# node_fields为列表,表示用到的节点图层的字段,0-ids字段,1-id字段
# line_fields为列表,表示用到的线图层的字段,0-id字段,1-起点字段,2-终点字段
flag, elevation = self.is_reset_elevation(line_features, line_directions, line_layer_name, line_fields[0])
if not flag:
return
if not self.reset_line_elevation(line_layer, line_features[0], line_directions[0],
elevation, source_srs, line_layer_name, line_fields[0]) or \
not self.reset_line_elevation(line_layer, line_features[1], line_directions[1],
elevation, source_srs, line_layer_name, line_fields[0]):
return
self.reset_node_elevation(node_layer, node_features[0], elevation)
self.reset_node_elevation(node_layer, node_features[1], elevation)
@staticmethod
def remove_near_points(line_layer, line_features, source_srs, alpha=0.01):
# 删除线上过近的点,远近的描述为alpha
for line_feature in line_features:
blh_points = line_feature.GetGeometryRef().GetPoints()
utm_points, utm_srs = transform_to_utm(blh_points, source_srs)
new_utm_points = remove_near_data(utm_points, alpha)
new_blh_points = transform_to_src(new_utm_points, source_srs, utm_srs)
geom = create_line_geometry(new_blh_points)
line_feature.SetGeometry(geom)
line_layer.SetFeature(line_feature)
def models(self, input_info):
node_layer, line_layer, node_features, line_features, \
line_directions, source_srs, node_fields, line_fields = input_info
line_layer_name = line_layer.GetName()
try:
self.remove_near_points(line_layer, line_features, source_srs, 0.01)
except Exception as e:
message = str(line_layer_name) + "在高程平滑之前的去除过近点的过程中出现异常:" + str(e)
self.info.append(message)
return
self.reset_elevation(node_layer, line_layer, line_layer_name, node_features,
line_features, line_directions, line_fields, source_srs)
return
......@@ -60,9 +60,9 @@ def three_point_dot_product(point1, point2, point3):
return result
def is_same_point_2d(point1, point2):
def is_same_point_2d(point1, point2, alpha=DEF_ALPHA):
# 判断两个点是不是一个点,仅选择x, y
flag = np.array([np.abs(point1[i] - point2[i]) < 0.000001 for i in range(len(point1))])
flag = np.array([np.abs(point1[i] - point2[i]) < alpha for i in range(len(point1))])
flag = flag[:2] # 当为三元数值得时候,仅取前两位进行比较
if flag.all():
return True
......@@ -105,3 +105,23 @@ def nearest_m_point_of_dataset_2d(point, dataset, m):
dist, ind = tree.query([np.array(point[:2])], k=m)
index_list = list(ind[0])
return index_list
def remove_near_data(dataset, alpha=DEF_ALPHA):
# 将数据集中非常接近的点去掉,每一个feature至少有起点和终点,并且这两个点不同
dataset1 = dataset[::-1] + []
target_dataset = []
point1 = dataset1.pop()
point2 = dataset1[0]
target_dataset.append(point1)
while dataset1:
point2 = dataset1.pop()
if is_same_point_2d(point1, point2, alpha):
continue
target_dataset.append(point2)
point1 = point2
if not is_same_point_3d(target_dataset[-1], point2):
if two_point_distance_2d(target_dataset[-1], point2) < alpha:
target_dataset.pop()
target_dataset.append(point2)
return target_dataset
# _*_ coding: utf-8 _*_ #
# @Time: 2020/6/1 11:40
# @Author: XiongChao
from osgeo import ogr, osr
import math
def transform_to_utm(blh_dataset, source_srs):
lon = blh_dataset[0][0]
utm_zone = 31 + math.floor(lon / 6.0)
utm_srs = osr.SpatialReference()
utm_srs.SetUTM(utm_zone, 1)
src_to_utm_trans = osr.CreateCoordinateTransformation(source_srs, utm_srs)
utm_dataset = []
for blh_point in blh_dataset:
temp_point = blh_point
utm_point = list(src_to_utm_trans.TransformPoint(temp_point[0], temp_point[1]))
utm_point[2] = blh_point[2]
utm_dataset.append(utm_point)
return utm_dataset, utm_srs
def transform_to_src(utm_dataset, source_srs, utm_srs):
utm_to_src_trans = osr.CreateCoordinateTransformation(utm_srs, source_srs)
src_dataset = []
for utm_point in utm_dataset:
temp_point = utm_point
src_point = list(utm_to_src_trans.TransformPoint(temp_point[0], temp_point[1]))
src_point[2] = utm_point[2]
src_dataset.append(src_point)
return src_dataset
def create_line_geometry(dataset):
geo = ogr.Geometry(ogr.wkbLineStringZM)
for i in range(len(dataset)):
geo.SetPoint(i, dataset[i][0], dataset[i][1], dataset[i][2])
return geo
def create_point_geometry(point):
geo = ogr.Geometry(ogr.wkbPointZM)
geo.SetPoint(0, point[0], point[1], point[2])
return geo
# _*_ coding: utf-8 _*_ #
# @Time: 2020/6/1 11:40
# @Author: XiongChao
from osgeo import ogr, osr
import math
def transform_to_utm(blh_dataset, source_srs):
lon = blh_dataset[0][0]
utm_zone = 31 + math.floor(lon / 6.0)
utm_srs = osr.SpatialReference()
utm_srs.SetUTM(utm_zone, 1)
src_to_utm_trans = osr.CreateCoordinateTransformation(source_srs, utm_srs)
utm_dataset = []
for blh_point in blh_dataset:
temp_point = blh_point
utm_point = list(src_to_utm_trans.TransformPoint(temp_point[0], temp_point[1]))
utm_point[2] = blh_point[2]
utm_dataset.append(utm_point)
return utm_dataset, utm_srs
def transform_to_src(utm_dataset, source_srs, utm_srs):
utm_to_src_trans = osr.CreateCoordinateTransformation(utm_srs, source_srs)
src_dataset = []
for utm_point in utm_dataset:
temp_point = utm_point
src_point = list(utm_to_src_trans.TransformPoint(temp_point[0], temp_point[1]))
src_point[2] = utm_point[2]
src_dataset.append(src_point)
return src_dataset
def create_line_geometry(dataset):
geo = ogr.Geometry(ogr.wkbLineStringZM)
for i in range(len(dataset)):
geo.SetPoint(i, dataset[i][0], dataset[i][1], dataset[i][2])
return geo
def create_point_geometry(point):
geo = ogr.Geometry(ogr.wkbPointZM)
geo.SetPoint(0, point[0], point[1], point[2])
return geo
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment