SpaceX 公司發射starlink衛星243顆星鏈衛星，衛星會很大程度影響天文觀測。但是有愛好者又想觀測到starlink衛星過境，所以我寫了以下代碼方便大家使用。

需要準備python的輪子包括：urllib，sgp4，spiceypy。

首先從celestrak下載starlink兩行報，并將兩行報保存在列表嵌套的字典中（不知道衛星兩行報的自行百度）。

然后指定衛星名稱，讀取兩行報文件進行預報（或者遍歷全部,遍歷的話需要并行處理。）。但是由于python模塊SGP4輸出位置位于TEME參考架，需要轉換至J2000進而轉換至WGS84坐標系（經度緯度坐標系），就得到衛星的星下點坐標。如果要計算某一位置的starlink觀測的方位角和天頂角需要計算衛星與此處的相對位置關系。（這一節先計算到衛星星下點，有需要的可以留言.... frame.py文件為坐標轉換文件。見后續blog說明）

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 23 19:54:02 2020

@author: panda
"""
import urllib.request
import os
import platform
from datetime import *
from sgp4.earth_gravity import wgs72,wgs84
from sgp4.io import twoline2rv
import frame
import spiceypy as spice
from math import atan2, cos, pi, sin, sqrt, tan
import numpy as np
def main():

? ? #星歷文件
? ? f_KERNELS='/home/panda/SPICE/gen/'
? ? spice.furnsh(f_KERNELS + 'pck/pck00010.tpc')
? ? spice.furnsh(f_KERNELS+'lsk/naif0011.tls')
? ? time_intev=60*60#時間間隔,單位：s
? ? loop=False#True#loop satellite or not（True=會循環所有starlink衛星）
? ? #star_name='STARLINK-53'#star_link_name
? ? rad2deg=180./pi
? ? file_nut80 = 'nut80.dat'地球旋轉矩陣
? ? url='http://www.celestrak.com/NORAD/elements/starlink.txt'
? ? tle_path='tle'
? ? result_path='result'
? ? today=datetime.utcnow().strftime('%Y%m%d')#time utc
? ? sys_p = platform.system()
? ? if sys_p == "Windows":
? ? ? ? print("OS is Windows!!!")
? ? ? ? path_tle=os.getcwd()+'\\'+tle_path+'\\'
? ? ? ? path_result=os.getcwd()+'\\'+result_path+'\\'
? ? elif sys_p == "Linux":
? ? ? ? print("OS is Linux!!!")
? ? ? ? path_tle=os.getcwd()+'/'+tle_path+'/'
? ? ? ? path_result=os.getcwd()+'/'+result_path+'/'+today+'/'
? ? else:
? ? ? ? print(sys_p)
? ? ? ? path_tle=os.getcwd()+'/'+tle_path+'/'
? ? ? ? path_result=os.getcwd()+'/'+result_path+'/'+today+'/'
? ??
? ? print(today)
? ? data = urllib.request.urlopen(url).read()
? ? save_tle(data,today,path_tle)?
? ? tles=split_tle((data))
? ? try:
? ? ? ? if not os.path.exists(path_result):
? ? ? ? ? ? os.makedirs(path_result)
? ? except IOError as e:
? ? ? ? print(e)
? ? for tle in tles:
? ? ? ? #print(tle['name'])
? ? ? ? #if tle['name']==star_name:
? ? ? ? new_list = [ x for x in tle['line1'].split(' ') if x != '' ]
? ? ? ? time_y=float(new_list[3])
? ? ? ? year_tle=str(int(time_y/1000)+2000)
? ? ? ? day_lte=(time_y-(int(time_y/1000)*1000))
? ? ? ? utc=year_tle+'-01-01T00:00:00'
? ? ? ? et=spice.spiceypy.utc2et(utc)+(day_lte-1)*24.*60.*60. ??
? ? ? ? end_et=et+24.*60.*60.*15
? ? ? ? file_name=path_result+tle['name']+'.txt'
? ? ? ? print(tle['name'])
? ? ? ? with open(file_name, 'w+') as fh:
? ? ? ? ? ? while et<= end_et:
? ? ? ? ? ? ? ? utc=spice.spiceypy.et2utc(et,'ISOC',0,24)

? ? ? ? ? ? ? ? #得到衛星地面星下點，的經緯度
? ? ? ? ? ? ? ? Earth2Sat, v_Earth2Sat,lat1,lon1,alt1=satellite_j2000(utc,et,tle['name'],tle['line1'],tle['line2'])
? ? ? ? ? ? ? ? #print(type(utc),type(lat1),type(lon1))
? ? ? ? ? ? ? ? fh.write('%s %-8.2f %-8.2f\n'%(str(utc),lat1,lon1))
? ? ? ? ? ? ? ? et = et + time_intev

? ? print(len(tles))
? ?
? ? spice.kclear()
def save_tle(data,today,path):
? ??
? ? #today=datetime.utcnow().strftime('%Y%m%d')#time utc
? ? file_name=path+today+'.txt'# filename
? ? try:
? ? ? ? if not os.path.exists(path):
? ? ? ? ? ? os.makedirs(path)
? ? ? ? if not os.path.exists(file_name):#if file not exist load file from url
? ? ? ? ? ? #url='http://www.celestrak.com/NORAD/elements/starlink.txt'
? ? ? ? ? ? #data = urllib.request.urlopen(url).read()
? ? ? ? ? ? with open(file_name, 'w+b') as fh:
? ? ? ? ? ? ? ? fh.write(data)
? ? ? ? else:
? ? ? ? ? ? fsize = os.path.getsize(file_name)
? ? ? ? ? ? if fsize==0:
? ? ? ? ? ? ? ? print('reloading:',file_name)
? ? ? ? ? ? ? ? #url='http://www.celestrak.com/NORAD/elements/starlink.txt'
? ? ? ? ? ? ? ? #data = urllib.request.urlopen(url).read()
? ? ? ? ? ? ? ? with open(file_name, 'w+b') as fh:
? ? ? ? ? ? ? ? ? ? fh.write(data) ? ? ?
? ? except IOError as e:
? ? ? ? print(e)
def split_tle(data):
? ? tle={}
? ? list_data=data.decode('utf-8').split('\n')
? ? tles=[]
? ? for i in range(len(list_data)-1):
? ? ? ? if list_data[i][0]=='1' or list_data[i][0]=='2':
? ? ? ? ? ? x=1
? ? ? ? else:
? ? ? ? ? ? #print(len(list_data[i]))
? ? ? ? ? ? tle['name']=list_data[i].strip()?
? ? ? ? ? ? tle['line1']=list_data[i+1].strip()?
? ? ? ? ? ? tle['line2']=list_data[i+2].strip()?
? ? ? ? ? ? tles.append(tle)
? ? ? ? ? ? tle={}
? ? return(tles)

def load_tle(url):
? ? mode = re.compile(r'\d+')
? ? data = urllib.request.urlopen(url).read()#
? ? data = data.decode('UTF-8')
? ? tle_lines=data.splitlines()
? ? for i in range(len(tle_lines)):
? ? ? ? line=tle_lines[i]
? ? ? ? if op.eq(line.strip(),celestrak_flag):
? ? ? ? ? ? TlE_name=(line)
? ? ? ? ? ? lines1=(tle_lines[i+1])
? ? ? ? ? ? lines2=(tle_lines[i+2])
? ? return(TlE_name,lines1,lines2)
def satellite_j2000(utc,et,TlE_name,lines1,lines2):
? ? yr=int(utc[0:4]) ? ? ? ? ?
? ? mon=int(utc[5:7]) ? ? ? ? ??
? ? day=int(utc[8:10]) ? ? ? ? ??
? ? hr=int(utc[11:13]) ? ? ? ??
? ? mins=int(utc[14:16]) ? ? ? ? ??
? ? sec=int(utc[17:19]) ?
? ? satellite = twoline2rv(lines1, lines2, wgs84)
? ? reci1,veci1= satellite.propagate(yr, mon, day, hr, mins, sec)
? ? jdtt, jdttfrac = jday(yr, mon, day, hr, mins, sec)
? ? ttt = (jdtt - 2451545.0) / 36525.0
? ? opt = '80'
? ? order = 106
? ? eqeterms = 2
? ? ateme = [0, 0, 0]
? ? file_nut80 = 'nut80.dat'
? ? reci, veci, aeci = frame.teme2eci(reci1, veci1, ateme, ttt, order, eqeterms, opt, file_nut80)
? ? lat1,lon1,alt1 = calculate(reci,et)
? ? return(reci,veci,lat1,lon1,alt1)
def calculate(options,et):
? ? x = options[0]#*1000.0
? ? y = options[1]#*1000.0
? ? z = options[2]#*1000.0
? ? xtrn=spice.spiceypy.pxform('J2000','IAU_EARTH',et)
? ? jstate=np.dot(options,np.transpose(xtrn))
? ? raddii=spice.spiceypy.bodvar(399,"RADII", 3)
? ? re=raddii[0]
? ? rp=raddii[2]
? ? flat=(re-rp)/re
? ? lon,lat,alt=spice.spiceypy.recpgr('EARTH',jstate,re,flat)
? ? lat = ((lat*180.)/pi)
? ? lon = ((lon*180.)/pi)
? ? return (lat, lon, alt)
def jday(yr, mon, day, hr, mins, sec):
? ? jd = 367.0 * yr - int( (7 * (yr +int( (mon + 9) / 12.0) ) ) * 0.25 ) + int( 275 * mon / 9.0 )+ day + 1721013.5
? ? jdfrac = (sec + mins * 60.0 + hr *3600.0) / 86400.0
? ? if jdfrac > 1.0:
? ? ? ? jd = jd + int(jdfrac)
? ? ? ? jdfrac = jdfrac - int(jdfrac)
? ? return(jd,jdfrac)
if __name__ == "__main__":
? ? main()

總結

以上是生活随笔為你收集整理的starlink卫星轨道预报的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。