python简单计算卫星位置

前言

一、思路

总的可分为两个部分：获取导航参数和计算卫星位置。
①获取导航参数：首先讲导航星历中的数据切片，存入csv文件中，再读取csv文件的数据并赋值给各参数
②计算卫星位置：首先要进行时间数据的换算，然后依次计算各步骤即可

二、自动读取导航星历并计算卫星位置


import csv
import math  as mwith open('D:\\学习课件\\GNSS\\brdc0320.22n\\brdc0320.22n', 'r') as f:if f == 0:print("不能打开文件！")else:print("导航文件打开成功！")nfile_lines = f.readlines()  # 按行读取N文件print(len(nfile_lines))f.close()def start_num():  #定义数据记录的起始行for i in range(len(nfile_lines)):if nfile_lines[i].find('END OF HEADER')!=-1:start_num=i+1return start_num
n_dic_list=[]n_data_lines_nums=int((len(nfile_lines)-start_num())/8)
print("一共%d组数据"%(n_data_lines_nums))#第j组，第i行
for j in range(n_data_lines_nums):n_dic = {}for i in range(8):data_content = nfile_lines[start_num() + 8 * j + i]n_dic['数据组数'] = j + 1if i == 0: n_dic['卫星PRN号'] = int(data_content.strip('\n')[0:2].strip(' '))n_dic['历元'] = data_content.strip('\n')[3:22]n_dic['卫星钟偏差(s)'] = float((data_content.strip('\n')[23:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))  # 利用字符串切片功能来进行字符串的修改n_dic['卫星钟漂移(s/s)'] = float((data_content.strip('\n')[42:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['卫星钟漂移速度(s/s*s)'] = float((data_content.strip('\n')[61:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))if i == 1:n_dic['IODE'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['C_rs'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['n'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['M0'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[60:79][-3:]).strip(' '))if i == 2:n_dic['C_uc'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['e'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['C_us'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['sqrt_A'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))if i == 3:n_dic['TEO'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['C_ic'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['OMEGA'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[ 42:60][-3:]).strip(' '))n_dic['C_is'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))if i == 4:n_dic['I_0'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['C_rc'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['w'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['OMEGA_DOT'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))if i == 5:n_dic['IDOT'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['L2_code'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['PS_week_num'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['L2_P_code'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))if i == 6:n_dic['卫星精度(m)'] = float((data_content.strip('\n')[4:22][0:-4] + 'e' + data_content.strip('\n')[4:22][-3:]).strip(' '))n_dic['卫星健康状态'] = float((data_content.strip('\n')[22:41][0:-4] + 'e' + data_content.strip('\n')[23:41][-3:]).strip(' '))n_dic['TGD'] = float((data_content.strip('\n')[41:60][0:-4] + 'e' + data_content.strip('\n')[42:60][-3:]).strip(' '))n_dic['IODC'] = float((data_content.strip('\n')[60:79][0:-4] + 'e' + data_content.strip('\n')[61:79][-3:]).strip(' '))n_dic_list.append(n_dic)with open( 'D:\\学习课件\\GNSS\\brdc0320.csv', 'w', newline='') as f:header=['数据组数', '卫星PRN号', '历元', '卫星钟偏差(s)', '卫星钟漂移(s/s)', '卫星钟漂移速度(s/s*s)','IODE', 'C_rs', 'n', 'M0', 'C_uc', 'e', 'C_us', 'sqrt_A', 'TEO', 'C_ic', 'OMEGA', 'C_is','I_0', 'C_rc', 'w', 'OMEGA_DOT', 'IDOT', 'L2_code', 'PS_week_num', 'L2_P_code','卫星精度(m)', '卫星健康状态', 'TGD', 'IODC']writer = csv.DictWriter(f, fieldnames=header)writer.writeheader()writer.writerows(n_dic_list)
f.close()with open('D:\\学习课件\\GNSS\\brdc0320.csv','rt') as csvfile:reader = csv.DictReader(csvfile)for row in reader:if row['数据组数']=='1':PRN=float(row["卫星PRN号"])TIME=row["历元"]year=int(TIME.strip('\n')[0:3])month=int(TIME.strip('\n')[3:6])day=int(TIME.strip('\n')[6:9])hour=int(TIME.strip('\n')[9:12])minute=int(TIME.strip('\n')[12:15])second=float(TIME.strip('\n')[15:20])a_0=float(row["卫星钟偏差(s)"]) a_1=float(row["卫星钟漂移(s/s)"])a_2=float(row["卫星钟漂移速度(s/s*s)"])IODE=float(row["IODE"])C_rs=float(row["C_rs"])δn=float(row["n"])M0=float(row["M0"])C_uc=float(row["C_uc"])e=float(row["e"])C_us=float(row["C_us"])sqrt_A=float(row["sqrt_A"])TEO=float(row["TEO"])C_ic=float(row["C_ic"])OMEGA=float(row["OMEGA"])C_is=float(row["C_is"])I_0=float(row["I_0"])C_rc=float(row["C_rc"])w=float(row["w"])OMEGA_DOT=float(row["OMEGA_DOT"])IDOT=float(row["IDOT"])L2_code=float(row["L2_code"])PS_week_num=float(row["PS_week_num"])L2_P_code=float(row["L2_P_code"])#卫星精度(m)=2#卫星健康状态=0TGD=float(row["TGD"])IODC=float(row["IODC"])print(IODC)  print(sqrt_A)          #卫星位置的计算
def CulLocation(PRN,year,month,day,hour,minute,second,a_0,a_1,a_2,IDOT,C_rs,δn,M0,C_uc,e,C_us,sqrt_A,TEO,C_ic,OMEGA,C_is,I_0,C_rc,w,OMEGA_DOT,t1):#1.计算卫星运行平均角速度 GM:WGS84下的引力常数 =3.986005e14，a:长半径GM = 398600500000000n_0 = m.sqrt(GM)/m.pow(sqrt_A,3)n = n_0 + δn   #2.计算归化时间t_k 计算t时刻的卫星位置  UT：世界时 此处以小时为单位  UT = hour + (minute / 60.0) + (second / 3600); #GPS时起始时刻1980年1月6日0点   year是两位数 需要转换到四位数  if year>=80:  if year==80 and month==1 and day<6 : year = year + 2000  else:year = year + 1900   else:year = year + 2000 if month<=2 :  year = year - 1month = month + 12 # 1，2月视为前一年13，14月#需要将当前需计算的时刻先转换到儒略日再转换到GPS时间  JD = (365.25 * year) + int(30.6001 * (month + 1)) + day + UT / 24 + 1720981.5;  WN = int((JD - 2444244.5) / 7);  #WN:GPS_week number 目标时刻的GPS周 t_oc = (JD - 2444244.5 - 7.0 * WN) * 24 * 3600.0; #t_GPS:目标时刻的GPS秒  #对观测时刻t1进行钟差改正,注意：t1应是由接收机接收到的时间if t1 is None:t_k=0else:δt = a_0+a_1(t1-t_oc)+a_2(t1-t_oc)^2 t = t1-δtt_k = t-TEO# if t_k > 302400:#     t_k -= 604800 # else :#     t_k += 604800 #3.平近点角计算M_k = M_0+n*t_kprint(M0)M_k = M0 + n * t_k #实际应该是乘t_k，但是没有接收机的观测时间，所以为了练手设t_k=0 #4.偏近点角计算 E_k  (迭代计算) E_k = M_k + e*sin(E_k)E= 0;  E1 = 1;  count = 0;  while abs(E1-E)>1e-10 :  count = count +1 E1 = E;  E = M_k +e*m.sin(E);  if count>1e8:  print("计算偏近点角时未收敛！") break  #5.计算卫星的真近点角 V_k = m.atan((m.sqrt(1 - e*e) * m.sin(E))/( m.cos(E) - e));  #6.计算升交距角 u_0(φ_k), ω:卫星电文给出的近地点角距u_0=V_k+w #7.摄动改正项 δu、δr、δi :升交距角u、卫星矢径r和轨道倾角i的摄动量δu=C_uc*m.cos(2*u_0)+C_us*m.sin(2*u_0)δr=C_rc*m.cos(2*u_0)+C_rs*m.sin(2*u_0)δi=C_ic*m.cos(2*u_0)+C_is*m.sin(2*u_0)#8.计算经过摄动改正的升交距角u_k、卫星矢径r_k和轨道倾角 i_ku=u_0+δur=m.pow(sqrt_A,2)*(1-e*m.cos(E))+δri=I_0+δi+IDOT * (t_k); #实际乘t_k=t-t_oe #9.计算卫星在轨道平面坐标系的坐标,卫星在轨道平面直角坐标系（X轴指向升交点）中的坐标为：x_k=r*m.cos(u)y_k=r*m.sin(u)#10.观测时刻升交点经度Ω_k的计算，升交点经度Ω_k等于观测时刻升交点赤经Ω与格林尼治恒星时GAST之差  Ω_k=Ω_0+(ω_DOT-omega_e)*t_k-omega_e*t_oeomega_e = 7.292115e-5 #地球自转角速度  OMEGA_k= OMEGA + (OMEGA_DOT - omega_e) *t_k - omega_e * TEO;  #星历中给出的Omega即为Omega_o=Omega_t_oe-GAST_w#11.计算卫星在地固系中的直角坐标lX_k=x_k*m.cos(OMEGA_k)-y_k*m.cos(i)*m.sin(OMEGA_k)Y_k=x_k*m.sin(OMEGA_k)+y_k*m.cos(i)*m.cos(OMEGA_k)Z_k=y_k*m.sin(i)#12.计算卫星在协议地球坐标系中的坐标，考虑级移#[X,Y,Z]=[[1,0,X_P],[0,1,-Y_P],[-X_p,Y_P,1]]*[X_k,Y_k,Z_k]if month>12 :  # 恢复历元year = year + 1month = month - 12 print("历元：",year,"年",month,"月",day,"日",hour,"时",minute,"分",second,"秒","卫星PRN号:",PRN,"平均角速度:",n,"卫星平近点角:",M_k,"偏近点角:",E,"真近点角:",V_k,"升交距角:",u_0,"摄动改正项:",δu,δr,δi ,"经摄动改正后的升交距角、卫星矢径和轨道倾角:",u,r,i,"轨道平面坐标X,Y:",x_k,y_k,"观测时刻升交点经度:",OMEGA_k,"地固直角坐标系X:",X_k,"地固直角坐标系Y:",Y_k,"地固直角坐标系Z:",Z_k)t1=None
CulLocation(PRN,year,month,day,hour,minute,second,a_0,a_1,a_2,IDOT,C_rs,δn,M0,C_uc,e,C_us,sqrt_A,TEO,C_ic,OMEGA,C_is,I_0,C_rc,w,OMEGA_DOT,t1)

结果：

注：1，2月视为前一年13，14月,所以最后有一个判断语句，恢复成历元
参考文章：

https://blog.csdn.net/m0_56180742/article/details/120942885
https://blog.csdn.net/weixin_39781323/artice/details/113371222

python读取导航电文并计算卫星位置相关推荐

matlab计算空间坐标,通过matlab计算卫星位置
卫星星历是描述卫星运动轨道的信息.也可以说卫星星历就是一组对应某一时刻的轨道参数及其变率.有了卫星星历就可以计算出任意时刻的卫星位置及其速度.GPS卫星星历分为预报星历和后处理星历.预报星历又称广播星 ...
基于MATLAB计算卫星位置
matlab卫星定位认识星历文件观测文件及位置计算星历文件观测文件 matlab面对对象卫星位置计算卫星高度角方位角计算(待更新...) 认识星历文件观测文件及位置计算初学者,欢迎指正批评. ...
用c++根据轨道六根数计算卫星位置
轨道六根数是描述卫星轨道的一组参数,包括: 轨道长半径(a):卫星轨道的半径,表示卫星到地球中心的平均距离. 轨道偏心率(e):卫星轨道的偏心率,表示轨道的椭圆程度. 轨道倾角(i):卫星轨道与地球赤 ...
python计算csv列平均值_利用Python读取CSV文件并计算某一列的均值和方差
近日需要对excel的csv文件进行处理,求取某银行历年股价的均值方差等一系列数据文件的构成很简单,部分如下所示总共有接近七千行数据,主要的工作就是将其中的股价数据提取出来,放入一个数组之中,然后 ...
python读取csv求平均数_利用Python读取CSV文件并计算某一列的均值和方差
近日需要对excel的csv文件进行处理,求取某银行历年股价的均值方差等一系列数据文件的构成很简单,部分如下所示总共有接近七千行数据,主要的工作就是将其中的股价数据提取出来,放入一个数组之中,然后 ...
根据轨道根数来计算卫星位置
RTKLIB学习总结（六）导航电文、卫星位置计算
文章目录一.导航电文 1.GNSS卫星信号的组成 2.导航电文的编排 3.遥测字(TLW) 4.交接字(HOW) 5.第一数据块 6.第二数据块 7.第三数据块二.卫星钟差钟漂改正 1.时钟校正参 ...
【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）
最近的卫星导航数据处理,老师让我们进行卫星位置的计算,从而使用绘图工具进行对卫星星下点的轨迹进行绘图,这里首先的步骤是读取卫星星历数据,计算卫星位置. 这次的课程目标主要是针对北斗卫星,进行对卫星位置 ...
根据卫星电文计算GPS卫星位置
n文件提供的轨道参数: PS:需要对文件内容格式做一定修改:D改成E/e,各个数据间以空格隔开计算步骤参考武汉大学出版的<GPS测量原理及应用>(第四版)的4.3部分. 使用Java编写 ...

python读取导航电文并计算卫星位置

python简单计算卫星位置

前言

一、思路

二、自动读取导航星历并计算卫星位置

python读取导航电文并计算卫星位置相关推荐

最新文章

热门文章

python读取导航电文并计算卫星位置

python简单计算卫星位置

前言

一、思路

二、 自动读取导航星历并计算卫星位置

python读取导航电文并计算卫星位置相关推荐

最新文章

热门文章

二、自动读取导航星历并计算卫星位置