大地高转正常高

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
# !/usr/bin/python3
# -*- coding: utf-8 -*-
"""
@Author
: yuong
@Version
: python3
@Organization
: SWJTU
@Email
: lych@my.swjtu.edu.cn
------------------------------------
@File
:
straight_fitting.py
@Description
:
@CreateTime
:
2021/12/7 17:00
------------------------------------
@ModifyTime
:
"""
import math
import re
import os
rootP = os.getcwd()
def getXYHh():
'''获取起始点的位置数据'''
knownID = open(".\SL_known.txt", encoding='utf-8').read().splitlines()
a = knownID[1].split(',')
b = knownID[2].split(',')
a = [float(a[0]), float(a[1]), float(a[2]), float(a[3])]
b = [float(b[0]), float(b[1]), float(b[2]), float(b[3])]
return [a, b]
def getInert():
'''获取插入点的位置数据'''
# rootPath=os.getcwd()
insertID = open(".\h_file.XYH", encoding='utf-8')
allID0 = insertID.read().splitlines()
# The remaining line of remarks
allID1 = []
for IDtxt in allID0[1:]:
IDtxt = re.split(",", IDtxt)
allID1.append(IDtxt)
return allID1
def fitting():
'''算法实现直线拟合的关键参数：高程异常变化率'''
t = getXYHh()
startPiont = t[0]
endPiont = t[1]
Sab = math.sqrt((endPiont[0] - startPiont[0]) * (endPiont[0] - startPiont[0]) + (
endPiont[1] - endPiont[1]) * (endPiont[1] - endPiont[1])) / 1000
# unit to km
gs = startPiont[2] - startPiont[3]
# Starting point elevation outlier
g = ((endPiont[2] - endPiont[3]) - gs) * 100
# unit to cm
gi = float(g / Sab)
if gi < 5:
print("满足直线拟合条件！")
return [gs, gi]
else:
print("不满足每km高程异常变化梯度小于5cm的条件，建议更换其它拟合模型！")
def distanceP():
'''判断插入点到拟合平面直线的垂向距离是否满足1km的要求；求取满足要求点位到起始点
距离在S的投影长度；返回点位的H、在S的投影长'''
insertID = getInert()
t = getXYHh()
a = t[0]
b = t[1]
agreeID = []
x = []
for x in insertID:
# L = abs((y1-y2)x3+(x2-x1)y3+x1y2-y1x2)/sqrt((y1-y2) ^ 2 +(x1-x2) ^ 2)
L = abs((a[1] - b[1]) * float(x[1]) + (b[0] - a[1]) * float(x[2]) + a[0] * b[1] - a[1]
* b[0]) / (math.sqrt((a[1] - b[1]) * (a[1] - b[1]) + (a[1] - b[1]) * (a[1] - b[1])))
hypo = math.sqrt((a[0] - float(x[1])) * (a[0] - float(x[1])) +
(a[1] - float(x[2])) * (a[1] - float(x[2])))
j = math.asin(L / hypo)
zs = []
if L < 1000:
z1 = hypo * math.cos(j)
zs.append(x[0])
zs.append(float(x[1]))
zs.append(float(x[2]))
zs.append(float(x[3]))
zs.append(z1)
agreeID.append(zs)
else:
print(x[0] + "距离超过1km,该点不适直线拟合法，建议更换其它拟合模型")
continue
# print(agreeID)
return agreeID
def answerh():
insertID = distanceP()
key = fitting()
result = ''
print(insertID)
print(key)
for hi in insertID:
hj = key[0] + key[1] * hi[4]
h = hi[3] - hj
result = str(hi[0]) + ',' + str(hi[1]) + ',' + 
str(hi[2]) + ',' + str(hi[3]) + ',' + str(h)
with open(rootP + "\h_fittingResult_SL.txt", "a") as f:
f.write(result + 'n')
f.close()
return
answerh()

1