鼎城网站建设,上海服饰网站建设,重庆渝中区企业网站建设联系电话,网站建设dw 什么软件Python遥感开发之时序数据的线性插值 0 历史博客1 实现思路2 代码实现3 效果展示 前言#xff1a;在遇到空间数据的时候#xff0c;尤其是哨兵、Landsat或者MODIS数据会出现局部值的空缺#xff0c;为了解决这些值的空缺#xff0c;通常采用插值的方法#xff0c;本博客使… 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))) .tifout_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 rE:\AAWORK\work\研究方向\rsei\data\月\inoutput_dir rE:\AAWORK\work\研究方向\rsei\data\月\outrsei_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(fSaved: {out_file})3 效果展示
没有插值之前的TIF数据 插值之后的TIF数据
