如何在.dat文件中寫入數據Python新行多個列表衛星軌道

Question

我創建了一個函數來生成和傳播衛星軌道。現在，我想將所有數據保存在.dat文件中。我不確定需要多少循環，或者完全可以如何操作。

我希望每個傳播步驟的時間，緯度，經度和高度都在一行上。

代碼數據：

myOrbitJ2000Time = [1085.0, 2170.0, 3255.0, 4340.1, 5425.1] 

lat = [48.5, 26.5, -28.8, -48.1, 0.1] 

lon = [-144.1, -50.4, -1.6, 91.5, 152.8] 

alt = [264779.5, 262446.1, 319661.8, 355717.3, 306129.0] 

在.dat文件所需的輸出：

J2000時間，緯度，經度，Alt鍵

1085.0, 48.6, -144.2, 264779.5 

2170.0, 26.5, -50.4, 262446.1 

3255.0, -28.7, -1.6, 319661.8 

4340.1, -48.0, 91.5, 355717.2 

5425.1, 0.1, 152.8, 06129.0 

全碼：

個

import orbital 
from orbital import earth, KeplerianElements, plot 
import matplotlib.pyplot as plt 
import numpy as np 
from astropy import time 
from astropy.time import TimeDelta, Time 
from astropy import units as u 
from astropy import coordinates as coord 


#def orbitPropandcoordTrans(orbitPineapple_J2000time, _ecc, _inc, _raan, _arg_pe, _meananom, meanMotion): 
def orbitPropandcoordTrans(propNum, orbitPineapple_J2000time, _ecc, _inc, _raan, _arg_pe, _meananom, meanMotion): 
     ''' 
     Create original orbit and run for 100 propagations (in total one whole orbit) 
     in order to get xyz and time for each propagation step. 
     The end goal is to plot the lat, lon, & alt data to see if it matches ISS groundtrace. 
     ''' 
    import orbital 
    from orbital import earth, KeplerianElements, plot 
    import matplotlib.pyplot as plt 
    import numpy as np 
    from astropy import time 
    from astropy.time import TimeDelta, Time 
    from astropy import units as u 
    from astropy import coordinates as coord 

    'Calculate Avg. Period from Mean Motion' 
    _avgPeriod = 86400/meanMotion 
    #print('_avgPeriod', _avgPeriod) 

    ###### 
    #propNum = int(_avgPeriod) 

    'Generate Orbit' 
    #orbitPineapple = KeplerianElements.with_period(5560, body=earth, e=0.0004525, i=(np.deg2rad(51.6414)), raan=(np.deg2rad(247.1662))) 
    orbitPineapple = KeplerianElements.with_period(_avgPeriod, body=earth, e=_ecc, i=(np.deg2rad(_inc)), raan=(np.deg2rad(_raan)), arg_pe=(np.deg2rad(_arg_pe)), M0=(np.deg2rad(_meananom))) #ref_epoch= 
    #plot(orbitPineapple) 
    #plt.show() 
    #print('') 
    #print('') 

    'Propagate Orbit and retrieve xyz' 
    myOrbitX = []   #X Coordinate for propagated orbit step 
    myOrbitY = []   #Y Coordinate for propagated orbit step 
    myOrbitZ = []   #Z Coordinate for propagated orbit step 
    myOrbitTime = []  #Time for each propagated orbit step 
    myOrbitJ2000Time = [] #J2000 times 
    #propNum = 100  #Number of propagations and Mean Anomaly size (one orbit 2pi/propNum) 

    for i in range(propNum): 
     orbitPineapple.propagate_anomaly_by(M=(2.0*np.pi/propNum)) #Propagate the orbit by the Mean Anomaly 
     myOrbitX.append(orbitPineapple.r.x)      #x vals 
     myOrbitY.append(orbitPineapple.r.y)      #y vals 
     myOrbitZ.append(orbitPineapple.r.z)      #z vals 
     myOrbitTime.append(orbitPineapple_J2000time)    #J2000 time vals 
     myOrbitJ2000Time.append(orbitPineapple.t) 

     #plot(orbitPineapple) 
    #print('X:', 'Length:', len(myOrbitX)) 
    #print(myOrbitX) 
    #print('Y:', 'Length:',len(myOrbitY)) 
    #print(myOrbitY) 
    #print('Z:', 'Length:', len(myOrbitZ)) 
    #print(myOrbitZ) 
    #print('J2000 Time:', 'Length:',len(myOrbitTime)) 
    #print(myOrbitTime) 


    '''Because the myOrbitTime is only the time between each step to be the sum of itself plus 
    all the previous times. And then I need to convert that time from seconds after J2000 to UTC.''' 
    myT = [] #UTC time list 

    for i in range(propNum): 
     myT.append((Time(2000, format='jyear') + TimeDelta(myOrbitTime[i]*u.s)).iso) #Convert time from J2000 to UTC 
    #print('UTC Time List Length:', len(myT)) 
    #print('UTC Times:', myT) 

    '''Now I have xyz and time for each propagation step and need to convert the coordinates from 
    ECI to Lat, Lon, & Alt''' 
    now = []  #UTC time at each propagation step 
    xyz =[]  #Xyz coordinates from OrbitalPy initial orbit propagation 
    cartrep = [] #Cartesian Representation 
    gcrs = [] #Geocentric Celestial Reference System/Geocentric Equatorial Inertial, the default coord system of OrbitalPy 
    itrs =[]  #International Terrestrial Reference System coordinates 
    lat = []  #Longitude of the location, for the default ellipsoid 
    lon = []  #Longitude of the location, for the default ellipsoid 
    alt = []  #Height of the location, for the default ellipsoid 


    for i in range(propNum): 
     xyz = (myOrbitX[i], myOrbitY[i], myOrbitZ[i])     #Xyz coord for each prop. step 
     now = time.Time(myT[i])           #UTC time at each propagation step 
     cartrep = coord.CartesianRepresentation(*xyz, unit=u.m)   #Add units of [m] to xyz 
     gcrs = coord.GCRS(cartrep, obstime=time.Time(myT[i]))   #Let AstroPy know xyz is in GCRS 
     itrs = gcrs.transform_to(coord.ITRS(obstime=time.Time(myT[i]))) #Convert GCRS to ITRS 
     loc = coord.EarthLocation(*itrs.cartesian.xyz)     #Get lat/lon/height from ITRS 
     lat.append(loc.lat.value)          #Create latitude list 
     lon.append(loc.lon.value)          #Create longitude list 
     alt.append(loc.height.value)   

    #print('Current Time:', now) 
    #print('') 
    #print('UTC Time:') 
    #print(myT) 
    print('myOrbitJ2000Time', myOrbitJ2000Time) 
    print('') 
    print('Lat:') 
    print('') 
    print(lat) 
    print('') 
    print('Lon:') 
    print(lon) 
    print('') 
    print('Alt:') 
    print(alt) 

