用R语言读取Modis时间序列图像

由于直接用R语言处理Modis影像比较麻烦，我的方法是用R语言的MODIS和Raster Packages读取Modis序列影像，并存储为.txt文件，然后用熟悉的MATLAB处理，这样比较方便。

1、打开安装了MODIS和Raster Packages的RStudio，启动这两个Packages。

2、输入如下代码，读取Modis NDVI图像（我使用的数据集是16天合成的MOD13Q1的植被产品，所有影像用MRT提前处理好，放在同一个文件夹下）。

/*********************************************代码开始了*******************************************/

获取所有文件名
fNames<-preStack(pattern='MOD13Q1.*.tif',path='I:/MODIS/h26v04/Mod13Q1/NDVI/')

获取文件时间
times<-orgTime(fNames)
timeV<-times$inputLayerDates

堆叠所有影像
layers<-stack(fNames)

取局部序列数据，testP为行向量（cell类型）。

注意：R语言是行优先存储的，和MATLAB的列优先方式有所不同，后面必须进行转置。
testP<-getValuesBlock(layers,row=330,nrows=100,col=9340,ncols=100)

存储为TXT文件
write.table(testP,file="E:/NDVI.txt",FALSE,quote=FALSE,sep=" ");

/*********************************************代码结束了*******************************************/

3、打开MATLAB，读入生成好的NDVI.txt数据，并自动生成如下读取代码，存储程序为importfile.m

这个代码是自动生成的，没什么好介绍的，最后生成两个变量：

data：10000*70维（70维是因为我们用了三年的影像，共70景，10000为文件大小为100*100，图像按行优先存储）

textdata：存储每幅影像的文件名，是一个1×1的cell类型，没有根据空格把文件名分开

/*********************************************importfile.m*******************************************/

function importfile(fileToRead1)
%IMPORTFILE(FILETOREAD1)
%  Imports data from the specified file
%  FILETOREAD1:  file to read
%  Auto-generated by MATLAB on 25-Jul-2014 11:41:38

% Import the file
newData1 = importdata(fileToRead1);

% Create new variables in the base workspace from those fields.
vars = fieldnames(newData1);
for i = 1:length(vars)
    assignin('base', vars{i}, newData1.(vars{i}));
end

/*********************************************代码结束了*******************************************/

4、为了能够以比较舒服的方式处理读取Modis数据，我们要将读入的数据转换为MATLAB处理起来比较方便的形式。

建立如下m文件：LoadModis.m

/*********************************************LoadModis.m*******************************************/

importfile('E:\NDVI.txt');  %自动读取数据

%将图像名分开存储
str=cell2mat(textdata);
 str=mat2str(str);
 FileNames=regexp(str,'\s+','split');   %以空格为分隔符读取各图像名
 
 %从data中将每幅图像恢复为2维存储，并进行行列转置
 data=data(:,2:71);
 ModisData=reshape(data,100,100,70); %从一维存储中恢复二维图像
 for i=1:70
     ModisData(:,:,i)=ModisData(:,:,i)';    %每一维转置
 end
 
 %显示第一幅图像
 img1=ModisData(:,:,1);
 mi=min(img1(:));
 ma=max(img1(:));
 imshow((img1-mi)./(ma-mi));

/*********************************************代码结束了*******************************************/

为了试验代码的正确性，我们显示了时间序列中的第一幅图像，并将MATLAB显示结果和ERDAS截图的结果进行了对比，证明了本程序的正确性

上：MATLAB截图     下：ERDAS aoi的显示结果