概述
# !/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()
在这里插入代码片
最后
以上就是靓丽春天为你收集整理的大地高转正常高的全部内容,希望文章能够帮你解决大地高转正常高所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
发表评论 取消回复