Python遥感开发之时序数据的线性插值

0 历史博客
1 实现思路
2 代码实现
3 效果展示

前言：在遇到空间数据的时候，尤其是哨兵、Landsat或者MODIS数据会出现局部值的空缺，为了解决这些值的空缺，通常采用插值的方法，本博客使用时序数据，采用线性插值的方法，完成遥感数据空缺值的处理。

0 历史博客

《Python开发之二维数组空缺值的近邻填充》
《Python开发之手动实现一维线性插值》

1 实现思路

1.首先读取时序数据，对读取的数据拼接在一起，然后转换成3D数组（时间、行和列），然后进行线性插值。
2.线性插值：由于输入的是一个三维数组，需要对输入的三维数据进行逐列处理，填充每列中的NaN值。如果一整列都是NaN，则保持不变。否则，对每列的NaN值进行线性插值，使得NaN值由其相邻的非NaN值来估计。最终得到一个插值完成的3D数组。
3.对三维数组的时间列进行遍历，挨个保存每一个时间对应的数据为TIF文件。

2 代码实现

from osgeo import gdal, gdalnumeric, gdalconst
import numpy as np
import osdef save_tif(data, file, output):ds = gdal.Open(file)shape = data.shapedriver = gdal.GetDriverByName("GTiff")dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32)  # 以float类型进行存储dataset.SetGeoTransform(ds.GetGeoTransform())dataset.SetProjection(ds.GetProjection())dataset.GetRasterBand(1).WriteArray(data)def read_tif02(filepath):data = gdalnumeric.LoadFile(filepath)data = data.astype(np.float32)a = data[0][0]data[data == a] = np.nanreturn data# 进行线性插值
def linear_interpolation(data):x = np.arange(data.shape[0])data_interpolated = np.zeros_like(data)for i in range(data.shape[1]):for j in range(data.shape[2]):y = data[:, i, j]if np.all(np.isnan(y)):data_interpolated[:, i, j] = yelse:mask = np.isnan(y)y[mask] = np.interp(x[mask], x[~mask], y[~mask])data_interpolated[:, i, j] = yreturn data_interpolated"""
通过传入文件夹的路径，输出文件的完整路径
通过传入文件夹的路径，传入输出的路径和名字，输出传出文件的完整路径
"""
def get_data_list(file_path, out = "", out_name = ""):list1 = []  # 文件的完整路径if os.path.isdir(file_path):fileList = os.listdir(file_path)if out_name != "" and out != "":for f in fileList:file_name = "".join(list(filter(str.isdigit, f))) + ".tif"out_data = out + "\\" + out_name + file_namelist1.append(out_data)else:for f in fileList:pre_data = file_path + '\\' + f  # 文件的完整路径list1.append(pre_data)return list1if __name__ == '__main__':rsei_file_path = r"E:\AAWORK\work\研究方向\rsei\data\月\in"output_dir = r"E:\AAWORK\work\研究方向\rsei\data\月\out"rsei_file_path_list = get_data_list(rsei_file_path)all_data_list = []# 读取数据for tif_file in rsei_file_path_list:data = read_tif02(tif_file)all_data_list.append(data)# 转换为3D数组 (时间, 行, 列)data_array = np.array(all_data_list)# 进行线性插值data_interpolated = linear_interpolation(data_array)for i in range(data_interpolated.shape[0]):file_name = os.path.basename(rsei_file_path_list[i])file_name = file_name.replace(".tif","_interpolated.tif")out_file = os.path.join(output_dir, file_name)save_tif(data_interpolated[i], rsei_file_path_list[i], out_file)print(f"Saved: {out_file}")