使用Matlab对MODIS数据进行SG滤波

需求

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

补充说明

1. 本例简单粗暴直接对原数据进行了全部滤波，实际中如果是MODIS数据可以考虑根据其质量控制数据，只将其不可信的像元进行滤波处理，可信像元保持原值。如考虑质量控制可参考其他大神的分享：
   基于matlab的MOD13A2-NDVI的植被指数重建-SG滤波与质量控制文件
2. 批量设NAN值，可方便在ArcGIS中可视化展示
   基于Python的SetNull批量处理