如何在.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)

所以“数据”将有你想要序列化的格式。不要忘记序列化你想要的标题。