目的
pythonで時系列データの補間処理をしたい。下の図はサンプルプログラムの実行結果。1時間間隔のデータを30分間隔に補間している。
サンプルコード
numpyとscipyのそれぞれの補間関数を使ったサンプルコード。
interp_sample.py
import numpy
from scipy.interpolate import interp1d
from datetime import datetime, timedelta
import math
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
def plot_data(x1, y1, x2, y2, x3, y3):
f, (ax1, ax2, ax3) = plt.subplots(3)
ax1.plot(x1, y1, "o-")
ax2.plot(x2, y2, "o-")
ax3.plot(x3, y3, "o-")
ax1.legend(['original data'])
ax2.legend(['interp data (linear) : numpy'])
ax3.legend(['interp data (cubic) : scipy'])
format = '%H:%M'
ax1.xaxis.set_major_formatter(mdates.DateFormatter(format))
ax2.xaxis.set_major_formatter(mdates.DateFormatter(format))
ax3.xaxis.set_major_formatter(mdates.DateFormatter(format))
# グラフの表示
plt.show()
if __name__ == "__main__":
# サンプル用の元データ作成。ここでは1時間間隔のデータを作る。
t0 = datetime(2019, 6, 1, 12, 0) # 初期時刻
dt = timedelta(minutes=60) # 時間間隔
num = 6 # 要素数
t = [t0 + dt * x for x in range(num)]
x = [math.sin(2 * math.pi * x/(num-1)) for x in range(num)]
# 補間用の時刻列作成。
# ここでは元の時間間隔の半分の時刻列を作る
dt_new = dt/2
num_new = 2*len(t)-1
t_new = [t[0] + x * dt_new for x in range(num_new)]
# interp関数はdatetime型を受け付けない
# 時刻をdatetime型からunix時間（float）に変換する
t_unix = [x.timestamp() for x in t]
t_new_unix = [x.timestamp() for x in t_new]
# 方法1: numpy で補間
x_numpy = numpy.interp(t_new_unix, t_unix, x)
# 方法2: scipy で補間
# 補間方法選択
# kind = "linear", "nearest", "zero", "slinear", "quadratic", "cubic", "previous", "next"
kind = "cubic"
f = interp1d(t_unix, x, kind=kind)
x_scipy = f(t_new_unix)
# データ表示
plot_data(t, x, t_new, x_numpy, t_new, x_scipy)
実行方法
python interp_sample.py
とすると，冒頭の図が出てくる。Python 3.7.3 で動作確認済み。パッケージが足りない場合は，以下のpipで入るはず。
# scipyのインストール（numpyも同時にインストールされるはず）
pip install scipy
# 可視化に使用するmatplotlibをインストール
pip install matplotlib
解説というか自分用メモ
- 補間処理ができるpythonパッケージとして，numpyとscipyが存在した。
- numpyは線形補間ぐらいしかできなさそうだが，scipyのほうは，線形補間以外にも
previous，cubicなど補間方法が豊富でこっちを使うほうがよさそう。 - scipyの
interp1d関数はいわゆる高階関数で，戻り値が関数になっている。この実装は面白いと思った。確かにこっちのほうが使いやすい。 - どっちの補間関数も，横軸（x軸，時間軸）の変数が
datetime型だとエラーが起きた。本サンプルではtimestampによって，いったんunixタイムに変換している。