python实现全站仪三角高程测量计算

一、介绍

在大比例尺地形图测绘过程中,通过三角高程测量来获取控制点的高程。

已知起始点A的高程,首先测站定为A,目标站定为B,盘左盘右分别读取竖直角读数(建议设置为垂直角90)和斜距,量取仪器高目标高,由此A到B往测结束。之后在B点架设全站仪,返测A点,重复上述观测方法,完成A到B往返观测。以此类推。

 

 二、组织数据

 

刚刚读取的盘左盘右竖直角即为半测回值,取平均求得一测回值。

现已知往测返测的竖直角读数(一测回值),目标高,仪器高,斜距,创建一个记事本:

输入格式:

  1. 高程起始点
  2. J1(往)竖直角的度,分,秒,J1与J2测站间斜距,仪器高,目标高,竖直角正负值(+/-),J1(返)竖直角度,分,秒,J1与J2斜间距,仪器高

 三、python计算三角高程

import math
def dms_to_radians(degree):
    #获取度分秒
    degrees_str=degree[0]
    minutes_str=degree[1]
    seconds_str=degree[2]
    #转化为数字
    degrees=float(degrees_str)
    minutes=float(minutes_str)
    seconds=float(seconds_str)
    #计算弧度
    radians=math.radians(degrees+(minutes/60)+(seconds/3600))
    return radians

#数据准备
v_degree_list=[] #竖直角弧度
D_list=[] #测站斜距
h_instrument_list=[]
h_target_list=[]
v_degree_list_back=[]#返回观测竖直角弧度
D_list_back=[] #测站返回斜距
h_instrument_list_back=[] #返回测站仪器高

filename='/home/z/hight_02'
with open(filename,'r') as f:
    fist_line=f.readline().strip() #起始点高程
    for line in f:
        items_list=[] #存储竖直角
        data=line.strip().split(',')
        item1,item2,item3,item4,item5,item6,item7=data[0],data[1],data[2],data[3],data[4],data[5],data[6]
        D_list.append(float(item4))
        h_instrument_list.append(float(item5))
        h_target_list.append(float(item6))
        items_list.append(item1)
        items_list.append(item2)
        items_list.append(item3)
        if item7=='+':
            v_degree_list.append(dms_to_radians(items_list))
        elif item7=='-':
            v_degree_list.append(dms_to_radians(items_list)*(-1))
        #算一下返回的高差
        a_list=[] #存储返回竖直角
        a1,a2,a3,a4,a5=data[7],data[8],data[9],data[10],data[11]
        D_list_back.append(float(a4))
        h_instrument_list_back.append(float(a5))
        a_list.append(a1)
        a_list.append(a2)
        a_list.append(a3)
        v_degree_list_back.append(dms_to_radians(a_list))

print("\n竖直角弧度:",v_degree_list)
print("竖直角弧度(返回):",v_degree_list_back)
print("\n测站斜距:",D_list)
print("测站斜距(返回):",D_list_back)
print("\n仪器高:",h_instrument_list)
print("仪器高(返回):",h_instrument_list_back)
print("\n目标高:",h_target_list)

#计算高差
H_detla_list=[]
for a,b,c,d in zip(v_degree_list,D_list,h_instrument_list,h_target_list):
    h=math.sin(a)*b+c-d
    H_detla_list.append(h)
print("\n高差:",H_detla_list)
H_detla_list_back=[] #返回测回高差
for a,b,c,d in zip(v_degree_list_back,D_list_back,h_instrument_list_back,h_target_list):
    h=math.sin(a)*b+c-d
    H_detla_list_back.append(h)
print("高差(返回):",H_detla_list_back)

H_detla_list_average_abs=[]
H_detla_list_average=[]
for m,n in zip(H_detla_list,H_detla_list_back):
    h=(abs(m)+abs(n))/2
    if m>0:
        H_detla_list_average.append(h)
    elif m<0:
        H_detla_list_average.append(h*(-1))
print("\n往返高差均值:",H_detla_list_average)


#计算高差改正数
H_sum=sum(H_detla_list_average)
print("\n理论之差:",H_sum)

D_hor_list=[]
for a,b in zip(v_degree_list,D_list): #计算平距
    D_hor_list.append(math.cos(a)*b)
print("\n平距:",D_hor_list)
D_sum=sum(D_hor_list)  #在此只计算往测平距,未取平均,可修改
detla=[] #高差改正数
for item in D_hor_list:
    detla.append((item/D_sum)*H_sum)
sign=-1
detla=list(map(lambda x:x*sign,detla))
print("\n高差改正数:",detla)

H_detla_new_list=[] #改正后高差
for m,n in zip(H_detla_list_average,detla):
    H_detla_new_list.append(m+n)
print("\n改正后高差:",H_detla_new_list)

#高程计算
Hight_list=[]
H=float(fist_line)
for item in H_detla_new_list:
    Hight_list.append(H)
    H+=item
print("\n高程:",Hight_list)

#高程检核
fist_line_check=Hight_list[-1]+H_detla_new_list[-1]
print("\n起始点高程检核值:",fist_line_check)

运行结果:

 

 

  • 3
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
高分一号影像一般是指GF-1卫星拍摄的影像。要计算高程需要用到DEM数据,这里以读取GeoTiff格式的DEM文件为例,使用GDAL库实现。 以下是一个基于GDAL库的Python函数,实现对GF-1影像进行高程计算: ```python import gdal def calculate_elevation(gf1_image_path, dem_path): """ 计算GF-1影像的高程,返回一个和GF-1影像分辨率一致的高程数组 :param gf1_image_path: GF-1影像路径 :param dem_path: DEM数据路径 :return: 高程数组 """ # 读取GF-1影像数据 gf1_image = gdal.Open(gf1_image_path) width = gf1_image.RasterXSize height = gf1_image.RasterYSize transform = gf1_image.GetGeoTransform() projection = gf1_image.GetProjection() gf1_band = gf1_image.GetRasterBand(1) gf1_data = gf1_band.ReadAsArray(0, 0, width, height) # 读取DEM数据 dem = gdal.Open(dem_path) dem_band = dem.GetRasterBand(1) dem_data = dem_band.ReadAsArray(0, 0, width, height) # 计算高程 gf1_elevation = gf1_data - dem_data # 输出高程数据 driver = gdal.GetDriverByName("GTiff") output_path = gf1_image_path[:-4] + "_elevation.tif" output_image = driver.Create(output_path, width, height, 1, gf1_band.DataType) output_image.SetGeoTransform(transform) output_image.SetProjection(projection) output_image.GetRasterBand(1).WriteArray(gf1_elevation) output_image.FlushCache() return gf1_elevation ``` 这个函数接受两个参数,`gf1_image_path`是GF-1影像的路径,`dem_path`是DEM数据的路径。函数首先使用GDAL库读取GF-1影像和DEM数据,并获取影像和DEM数据的宽度、高度、仿射变换参数和投影信息。然后,函数将GF-1影像和DEM数据作差,得到高程数据。最后,函数将高程数据保存为一个GeoTiff格式的文件,并返回高程数据。 需要注意的是,函数中的计算结果是一个和GF-1影像分辨率一致的高程数组,每个元素代表一个像素的高程值。如果需要得到高程值的物理单位,需要根据DEM数据的信息进行换算。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值