基于MATLAB批量读取栅格影像及处理思路

我个人本是要处理60景遥感影像，需逐像元处理，但是由于刚开始我对matlab掌握不深，所以不会编程，就全部用最低级的方法也就是ArcGIS来处理，two thousand years later…我用arcgis处理完毕。

后来又研究了一番，写出了完整的代码，想想都是泪/(ㄒoㄒ)/~~

以下代码是批量读取数据后，求了一下均值，其他针对多景影像逐像元栅格计算等处理的都是同样的思路。

代码思路：
（1）第一部分：批量读取文件夹下的所有文件，通用；
（2）第二部分（核心）：利用for循环创建三维矩阵。一景栅格影像是二维m×n矩阵，那么对于多景影像的话，就可以将其排列为三维矩阵，新加的影像位于三维矩阵的页上，也就是说有多少景影像，三维矩阵就有多少页。进度条可要可不要。
（3）将处理得到的结果栅格影像输出，也是通用的。

具体代码如下（没有任何删减）：

%%选择目标文件夹,注意文件夹里的数据名要按顺序排列
%特别是进行时间序列分析时
path=uigetdir('','Select a location'); %选择文件夹位置，获取路径
fileList=dir(fullfile(path,'*.tif')); %获取文件夹下所有tif文件
n=length(fileList); %获取文件个数

%%创建三维矩阵，并获取参考矩阵和投影信息
h=waitbar(0,'Running','name','3DMAT'); %进度条1
for i=1:n
    filePath=fullfile(path,fileList(i).name);
    data(:,:,i)=geotiffread(filePath);
    info=geotiffinfo(filePath);
    R=info.RefMatrix; %参考矩阵
    proj=info.GeoTIFFTags.GeoKeyDirectoryTag; %投影信息
    str=['Completed:','',num2str(round(i*100/n)),'%'];
    waitbar(i/n,h,str); %进度条2
end
delete(h); %删除进度条

%%进行所需各项计算
[m,n,p]=size(data); %行、列、页
for i=1:m %行循环
    for j=1:n %列循环
        pixel=reshape(data(i,j,:),[p,1]); %重点！每页对应的像元转为一列
        pixelMean(i,j)=mean(pixel); %求对应像元的均值,并将其放在原始(i,j)位置
    end
end

%%将计算结果写入到新的TIFF影像数据中
[outName,outPath]=uiputfile('*.tif','Save a new tif file');
geotiffwrite(fullfile(outPath,outName),pixelMean,...
    R,'GeoKeyDirectoryTag',proj); %使输出的投影与原始一致

注意：
（1）for循环中info信息的读取必不可少，参考矩阵信息和投影信息是为了最后输出的栅格影像与原始影像保持一致，记住就好。有matlab基础的话，每行代码的意思不难理解，重要的部分注释已标注。
（2）对遥感影像一般都是转成tiff格式进行处理，但要注意，有些软件转换为tiff格式后仍会会带有**.hdr文件**（头文件），这种情况下是无法读取的。建议采用IDL将文件格式批量转为tiff后再处理。

希望本篇内容可以对大家有帮助~