MODIS、Landsat等遥感影像批量空间插值Matlab代码

clear
fname=‘D:\2.tif‘;%用于提取区域坐标系的模板tif图像，以便写出tif格式
 [dataR] = geotiffread（fname）;%提取坐标系
 info = geotiffinfo（fname）;
 mydir=‘D:\modis_tif\‘;%路径，批处理
 temp=dir（[mydir‘*.tif‘]）;%找到路径下所有TIF文件
 num_temp=length（temp）;
 %循环处理所有TIF文件
 for ii=1:num_temp
     filename=[mydirtemp（ii）.name];%文件名
      modis_tif=imread（filename）;%读取tif数据
       modis_tif=double（modis_tif）;
        re_modis_tif=reshape（modis_tif1[]）;%重排列
        inx=find（re_modis_tif==32767）;%我这里要除去的是32767，可以把它换成其他值
        %inx2=find（re_modis_tif==写入所有要替换和插值的值）
        re_modis_tif（inx）=nan;%将所有要插值的值赋为nan
         new_modis_tif=reshape（re_modis_tif10001000）;
         %找到所有的nan，循环填补，考虑到窗口大小对边界的适应情况，分各边界处理
 [mn]=find（isnan（new_modis_tif））;
 for j=1:length（m）
    m1=m（j）;
    n1=n（j）;
    if m1<=5&n1<=5%左上角边界，我这里窗口大小设置为5*5
     fill_value=nanmean（new_modis_tif（1:m1+51:n1+5））; %以该像元为中心5*5窗口内可用值均值填补该像元  
new_modis_tif（m1n1）=nanmean（fill_value）;
    else
if m1<=5&n1>=5&n1<=995%上边界，后面依次是各个边界，不挨个注释了，我的总的图像是1000*1000个像元
    fill_value=nanmean（new_modis_tif（1:m1+5n1-5:n1+5））;
    new_modis_tif（m1n1）=nanmean（fill_value）;
else
    if m1>=5&m1<=995&n1<=5
      fill_value= nanmean（new_modis_tif（m1-5:m1+51:n1+5））;  
 new_modis_tif（m1n1）=nanmean（fill_value）;
    else
        if m1<=5&n1>=95
            fill_value=nanmean（new_modis_tif（1:m1+5n1-5:1000））; 
      new_modis_tif（m1n1）=nanmean（fill_value）;
        else
            if m1>=5&n1>=995&m1<=995
                fill_value=nanmean（new_modis_tif（m1-5:m1+5n1-5:1000））;  
         new_modis_tif（m1n1）=nanmean（fill_value）;  
            else
                if m1<=5&n1<=995&n1>=5
                    fill_value=nanmean（new_modis_tif（1:m1+5n1-5:n1+5））; 
                   new_modis_tif（m1n1）=nanmean（fill_value）; 
                else
                    if m1<=995&n1<=5&m1>=5
                       fill_value= nanmean（new_modis_tif（m1-5:m1+51:n1+5））; 
                     new_modis_tif（m1n1）=nanmean（fill_value）;