测绘程序设计试题集——基础篇（编写内容仅作参考）
参考教材——李英冰老师的《测绘程序设计试题集》（2017版）
一、 原文
仅展示试题1的原文，如需了解其余内容，请自行购买图书
李老师著作主页-http://ybli.users.sgg.whu.edu.cn/book/
二、数据集
测站名,时间（YYYYMMDD）,经度（°）,纬度（°）,大地高（m）,高度角（°）
P01,20170728,116.407143,39.913607,45.234,20
Р02,20170728,125.460913,43.850217,102.432,30
РОЗ,20170728,82.650403,63.604725,251.672,40
Р04,20170728,107.37637,26.043428,35.277,45
POS,20170728,112.674792,1.514577,24.556,50
Р06,20170728,138.283829,-30.160507,176.889,55
P07,20170728,102.961019,-77.21885,435.279,60
三、代码
# 读取数据
import math
with open('坐标信息.txt','r',encoding='utf-8') as file:
data = file.read()
data = data.split('\n')
new_data = []
for i in data:
ii = i.split(',')
new_data.append(ii)
new_data1 = new_data[1:]
# 湿分量投影函数计算
aw15 = 0.00058021897
bw15 = 0.0014275268
cw15 = 0.043472961
aw30 = 0.0056794847
bw30 = 0.0015138625
cw30 = 0.046729510
aw45 = 0.0058118019
bw45 = 0.0014572752
cw45 = 0.043908931
aw60 = 0.0059727542
bw60 = 0.0015007428
cw60 = 0.044626982
aw75 = 0.0061641693
bw75 = 0.0017599082
cw75 = 0.54736038
MWE_list = []
for i in range(len(new_data1)):
weidu = float(new_data1[i][3])
E = float(new_data1[i][5])
E1 = math.radians(E)
if weidu < 15:
a = aw15
b = bw15
c = cw15
elif 15 <= weidu < 30:
a = aw15 + (weidu - 15) * (aw30 - aw15) / (30 - 15)
b = bw15 + (weidu - 15) * (bw30 - bw15) / (30 - 15)
c = cw15 + (weidu - 15) * (cw30 - cw15) / (30 - 15)
elif 30 <= weidu < 45:
a = aw30 + (weidu - 30) * (aw45 - aw30) / (45 - 30)
b = bw30 + (weidu - 30) * (bw45 - bw30) / (45 - 30)
c = cw30 + (weidu - 30) * (cw45 - cw30) / (45 - 30)
elif 45 <= weidu < 60:
a = aw45 + (weidu - 45) * (aw60 - aw45) / (60 - 45)
b = bw45 + (weidu - 45) * (bw60 - bw45) / (60 - 45)
c = cw45 + (weidu - 45) * (cw60 - cw45) / (60 - 45)
elif 60 <= weidu < 75:
a = aw60 + (weidu - 60) * (aw75 - aw60) / (75 - 60)
b = bw60 + (weidu - 60) * (bw75 - bw60) / (75 - 60)
c = cw60 + (weidu - 60) * (cw75 - cw60) / (75 - 60)
else:
a = aw75
b = bw75
c = cw75
MWE = (1 / (1 + a / (1 + b / (1 + c)))) / (1 / (math.sin(E1) + a / (math.sin(E1) + b / (math.sin(E1) + c))))
MWE_list.append(MWE)
# 干分量投影函数计算
ahg15 = 0.0012769934
bhg15 = 0.0029153695
chg15 = 0.062610505
ahg30 = 0.0012683230
bhg30 = 0.0029152299
chg30 = 0.062837393
ahg45 = 0.0012465397
bhg45 = 0.0029288445
chg45 = 0.063721774
ahg60 = 0.0012196049
bhg60 = 0.0029022565
chg60 = 0.063824265
ahg75 = 0.0012045996
bhg75 = 0.0029024912
chg75 = 0.064258455
ahp15 = 0.0
bhp15 = 0.0
chp15 = 0.0
ahp30 = 0.000012709626
bhp30 = 0.000021414979
chp30 = 0.000090128400
ahp45 = 0.000026523662
bhp45 = 0.000030160779
chp45 = 0.000043497037
ahp60 = 0.000034000452
bhp60 = 0.000072562722
chp60 = 0.0084795348
ahp75 = 0.000041202191
bhp75 = 0.00011723375
chp75 = 0.0017037206
MDE_list = []
for i in range(len(new_data1)):
weidu = float(new_data1[i][3])
E = float(new_data1[i][5])
E1 = math.radians(E)
aht = 2.53e-5
bht = 5.49e-3
cht = 1.14e-3
H = float(new_data1[i][4])
t0 = 28
tt = [31,28,31,30,31,30,31,31,30,31,30,31]
year = int(new_data1[i][1][0:4])
month = int(new_data1[i][1][4:6])
day = int(new_data1[i][1][6:8])
# 闰年判断
if (year % 4 ==0 and year % 100 != 0) or (year % 400 == 0):
tt[1] = 29
nianjiri = sum(tt[0:month-1]) + day
else:
nianjiri = sum(tt[0:month-1]) + day
if weidu < 15:
a = ahg15 + ahg15 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg15 + bhg15 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg15 + chg15 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
elif 15 <= weidu < 30:
a = ahg15 + (weidu - 15) * (ahg30 - ahg15) / (30 - 15) \
+ ahg15 + (weidu - 15) * (ahg30 - ahg15) / (30 - 15) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg15 + (weidu - 15) * (bhg30 - bhg15) / (30 - 15) \
+ bhg15 + (weidu - 15) * (bhg30 - bhg15) / (30 - 15) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg15 + (weidu - 15) * (chg30 - chg15) / (30 - 15) \
+ chg15 + (weidu - 15) * (chg30 - chg15) / (30 - 15) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
elif 30 <= weidu < 45:
a = ahg30 + (weidu - 30) * (ahg45 - ahg30) / (45 - 30) \
+ ahg30 + (weidu - 30) * (ahg45 - ahg30) / (45 - 30) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg30 + (weidu - 30) * (bhg45 - bhg30) / (45 - 30) \
+ bhg30 + (weidu - 30) * (bhg45 - bhg30) / (45 - 30) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg30 + (weidu - 30) * (chg45 - chg30) / (45 - 30) \
+ chg30 + (weidu - 30) * (chg45 - chg30) / (45 - 30) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
elif 45 <= weidu < 60:
a = ahg45 + (weidu - 45) * (ahg30 - ahg45) / (60 - 45) \
+ ahg45 + (weidu - 45) * (ahg30 - ahg45) / (60 - 45) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg45 + (weidu - 45) * (bhg30 - bhg45) / (60 - 45) \
+ bhg45 + (weidu - 45) * (bhg30 - bhg45) / (60 - 45) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg45 + (weidu - 45) * (chg30 - chg45) / (60 - 45) \
+ chg45 + (weidu - 45) * (chg30 - chg45) / (60 - 45) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
elif 60 <= weidu < 75:
a = ahg60 + (weidu - 60) * (ahg75 - ahg60) / (75 - 60) \
+ ahg60 + (weidu - 60) * (ahg75 - ahg60) / (75 - 60) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg60 + (weidu - 60) * (bhg75 - bhg60) / (75 - 60) \
+ bhg60 + (weidu - 60) * (bhg75 - bhg60) / (75 - 60) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg60 + (weidu - 60) * (chg75 - chg60) / (75 - 60) \
+ chg60 + (weidu - 60) * (chg75 - chg60) / (75 - 60) * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
else:
a = ahg75 + ahg75 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
b = bhg75 + bhg75 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
c = chg75 + chg75 * math.cos(2 * math.pi * (nianjiri - t0) / 365.25)
ss = (1 / (1 + a / (1 + b / (1 + c)))) / (1 / (math.sin(E1) + a / (math.sin(E1) + b / (math.sin(E1) + c))))
MDE = ss + (1 / math.sin(E1) - (1 / (1 + aht / (1 + bht / (1 + cht)))) / (1 / (math.sin(E1) + aht /
(math.sin(E1) + bht / (math.sin(E1) + cht))))) * (H / 1000)
MDE_list.append(MDE)
# 延迟改正计算
ZHD_list = []
for i in range(len(new_data1)):
H = float(new_data1[i][4])
ZHD = 2.29951 * math.e ** (-0.000116 * H)
ZHD_list.append(round(ZHD,3))
ZWD = 0.1
SS_list = []
for i in range(len(new_data1)):
SS = ZHD_list[i] * MDE_list[i] + ZWD * MDE_list[i]
SS_list.append(round(SS, 3))
# 输出
print('测站名 高度角 ZHD MDE ZWD MWE SS')
for i in range(len(new_data1)):
print(new_data1[i][0], new_data1[i][-1], ZHD_list[i], MDE_list[i], ZWD, MWE_list[i],SS_list[i])
# 保存文件
with open('8.txt', 'w', encoding='utf-8') as file1:
file1.write('测站名 高度角 ZHD MDE ZWD MWE SS')
file1.write('\n')
for i in range(len(new_data1)):
file1.write(str(new_data1[i][0]) + ' ' + str(new_data1[i][-1]) + ' ' +
str(ZHD_list[i]) + ' ' + str(MDE_list[i]) + ' ' +
str(ZWD) + ' ' + str(MWE_list[i])+' ' + str(SS_list[i]))
file1.write('\n')
四、代码问题
- 因作者水平有限，且时间冲突等原因，以上部分代码中未使用方便阅读的变量名，还请读者包含，自行理解；
- 第二问以后计算结果存在问题；
- 作者研究方向不涉及电离层和对流层，因此对其算法不理解，以上代码仅可参考。