2017-10-17 146 views
0

我做了numpy.polyfit()緯度,經度,&海拔高度數據的衛星軌道和插值(50分)與numpy.polyval()如何內插numpy.polyval和numpy.polyfit python

現在,我只需要一個窗口(0-4.5度經度)並進行更高分辨率插值(6000點)。我認爲我需要使用第一個低分辨率擬合的擬合係數來插入我的經度窗口,而我不太確定如何做到這一點。

輸入:

lat = [27.755611104020687, 22.50661883405905, 17.083576087905502, 11.53891099628959, 5.916633366002468, 0.2555772624429494, -5.407902834141322, -11.037514984810027, -16.594621304857206, -22.03556688048686, -27.308475759820045, -32.34927891621322, -37.07690156937186, -41.38803163295967, -45.15306971601912, -48.21703193866987, -50.41165326774015, -51.58419672864487, -51.63883932997542, -50.57025116952513, -48.46557920053242, -45.47329014246061, -41.76143266388077, -37.48707787049647, -32.782653540783, -27.754184631685046, -22.48503337048438, -17.041097574740743, -11.475689837873944, -5.833592289780744, -0.1543286595142316, 5.525119007560692, 11.167878192881306, 16.73476477885508, 22.18160021405449, 27.455997555900108, 32.493386953033685, 37.21222272985329, 41.508824407948275, 45.25350232626601, 48.291788915858554, 50.45698534747271, 51.59925055739275, 51.62660832560593, 50.53733379179681, 48.420673231121725, 45.42531420150485, 41.71819693220144, 37.45473807165676, 32.76569228387106] 

lon = [-109.73105744378498, -104.28690174554579, -99.2435132929552, -94.48533149079628, -89.91054414962821, -85.42671400689177, -80.94616150449806, -76.38135021210172, -71.6402674905218, -66.62178379632216, -61.21120467960157, -55.27684029674759, -48.66970878028004, -41.23083703244677, -32.813881865289346, -23.332386757370532, -12.832819226213942, -1.5659455609661785, 10.008077792630402, 21.33116444634303, 31.92601575632583, 41.51883213364072, 50.04498630545507, 57.58103957109249, 64.26993028992476, 70.2708323505337, 75.73441871754586, 80.7944079829813, 85.56734813043659, 90.1558676264546, 94.65309120129724, 99.14730128118617, 103.72658922048785, 108.48349841714494, 113.51966824008079, 118.95024882101737, 124.9072309203375, 131.5395221402974, 139.00523971191907, 147.44847902856114, 156.95146022590976, 167.46163867248032, 178.72228750873975, -169.72898181991064, -158.44642409799974, -147.8993300787564, -138.35373014113995, -129.86955508919888, -122.36868103811106, -115.70852432245486] 

alt = [374065.49207488785, 372510.1635949105, 371072.75959230476, 369836.3092635453, 368866.7921820211, 368209.0950216997, 367884.3703536549, 367888.97894243425, 368195.08833668986, 368752.88080031495, 369494.21701128664, 370337.49662954226, 371193.3839051864, 371971.0136622536, 372584.272228585, 372957.752022573, 373032.0104747458, 372767.8112563471, 372149.0940816824, 371184.49208500446, 369907.2992362557, 368373.8795969478, 366660.5935723809, 364859.4071422184, 363072.42955020745, 361405.69765685993, 359962.58417682414, 358837.24421522504, 358108.5277743581, 357834.7679493668, 358049.8054538341, 358760.531463618, 359946.1257064284, 361559.04646970675, 363527.70518032915, 365760.6377191965, 368151.8843206526, 370587.2165838985, 372950.8014553002, 375131.8814988529, 377031.06540952163, 378565.8596562773, 379675.13241518533, 380322.2707576381, 380496.8682141012, 380214.86538256245, 379517.14674525027, 378466.68079100474, 377144.36811517406, 375643.83731560566] 

myOrbitJ2000Time =[ 20027712., 20027713., 20027714., 20027715., 20027716., 
     20027717., 20027718., 20027719., 20027720., 20027721., 
     20027722., 20027723., 20027724., 20027725., 20027726., 
     20027727., 20027728., 20027729., 20027730., 20027731., 
     20027732., 20027733., 20027734., 20027735., 20027736., 
     20027737., 20027738., 20027739., 20027740., 20027741., 
     20027742., 20027743., 20027744., 20027745., 20027746., 
     20027747., 20027748., 20027749., 20027750., 20027751., 
     20027752., 20027753., 20027754., 20027755., 20027756., 
     20027757., 20027758., 20027759., 20027760., 20027761.] 

代碼:

deg = 30 #polynomial degree for fit 
fittime = myOrbitJ2000Time - myOrbitJ2000Time[0] 

'Latitude Interpolation'  
fitLat = np.polyfit(fittime, lat, deg) 
polyval_lat = np.polyval(fitLat,fittime) 

'Longitude Interpolation' 
fitLon = np.polyfit(fittime, lon, deg) 
polyval_lon = np.polyval(fitLon,fittime) 

'Altitude Interpolation' 
fitAlt = np.polyfit(fittime, alt, deg) 
polyval_alt = np.polyval(fitAlt,fittime) 


'Get Lat, Lon, & Alt values for a window of 0-4.5 deg Longitude' 
lonwindow =[] 
latwindow = [] 
altwindow = [] 
for i in range(len(polyval_lat)): 
    if 0 < polyval_lon[i] < 4.5:   # get lon vals in window 
     lonwindow.append(polyval_lon[i]) #append lon vals 
     latwindow.append(polyval_lat[i]) #append corresponding lat vals 
     altwindow.append(polyval_alt[i]) #append corresponding alt vals 

lonwindow = np.array(lonwindow) 

只要是明確的 - 問題是我只有在窗口範圍內的一個點,我想用插值/等式/來自上一步的曲線。那麼我可以使用它再次插值並在我的窗口範圍內生成6000個點。

+0

您對'myOrbitJ2000Time'使用了什麼值?我試圖猜測'np.arange(50)',我不知道是不是因爲這個原因,但是它說這個polyfit可能條件不好。 – Reti43

+0

抱歉,忘記了。我用時間數據更新了我的問題。它總是說,如果你使用一個高於〜10的數字作爲度數輸入,polyfit可能條件不好。 – Rose

+0

因爲你的'myOrbitJ2000Time'值是連續的,所以'fittime'最終變成'np.arange(50)'。更重要的是,我剛剛意識到你的度數是一個角度值,而polyfit的度數是指多項式的度數。你究竟想要適合什麼?你的'lat'和'lon'值是以度數表示的。 – Reti43

回答

0

原來的答覆張貼here

首先,使用舊的時間(x軸)的值,和內插的經度(y軸)的值生成多項式擬合係數。

import numpy as np 
import matplotlib.pyplot as plt 

poly_deg = 3 #degree of the polynomial fit 
polynomial_fit_coeff = np.polyfit(original_times, interp_lon, poly_deg) 

接下來,使用np.linspace()根據窗口中慾望點的數量生成任意時間值。

start = 0 
stop = 4 
num_points = 6000 
arbitrary_time = np.linspace(start, stop, num_points) 

最後,使用擬合係數和任意時間來獲取實際插值經度(y軸)值和圖。

lon_intrp_2 = np.polyval(polynomial_fit_coeff, arbitrary_time) 

plt.plot(arbitrary_time, lon_intrp_2, 'r') #interpolated window as a red curve 
plt.plot(myOrbitJ2000Time, lon, '.') #original data plotted as points