标准化降水SPI计算

HFUT1202

已于 2023-08-09 08:54:48 修改

阅读量614

点赞数 2

文章标签：算法

于 2023-08-09 08:53:55 首次发布

本文链接：https://blog.csdn.net/weixin_46982703/article/details/132180208

版权

clear;clc;
scale = 12;
Rm = xlsread('mon.xlsx');
Rspi = [-0.5 -1 -1.5 -2];

%%%%%%%%%%%%scale12=12;scale12=12;%%%%%%%%%
SPI=length(Rm);
L = scale-1;
R=zeros(SPI-L,1);
for j=1:(SPI-L)
    R(j)=sum(Rm(j:j+L));
end
Rn=R(1:12:length(R),:);
[x]=find(Rn==0);
Rn_nozero=Rn;
Rn_nozero(x)=[];
q=length(x)/length(Rn);
parm=gamfit(Rn_nozero);
G=q+(1-q)*gamcdf(Rn,parm(1),parm(2));
SPIn12=norminv(G,0,1);
SPI12 = SPIn12;



%%%%%%%%%%%%scale12=12;scale12=12;%%%%%%%%%

%%%%%%%%%%%%scale1 = 1;scale1 = 1;%%%%%%%%%
Dm=Rm;
for i1=1:12                  %月份1、2.。。12
    k1=i1:12:length(Dm);     %60到19年的 1、2.。12月
    Dn=Dm(k1);               %对应所有年某月的dm
    [x]=find(Dn==0);
    Rn_nozero=Dn;
    Rn_nozero(x)=[];
    q=length(x)/length(Dn);
    parm=gamfit(Rn_nozero);
    G=q+(1-q)*gamcdf(Dn,parm(1),parm(2));
    SPI1(:,i1)=norminv(G,0,1);
end
SPI1 = SPI1';
SPI1 = SPI1(:);
%%%%%%%%%%%%scale1 = 1;scale1 = 1;%%%%%%%%%



%%%%%%%%%%%%scale3 = 3;scale3 = 3;%%%%%%%%%
scale3= 3;
L3 = scale3-1;
R3=zeros(SPI-L3,1);
for j=1:(SPI-L3)
    R3(j)=sum(Rm(j:j+L3));
end
%%%%%%%%%%%%scale3 = 3;scale3 = 3;%%%%%%%%%

%%%%%%%%%%%%春春春春春春春春春春春春春;%%%%%%%%%
Rspr = R3(3:12:length(R3));
[x]=find(Rspr==0);
Rn_nozero=Rspr;
Rn_nozero(x)=[];
q=length(x)/length(Rspr);
parm=gamfit(Rn_nozero);
G=q+(1-q)*gamcdf(Rspr,parm(1),parm(2));
SPInspr=norminv(G,0,1);
SPIspr = SPInspr;
%%%%%%%%%%%%春春春春春春春春春春春春春;%%%%%%%%%

%%%%%%%%%%%%夏夏夏夏夏夏夏夏夏夏夏夏夏夏夏;%%%%%%%%%
Rsum = R3(6:12:length(R3));
[x]=find(Rsum==0);
Rn_nozero=Rsum;
Rn_nozero(x)=[];
q=length(x)/length(Rsum);
parm=gamfit(Rn_nozero);
G=q+(1-q)*gamcdf(Rsum,parm(1),parm(2));
SPInsum=norminv(G,0,1);
SPIsum = SPInsum;
%%%%%%%%%%%%夏夏夏夏夏夏夏夏夏夏夏夏夏夏夏;%%%%%%%%%

%%%%%%%%%%%%秋秋秋秋秋秋秋秋秋秋秋秋;%%%%%%%%%
Raut = R3(9:12:length(R3));
[x]=find(Raut==0);
Rn_nozero=Raut;
Rn_nozero(x)=[];
q=length(x)/length(Raut);
parm=gamfit(Rn_nozero);
G=q+(1-q)*gamcdf(Raut,parm(1),parm(2));
SPInaut=norminv(G,0,1);
SPIaut= SPInaut;
%%%%%%%%%%%%秋秋秋秋秋秋秋秋秋秋秋秋;%%%%%%%%%

%%%%%%%%%%%%冬冬冬冬冬冬冬冬冬冬冬冬冬冬;%%%%%%%%%
Rwin = R3(12:12:length(R3));
[x]=find(Rwin==0);
Rn_nozero=Rwin;
Rn_nozero(x)=[];
q=length(x)/length(Rwin);
parm=gamfit(Rn_nozero);
G=q+(1-q)*gamcdf(Rwin,parm(1),parm(2));
SPInwin=norminv(G,0,1);
SPIwin = SPInwin;
SPIwin(length(SPIaut),1) = nan;
%%%%%%%%%%%%冬冬冬冬冬冬冬冬冬冬冬冬冬冬;%%%%%%%%%

SP = [SPI12 SPIspr SPIsum SPIaut SPIwin];
SP0 = nan(length(SPI1)-length(SPI12),5);
SP = [SP;SP0];
SP = [SPI1 SP];

[m,n] = size(SP);
for i=1:n
    han(1,i) = length(find(SP(:,i)>Rspi(1)));
    han(2,i) = length(find(SP(:,i)>Rspi(2)&SP(:,i)<=Rspi(1)));
    han(3,i) = length(find(SP(:,i)>Rspi(3)&SP(:,i)<=Rspi(2)));
    han(4,i) = length(find(SP(:,i)>Rspi(4)&SP(:,i)<=Rspi(3)));
    han(5,i) = length(find(SP(:,i)<=Rspi(4)));
end
han1 = han ./sum(han);
[m,n] = size(han1);
han1(m+1:length(SPI1),:)=nan;
SP = [han1 SP];
xlswrite('SPEI-SPI.xlsx',SP,'SPI年春夏秋冬');

月降水数据，计算SPI1 SPI3 SPI12