需求
Savitzky-Golay滤波器（简称为S-G滤波器），是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。
因为云覆盖、气溶胶等大气因素会对MODIS数据等造成噪声影响，有必要用滤波等方法来减小或消除这种影响。
以处理植被指数EVI数据为例
%% 读取一个现有文件，获取EVI数据的矩阵行列大小
[a,R] = geotiffread('D:\MODIS\EVI\evi20010101.tif');
info = geotiffinfo('D:\MODIS\EVI\evi20010101.tif');
[m,n] = size(a);
%% 读取EVI并进行重排列
evi_dz = dir('D:\MODIS\EVI\*.tif'); % 获取所有tif文件名
qs = length(evi_dz); % 总期数
evisum = zeros(m*n,1)+NaN; % 构建和EVI同大小的NaN矩阵
k = 1; % 初始化为1
for i = 1:length(evi_dz)
filename = strcat(evi_dz(i).folder,'\',evi_dz(i).name); % 文件路径和文件名
evi = double(importdata(filename)); % 获取evi数据
evi = reshape(evi,m*n,1); % 转为1列，方便处理
evi(evi<-0.2) = NaN; % 缺失值设为nan
evisum(:,k) = evi; % 第一列为第一个时间的evi
k = k+1;
end
%% 因为有部分缺失数据，所以先进行插值
evisum = fillmissing(evisum,'linear',2); % 对每行进行线性插值，可能会产生超出正常范围的值
evisum(evisum<-0.2) = NaN; % 将小于正常范围最小值-0.2的值设为NAN
evisum(evisum>1)= NaN; % 将大于正常范围最大值1的值设为NAN
%% SG滤波
for i = 1:m*n
data1 = evisum(i,:);
data_sg = sgolayfilt(data1,3,5); % sg滤波，参数设置为默认的3,5
evisum(i,:) = data_sg; % 重建EVI
end
%% 输出tif影像
for k = 1:qs
data2 = evisum(:,k);
evi = reshape(data2,m,n); % 重排列为原矩阵大小
evi(evi<-0.2) = NaN; % 超出正常范围的设为NaN
evi(evi>1) = NaN; % 超出正常范围的设为NaN
evi(isnan(evi)) = -9999; % 将NAN值设为-9999
filenem2 = ['D:\MODIS\EVI2\',evi_dz(k).name(4:15)]; % 输出文件夹和文件名
geotiffwrite(filenem2,evi,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag); % 输出
end
补充说明
- 本例简单粗暴直接对原数据进行了全部滤波，实际中如果是MODIS数据可以考虑根据其质量控制数据，只将其不可信的像元进行滤波处理，可信像元保持原值。如考虑质量控制可参考其他大神的分享：
基于matlab的MOD13A2-NDVI的植被指数重建-SG滤波与质量控制文件 - 批量设NAN值，可方便在ArcGIS中可视化展示
基于Python的SetNull批量处理
