np.polyfit y入力用の2D配列にリストを取得する方法

Question

衛星軌道の緯度、経度、高度データを生成しました。今、補間するために私のデータの多項式近似をしたいと思います。 numpy.polyfit()は、1列につき1つのデータセットを含むy座標の2D配列のみを取ります。

現在、経度と高度は2つの別々のリストにあり、それらを2次元配列に結合する必要があります。私はnp.matrix([lon] [alt])を試してみましたが、私はエラーを取得：

TypeError: list indices must be integers, not list

入力：

Lat: [-22.0, -50.5, -5.8, 45.2, 32.7] 

Lon: [-66.6, 21.3, 90.1, 147.4, -115.7] 

Alt: [368752.8, 371184.4, 357834.7, 375131.8, 375643.8] 

所望の出力：

lonAlt_2D_array = ([-66.6, 21.3, 90.1, 147.4, -115.7] [368752.8, 371184.4, 357834.7, 375131.8, 375643.8]) 

fit = np.polyfit(lat, lonAlt_2D_array, 2) #not sure if 2 is correct degree 

全コード：

''' 
Satellite Orbit Propagation Function 
''' 

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) 

    'Generate Orbit' 
    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() 

    '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 
    #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) 


    'Getting the correct J2000 Time' 
    times = [orbitPineapple_J2000time] * propNum 
    #print('times',times) 
    #print('') 

    myOrbitJ2000Time = [] #J2000 times 
    for i in range(propNum): 
     myOrbitJ2000Time.append(times[i] + i) 



    '''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('Lat:') 
    print(lat) 
    print('Lon:') 
    print(lon) 
    print('Alt:') 
    print(alt) 


    lonAlt_2D_array = np.matrix([lon] [alt]) 

    fit = np.polyfit(lat, lonAlt_2D_array, 2) 

orbitPropandcoordTrans（5,2007712.00,0.000939,51.5777,172.5018,323.1066,173.4358,15.68522506）

Answer 1

5x2マトリックスを作成するには、lonAltMatrix = [list(a) for a in zip(lon, alt)]が使用されます。

np.polyfit（）の場合、y座標の行数はx座標の長さと一致する必要があります。