orbitPropandcoordTrans（5，-34963095，0.0073662，51.5946，154.8079，103.6176，257.3038，15.92610159）

Answer 1

要回答您的一般問題，不要定義一堆Python列表（這些列表很慢添加和使用，特別是在您擴大解決方案時處理大量值），您可以先開始創建Numpy數組而不是列表。初始化Numpy數組時，通常需要指定數組的大小，但在這種情況下，這很容易，因爲您知道需要多少傳播。例如：

>>> orbitX = np.empty(propNum, dtype=float) 

創建的propNum空numpy的陣列浮點值（此處「空」指的是陣列只是未初始化的任何值 - 這是創建一個新的數組，因爲我們的最快方法「重新去填補它的所有值以後反正

然後在你的循環，而不是使用orbitX.append(x)你將會分配給數組中對應於當前刻度值：。orbitX[i] = x同爲其他案件

。

然後有幾種可能性，如何出將放入您的數據中，但使用Astropy Table可提供很大的靈活性。你可以創建一個包含你想輕鬆地像列的表：

>>> from astropy.table import Table 
>>> table = Table([J2000, lat, lon, alt], names=('J2000', 'lat', 'lon', 'alt')) 

關於具有表對象的好處是，有一噸的輸出格式選項。例如。在Python提示符下：

>>> table 
    J2000   lat   lon   alt 
    float64  float64  float64  float64 
------------- -------------- -------------- ------------- 
1085.01128806 48.5487129749 -144.175054697 264779.500823 
2170.02257613 26.5068122883 -50.3805485685 262446.085716 
3255.03386419 -28.7915478311 -1.6090285674 319661.817451 
4340.04515225 -48.0536526356 91.5416828221 355717.274021 
5425.05644032 0.084422392655 152.811717713 306129.02576 

要輸出到文件，首先必須考慮如何將數據格式化。目前已經有很多常用的數據格式可供您考慮使用，但這取決於數據的內容以及將要使用的格式（「。dat文件「並不意味着什麼文件格式;或者說，它可能意味着任何東西）。但在你的問題中，你指出你想要什麼是所謂的」逗號分隔值「（CSV），其中的數據被格式化的列下去的時候，用逗號分開的行內的每個值的Table類可以輸出CSV（和任何變體）很容易地：。

