2017-10-18 122 views
-1

我对numpy.polyfit()numpy.polyval()进行了第一次插值,获取了整个卫星轨道的50个经度值。如何在python中执行第二次插值

现在,我只想看一个0-4.5度经度的窗口,并进行第二次插值,以便在窗口中有6,000个经度点。

我需要使用第一次插值的方程/曲线来创建第二个,因为窗口范围只有一个点。我不知道如何做第二次插值。

输入:

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] 

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] 

'Longitude Interpolation' 
fitLon = np.polyfit(fittime, lon, deg) #gets fit coefficients 
polyval_lon = np.polyval(fitLon,fittime) #interp.s to get actual values 


'Get Longitude values for a window of 0-4.5 deg Longitude' 
lonwindow =[] 

for i in range(len(polyval_lon)): 
    if 0 < polyval_lon[i] < 4.5:   # get lon vals in window 
     lonwindow.append(polyval_lon[i]) #append lon vals 

lonwindow = np.array(lonwindow) 

回答

0

首先,生成使用旧的时间(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 
相关问题