如何在 .dat 文件中写入数据 Python new line multiple lists Satellite Orbit
How to write data in a .dat file Python new line multiple lists Satellite Orbit
我创建了一个函数来生成和传播卫星轨道。现在，我想将所有数据保存在一个 .dat 文件中。我不确定我需要多少个 for 循环，也不知道该怎么做。
我希望每个传播步骤的时间、纬度、经度和高度都在一条线上。
数据代码：
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)
可能最简单的方法是使用 zip():
data = zip(myOrbitJ2000Time, lat, lon, alt)
因此 "data" 将具有您要序列化的格式。不要忘记序列化你想要的header。
回答您的一般性问题，而不是定义一堆 Python 列表（追加和使用它们很慢，尤其是在您扩大分辨率时处理大量值时） ，您可以从创建 Numpy 数组而不是列表开始。在初始化 Numpy 数组时，通常需要指定数组的大小，但在这种情况下这很容易，因为您知道需要多少次传播。例如：
>>> orbitX = np.empty(propNum, dtype=float)
创建一个包含 propNum 个浮点值的空 Numpy 数组（此处 "empty" 表示该数组未使用任何值进行初始化——这是创建新数组的最快方法无论如何，我们稍后会填充其中的所有值。
然后在您的循环中而不是使用 orbitX.append(x)，您将分配给数组中对应于当前报价的值：orbitX[i] = x。其他情况也一样。
那么如何输出你的数据有几种可能性，但是使用 Astropy Table 提供了很大的灵活性。您可以创建一个包含您想要的列的 table，例如：
>>> from astropy.table import Table
>>> table = Table([J2000, lat, lon, alt], names=('J2000', 'lat', 'lon', 'alt'))
拥有 Table 对象的好处是有大量的输出格式选项。例如。在 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 文件”在文件格式方面没有任何意义；或者更确切地说，它可以意思是任何东西）。但是在您的问题中，您指出您想要的是所谓的 "comma-separated values" (CSV)，其中数据的格式为列向下，一行中的每个值用逗号分隔。 Table class 可以很容易地输出 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，它的优点是它输出额外的元数据，可以很容易地读回 Astropy Table）：
>>> 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 应该有。
我创建了一个函数来生成和传播卫星轨道。现在，我想将所有数据保存在一个 .dat 文件中。我不确定我需要多少个 for 循环，也不知道该怎么做。
我希望每个传播步骤的时间、纬度、经度和高度都在一条线上。
数据代码：
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)
可能最简单的方法是使用 zip():
data = zip(myOrbitJ2000Time, lat, lon, alt)
因此 "data" 将具有您要序列化的格式。不要忘记序列化你想要的header。
回答您的一般性问题，而不是定义一堆 Python 列表（追加和使用它们很慢，尤其是在您扩大分辨率时处理大量值时） ，您可以从创建 Numpy 数组而不是列表开始。在初始化 Numpy 数组时，通常需要指定数组的大小，但在这种情况下这很容易，因为您知道需要多少次传播。例如：
>>> orbitX = np.empty(propNum, dtype=float)
创建一个包含 propNum 个浮点值的空 Numpy 数组（此处 "empty" 表示该数组未使用任何值进行初始化——这是创建新数组的最快方法无论如何，我们稍后会填充其中的所有值。
然后在您的循环中而不是使用 orbitX.append(x)，您将分配给数组中对应于当前报价的值：orbitX[i] = x。其他情况也一样。
那么如何输出你的数据有几种可能性，但是使用 Astropy Table 提供了很大的灵活性。您可以创建一个包含您想要的列的 table，例如：
>>> from astropy.table import Table
>>> table = Table([J2000, lat, lon, alt], names=('J2000', 'lat', 'lon', 'alt'))
拥有 Table 对象的好处是有大量的输出格式选项。例如。在 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 文件”在文件格式方面没有任何意义；或者更确切地说，它可以意思是任何东西）。但是在您的问题中，您指出您想要的是所谓的 "comma-separated values" (CSV)，其中数据的格式为列向下，一行中的每个值用逗号分隔。 Table class 可以很容易地输出 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，它的优点是它输出额外的元数据，可以很容易地读回 Astropy Table）：
>>> 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 应该有。