>>> import sys 
>>> table.write(sys.stdout, format='ascii.csv') # Note: I'm just using sys.stdout for demonstration purposes; normally you would give a filename 
J2000,lat,lon,alt 
1085.011288063746,48.54871297493748,-144.17505469691633,264779.5008225624 
2170.022576127492,26.506812288280788,-50.38054856853237,262446.0857159357 
3255.0338641912376,-28.79154783108773,-1.6090285674024463,319661.81745081506 
4340.045152254984,-48.05365263557444,91.54168282208444,355717.2740210131 
5425.05644031873,0.08442239265500713,152.81171771323176,306129.0257600865 

有許多雖然其它選項例如，如果你想要的數據很好地格式化在對齊的列中，你也可以這樣做，你可以閱讀更多有關它here。（另外，我建議如果你想要一個純文本文件格式，你可以嘗試ascii.ecsv，它的優點是它輸出addit有理的元數據，允許它很容易讀回的Astropy表）：

>>> table.write(sys.stdout, format='ascii.ecsv') 
# %ECSV 0.9 
# --- 
# datatype: 
# - {name: J2000, datatype: float64} 
# - {name: lat, datatype: float64} 
# - {name: lon, datatype: float64} 
# - {name: alt, datatype: float64} 
# schema: astropy-2.0 
J2000 lat lon alt 
1085.01128806 48.5487129749 -144.175054697 264779.500823 
2170.02257613 26.5068122883 -50.3805485685 262446.085716 
3255.03386419 -28.7915478311 -1.6090285674 319661.817451 
4340.04515225 -48.0536526356 91.5416828221 355717.274021 
5425.05644032 0.084422392655 152.811717713 306129.02576 

我會注意的另一個無關的事情是，在Astropy許多對象可以在整個陣列的操作，除了單一的價值觀，而這通常可以更有效率，尤其是對於許多價值。特別是，你有這樣的Python循環：

for i in range(propNum): 
    xyz = (myOrbitX[i], myOrbitY[i], myOrbitZ[i])     #Xyz coord for each prop. step 
    now = time.Time(myT[i])           #UTC time at each propagation step 
    cartrep = coord.CartesianRepresentation(*xyz, unit=u.m)   #Add units of [m] to xyz 
    gcrs = coord.GCRS(cartrep, obstime=time.Time(myT[i]))   #Let AstroPy know xyz is in GCRS 
    itrs = gcrs.transform_to(coord.ITRS(obstime=time.Time(myT[i]))) #Convert GCRS to ITRS 
    loc = coord.EarthLocation(*itrs.cartesian.xyz)     #Get lat/lon/height from ITRS 
    lat.append(loc.lat.value)          #Create latitude list 
    lon.append(loc.lon.value)          #Create longitude list 
    alt.append(loc.height.value) 

這可以將它們視爲陣列座標而不是完全重寫，沒有明確的循環。例如：

>>> times = time.Time(myT) # All the times, not just a single one 
>>> cartrep = coord.CartesianRepresentation(orbitX, orbitY, orbitZ, unit=u.m) # Again, an array of coordinates 
>>> gcrs = coord.GCRS(cartrep, obstime=times) 
>>> itrs = gcrs.transform_to(coord.ITRS(obstime=ts)) 
>>> loc = coord.EarthLocation(*itrs.cartesian.xyz) # I'm not sure this is the most efficient way to do this but I'm not an expert on the coordinates package 

這樣我們就可以獲得所有的座標數組了。例如：

>>> loc.lat 
<Latitude [ 48.54871297, 26.50681229,-28.79154783,-48.05365264, 
      0.08442239] deg> 

所以這通常是更有效的方式來處理大量的座標值。同樣，在原始代碼中轉換myT時，您可以創建一個單一的TimeDelta數組，並將其添加到您的初始時間，而不是循環遍歷所有時間偏移量。

不幸的是，我不是orbital軟件包的專家，但它似乎並不像一個人想要在軌道上的不同點上獲得座標那樣容易。 ISTM就像應該有的那樣。

Answer 2

也許最簡單的方法是使用拉鍊（）：

data = zip(myOrbitJ2000Time, lat, lon, alt)

所以「數據」將有你想要序列化的格式。不要忘記序列化你想要的標題。