光谱响应函数线性内插

0 前言

国内卫星数据，官方给的光谱响应函数一般为间隔1nm的数据，但有些大气校正模型可能需要其他间隔的数据，比如6s需要2.5nm间隔的数据，因此需要对该数据进行插值处理。当然可以实现该功能的方法千千万，本文尝试用python进行线性插值处理，有其它方法欢迎留言补充。

1 内容

1.1 光谱响应函数

国内卫星大部分的光谱响应函数可通过中国资源卫星中心进行下载，下载地址

1.2 np.interp()插值处理

1.2.1 np.interp()简介

numpy 中有个函数interp可实现线性内插，下面简单介绍主要参数，具体内容参考官方资料

numpy.interp(x, xp, fp, left=None, right=None, period=None)
主要参数介绍：
x为待插数据的x轴，xp为原始数据的x轴数据，fp为原始数据的y轴数据

1.2.2 代码

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

wavelength = pd.read_excel('GF.xlsx', 'PMS1')
wavelength.head()

length_interp_bin = np.arange(450, 891, 2.5)  #截距段
x = wavelength.iloc[:,0].tolist()
wavelength_interp = []
for i in range(2,6):
    wavelength_interp_temp = np.interp(length_interp_bin, x, wavelength.iloc[:,i])
    wavelength_interp.append(wavelength_interp_temp)
    
wavelength_interp = np.array(wavelength_interp)

plt.figure(figsize=(20,10))
plt.plot(x, wavelength.iloc[:,2:6], '.')
plt.plot(length_interp_bin, wavelength_interp.T, 'x')
plt.show()

# 分别读取各波段范围内的插值数据
b1 = wavelength_interp[0, np.argwhere(length_interp_bin==450)[0][0]: np.argwhere(length_interp_bin==520)[0][0]+1]
b2 = wavelength_interp[1, np.argwhere(length_interp_bin==520)[0][0]: np.argwhere(length_interp_bin==590)[0][0]+1]
b3 = wavelength_interp[2, np.argwhere(length_interp_bin==630)[0][0]: np.argwhere(length_interp_bin==690)[0][0]+1]
b4 = wavelength_interp[3, np.argwhere(length_interp_bin==770)[0][0]: np.argwhere(length_interp_bin==890)[0][0]+1]

1.3 使用scipy插值

1.3.1 代码

接上面

# 导入包
from scipy import interpolate

plt.figure(figsize=(15,8))
for kind in ["nearest", "zero", "slinear", "quadratic", "cubic"]:
    # "nearest", "zero"为阶梯插值
    # slinear 线性插值
    # "quadratic", "cubic" 为2阶、3阶B样条插值
    for i in range(2,6):
        f = interpolate.interp1d(x,  wavelength.iloc[:,i], kind=kind)
        inter_wavelength = f(length_interp_bin)
        plt.plot(length_interp_bin, inter_wavelength,'x', label=str(kind+'_band'+str(i)))
plt.legend()
plt.show()

end