interp1d插值numpy datetime64时间序列报错解决方案

scipy.interpolate.interp1d 仅支持数值型输入作为插值的x轴，无法直接处理datetime64（即dtype为<M8[ns]的日期时间类型）数组，两次报错的核心原因都是类型转换逻辑错误：

• 第一次直接传入datetime64数组：插值器内部做数值运算时，无法对时间类型和浮点数做算术运算，触发类型错误
• 第二次直接用datetime64除以timedelta64：numpy中只有两个datetime64的差值才是timedelta64类型，直接拿绝对日期时间除以时间差单位，类型不匹配，同样报错

要实现不等间隔时间序列的插值绘图，按以下流程处理即可：

1. 将原始时间序列转为numpy数组，计算所有时间点相对于首个采样点的时间差，再转换为以秒/毫秒等为单位的数值型时间戳，作为插值器的输入x轴
2. 用数值型时间戳和对应测量值创建插值器，可根据需求选择插值类型：linear（线性，默认）、quadratic（二次样条）、cubic（三次样条）；如果是实测数据怕出现不符合实际的过冲，推荐用保形插值PchipInterpolator
3. 生成密度足够高的数值型时间点（仅在原始采样点插值无法体现平滑效果），传入插值器得到平滑后的y值
4. 将高密度数值时间点转回datetime64类型，即可和原始测量点一起绘图

import numpy as np
from scipy.interpolate import interp1d, PchipInterpolator
import matplotlib.pyplot as plt

T = [
    np.datetime64('2020-01-01T00:00:00.000000000'),
    np.datetime64('2020-01-02T00:00:00.000000000'),
    np.datetime64('2020-01-03T00:00:00.000000000'),
    np.datetime64('2020-01-05T00:00:00.000000000'),
    np.datetime64('2020-01-06T00:00:00.000000000'),
    np.datetime64('2020-01-09T00:00:00.000000000'),
    np.datetime64('2020-01-13T00:00:00.000000000'),
]
Z = [543, 234, 435, 765, 564, 235, 345]

# 转numpy数组方便向量计算
T = np.array(T)
Z = np.array(Z)

# 时间类型转数值：计算距首个采样点的秒数
t0 = T[0]
t_numeric = (T - t0) / np.timedelta64(1, 's')

# 创建插值器，这里用三次样条，替换成PchipInterpolator(t_numeric, Z)可实现无过冲保形插值
interp_func = interp1d(t_numeric, Z, kind='cubic')

# 生成高密度插值点：这里步长设为3600秒（1小时），可根据时间范围调整密度
t_dense_num = np.arange(t_numeric.min(), t_numeric.max(), 3600)
# 数值时间戳转回datetime64用于绘图
T_dense = t0 + t_dense_num * np.timedelta64(1, 's')
Z_dense = interp_func(t_dense_num)

# 绘图
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot()
ax.plot(T, Z, 'o', label='原始精确测量点')
ax.plot(T_dense, Z_dense, '-', label='三次样条插值曲线')
ax.legend()
ax.set_xlabel('采样时间')
ax.set_ylabel('测量值Z')
plt.gcf().autofmt_xdate() # 自动旋转x轴时间标签，避免重叠
plt.show()

• 不要仅在原始采样点位置计算插值，否则插值结果和原始点完全重合，和直线连接效果无区别，必须生成更高密度的采样点才能得到平滑曲线
• 三次样条插值在数据波动较大时可能出现超出原始值范围的过冲现象，对物理测量、工程数据推荐使用PchipInterpolator，可以保证曲线经过所有原始点的同时，不会出现不符合实际的过冲
• 高密度点的步长可根据数据总时长调整：跨年级别的数据可以用1天为步长，跨小时级别的数据可以用1分钟为步长，平衡平滑度和渲染性能

内容的提问来源于stack exchange，提问作者vdaiep