Matlab软件求时间序列栅格降水数据的水文年型

[a,R]=geotiffread('C:\Users\WH\Desktop\EPRE\PRE\Idw_PRE_1985.tif');%先导入投影信息
info=geotiffinfo('C:\Users\WH\Desktop\EPRE\PRE\Idw_PRE_1985.tif');
[m,n]=size(a);
cd=2014-1985+1;%时间跨度，根据需要自行修改
datasum=zeros(m*n,cd)+NaN; 
k=1;
for year=1985:2014 %起始年份
    filename=['C:\Users\WH\Desktop\EPRE\PRE\Idw_PRE_',int2str(year),'.tif'];
    data=importdata(filename);
    data=reshape(data,m*n,1);
    datasum(:,k)=data;
    k=k+1;
end
result=zeros(m*n,cd)+NaN;
for i=1:size(datasum,1)
    data=datasum(i,:);
    if min(data)>0 %判断是否是有效值,我这里的有效值必须大于0
        valuesum=[];
        for k1=1:cd
            value=length(find(data>=data(k1)))/(cd+1); %经验频率法求水文年型 
            if value <= 0.5 %这里我想把不同的水文年型赋值
               value1 = 189*0.92; 
            elseif value <= 0.75
                value1 = 316.5*0.92; 
            else 
                value1 = 400.5*0.92; 
            end
            valuesum=[valuesum;value1];
        end
        result(i,:)=valuesum;
    end
end
for k2 = 1:cd
    result2=zeros(m,n)+NaN;
    result2=reshape(result(:,k2),m,n);
    k3=k2+1984;
    filename=['C:\Users\WH\Desktop\EPRE\guanshuidinge\Guanshuidinge_',int2str(k3),'.tif'];
    geotiffwrite(filename,result2,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end