利用matlab进行20摄氏度中止旬计算

clc
clear

%读取tif文件
%[data_tiff, R] = geotiffread('D:\DATA\1979_avg\19790101_avg.tif');
%展示读取结果
%imshow(data_tiff);

filepath='D:\Test\19790727_avg.tif';                                        %%图像名称与路径
Info=imfinfo(filepath); 
% tif='tif';
% format=Info.Format;
% if  (strcmp(format ,tif)==0)
%     disp('载入的不是tif图像，请确认载入的数据');                  %%确保载入的图像是tiff图像
% end
%  
Slice=size(Info,1);                                            %%获取图片z向帧数
% Width=Info.Width;
% Height=Info.Height;
% Image=zeros(Height,Width,Slice);
 
% for i=1:Slice
%     Image(:,:,i)=imread(filepath,i);                                  %%一层一层的读入图像
% end
% Image
%%获取图片信息并判断是否为tifInfo=imfinfo(data_tiff);
%批量读取文件夹内的tif，并分别存贮
Path = 'D:\Test\';                   % 设置数据存放的文件夹路径
File = dir(fullfile(Path,'*.tif'));  % 显示文件夹下所有符合后缀名为.tif文件的完整信息
FileNames = {File.name}';            % 提取符合后缀名为.tif的所有文件的文件名，转换为n行1列

%获取lat lon 最终利用arcgis general获取
% file_name = 'D:\Test\19790727_avg.tif';
% info = geotiffinfo(file_name);
% [x,y] = pix2map(info.RefMatrix, 1, 1);
% [lat,lon] = projinv(info, x,y);

Length_Names = size(FileNames,1);    % 获取所提取数据文件的个数
%批量读取
data_tiff=cell(1,Length_Names);   %创建存储的空cell
for k = 1 : Length_Names
    K_Trace = strcat(Path, FileNames(k));              % 连接路径和文件名得到完整的文件路径
%     data_tiff=cell(1,Length_Names); 不能在for循环内，否则会一直重置
    [data_tiff{k}, R] = geotiffread(K_Trace{1,1});     %读取二维矩阵，由于data_tiff是1*n的元胞数组，可以用单个k表示第k个内容
%                                                        由于K_Trace是元胞数组格式，需要加{1,1}才能得到字符串
    info = geotiffinfo(K_Trace{1,1});
   
%     读取数据（因为这里是.txt格式数据，所以直接用load()函数)
%     eval(['Data',num2str(k),'=','load(K_Trace{1,1})',';']);
%     注意1：eval()函数是括号内的内容按照命令行执行，
%           即eval(['a','=''2','+','3',';'])实质为a = 2 + 3;
%     注意2：由于K_Trace是元胞数组格式，需要加{1,1}才能得到字符串

end

%创建一个用于标记的400*700的空矩阵
mark=ones(400,700)*0;

%当读到小于20时，满足未被标记的情况，且data_tiff并不是NaN的情况下，将该次读到的时间赋值到mark矩阵中
for k = 1 : Length_Names
    A= data_tiff{k} >= 20;
   for i = 1 : 400
       for j = 1 : 700
           if A(i,j) == 0 && mark(i,j) == 0 && data_tiff{k}(i,j) > -50
               mark(i,j) = k;
           end
       end
   end
end

%将mark矩阵中等于0的值赋值为ArcGIS中的NaN,发现不互通，需要在ArcGIA中操作，使用reclassify进行重新复制
% for i = 1 : 400
%        for j = 1 : 700
%            if mark(i,j) == 0
%                mark(i,j)=-999;
%            end
%        end
% end
%     

for i=1:Slice
    J=uint8(mark(:,:,i));           %%一层一层写出图像 
 
    %%imwrite(J,[dir,'.tif'],'WriteMode','Append');
    imwrite(J,[num2str(i,'%04d'),'.tif']);
end
%会发生倒置现象，需要交换一下
data=flipud(J); 
%输出为tif文件
R = georasterref('RasterSize', size(data),'Latlim', [double(15.0007629395) double(55.0007629395)] , 'Lonlim',[double(70.0000038147) double(140.000003815)]);
   geotiffwrite(['D:\1979a.tif'],data,R);
    disp(['done'])