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()