clear
fname=‘D:\2.tif‘;%用于提取區(qū)域坐標系的模板tif圖像，以便寫出tif格式
?[dataR]?=?geotiffread（fname）;%提取坐標系
?info?=?geotiffinfo（fname）;
?mydir=‘D:\modis_tif\‘;%路徑，批處理
?temp=dir（[mydir‘*.tif‘]）;%找到路徑下所有TIF文件
?num_temp=length（temp）;
?%循環(huán)處理所有TIF文件
?for?ii=1:num_temp
?????filename=[mydirtemp（ii）.name];%文件名
??????modis_tif=imread（filename）;%讀取tif數(shù)據(jù)
???????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，循環(huán)填補，考慮到窗口大小對邊界的適應(yīng)情況，分各邊界處理
?[mn]=find（isnan（new_modis_tif））;
?for?j=1:length（m）
????m1=m（j）;
????n1=n（j）;
????if?m1<=5&n1<=5%左上角邊界，我這里窗口大小設(shè)置為5*5
?????fill_value=nanmean（new_modis_tif（1:m1+51:n1+5））;?%以該像元為中心5*5窗口內(nèi)可用值均值填補該像元??
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）;
?