日数据 计算SPEI

最新推荐文章于 2024-03-01 10:23:38 发布
最新推荐文章于 2024-03-01 10:23:38 发布
阅读量512 收藏 7
点赞数 1
文章标签： 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_46982703/article/details/132183900
版权

 

 

clc;clear;
data = xlsread('成山头1958-2020.xlsx');
data(any((data(:,3)*100+data(:,4)==229),2),:)=[];
%%%计算PET= E0  注意表格数据气压单位百帕 风俗 两米风速
ELE = data(1,14); LAT = data(1,13);as = 0.25;bs = 0.5;
M = data(:,3);
Ri = data(:,4); 
JJ = [M Ri];
for i = 1:length(JJ)
    if JJ(i,1)==1
        J(i,1)=JJ(i,2);
    elseif JJ(i,1)==2
        J(i,1)=JJ(i,2)+31;
    elseif JJ(i,1)==3
        J(i,1)=JJ(i,2)+31+28;
    elseif JJ(i,1)==4
        J(i,1)=JJ(i,2)+31+28+31;
    elseif JJ(i,1)==5
        J(i,1)=JJ(i,2)+31+28+31+30;
    elseif JJ(i,1)==6
        J(i,1)=JJ(i,2)+31+28+31+30+31;
    elseif JJ(i,1)==7
        J(i,1)=JJ(i,2)+31+28+31+30+31+30;
    elseif JJ(i,1)==8
        J(i,1)=JJ(i,2)+31+28+31+30+31+30+31;
    elseif JJ(i,1)==9
        J(i,1)=JJ(i,2)+31+28+31+30+31+30+31+31;
    elseif JJ(i,1)==10
        J(i,1)=JJ(i,2)+31+28+31+30+31+30+31+31+30;
    elseif JJ(i,1)==11
        J(i,1)=JJ(i,2)+31+28+31+30+31+30+31+31+30+31;
    else J(i,1)=JJ(i,2)+31+28+31+30+31+30+31+31+30+31+30;
    end
end

R = data(:,5);            %降水量
Pre = data(:,6);            %气压
Tmean = data(:,7);        %平均温度
Ws = data(:,10);           %风速
Hum = data(:,11);         %湿度
Sun = data(:,12);         %日照
u2=Ws;%u2为2米高处的风速
%∆为饱和水汽压曲线的斜率
delta=4098*0.6108*exp(17.27*Tmean./(Tmean+237.3))./(Tmean+237.3).^2;
eS=0.6108*exp(17.27*Tmean./(Tmean+237.3));  %饱和水汽压
eA=eS.*Hum/100;                    %ea为实际水汽压
gamma1 = 0.665/10^3 *Pre/10;%r为湿度计常数，kPa/℃
Phi=LAT*pi/180;     %纬度 由度数化成弧度
%黄赤交角，是指地球公转轨道面（黄道面）与赤道面（天赤道面）的交角
Delt=0.408*sin(2*pi*J/365-1.39);
ws=acos(-tan(Phi)*tan(Delt));  %日落角
dr=1+0.033*cos(2*pi*J/365); %地 日 距离
Gsc=0.0820;  %太阳常数
%Ra为用相对应的蒸发单位表达的地球辐射值，mm/d
Ra=((24*60)/pi)*Gsc*dr.*(ws.*sin(Phi).*sin(Delt)+cos(Phi).*cos(Delt).*sin(ws));
Rso=(as+bs+2*ELE/100000).*Ra;%晴空辐射
N=24*ws/pi;                  %可能日照时数
Rs=(as+bs*(Sun./N)).*Ra;
alpha=0.23;          %作物反射率
Rns=(1-alpha)*Rs;    %短波辐射
%计算净长波辐射RnL
Tkmax=Tmean+273.16;Tkmin=Tmean+273.16;
fcd=1.35*Rs./Rso-0.35;
Rnl=(4.903*10^(-9))*fcd.*(0.34-0.14*sqrt(eA)).*((Tkmax.^4+Tkmin.^4)/2);
Rn=Rns-Rnl;
E0=(0.408 .* delta .*Rn+gamma1 *900 ./(Tmean+273) .*u2 .*(eS-eA))./(delta+gamma1 .*(1+0.34*u2));
PET1 = [data(:,2:4),E0 data(:,5)];
% PET = zeros(ceil(length(data)/365*12),31);
a=1;b=1;
for i = 1:length(PET1)-1
    if PET1(i,2)==PET1(i+1,2)
        PET(a,b)=PET1(i,4);P(a,b)=PET1(i,5);
        PET(a,b+1)=PET1(i+1,4);P(a,b+1)=PET1(i+1,5);
        b=b+1;
    else 
        a = a+1;
        b=1;
    end
end

my_E0=sum(PET,2);
Rm = sum(P,2);
%%%%%%%%%%%%%%%%%ET0  ET0%%%%%%%%%%%%%%%%%%%
D_month = Rm-my_E0;
%%%%%%%%%%%%scale12 = 12;scale12 = 12;%%%%%%%%%
scale12 = 12;
L=length(D_month);
D_monthly=zeros(L-(scale12-1),1);
for ii=1:(L-(scale12-1)) %数据累加成12个月
    D_monthly(ii)=sum(D_month(ii:ii+scale12-1));
end
Dm=D_monthly;
Dn=Dm(1:scale12:length(Dm));

Row=find(Dn); %序号
D2=[Row,Dn]; %序号+data
n1=length(D2); %个数
D3=sortrows(D2,2);%data排序
D4=D3(:,2);
for n2=1:n1
    F(n2,1)=(n2-0.35)/n1;
end
w0=mean(D4);
w1=mean((1-F).*D4);
w2=mean((1-F).^2.*D4);
a=(2*w1-w0)/(6*w1-w0-6*w2); %白塔
b=a*(w0-2*w1);
c=gamma(1+1/a)*gamma(1-1/a);
d=b/c;                      %阿尔法
r=w0-b;                     %r
for n3=1:n1
    P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
    if P>0.5
        t(n3)=sqrt(-2*log(1-P(n3)));

        SPEIm(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
    else
        t(n3)=sqrt(-2*log(P(n3)));

        SPEIm(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
    end
end
SPEI0(:,1)=D3(:,1);
SPEI0(:,2)=SPEIm;
mySPEI=sortrows(SPEI0,1);
SPEI12=mySPEI(:,2);
%%%%%%%%%%%%scale12 = 12;scale12 = 12;%%%%%%%%%

%%%%%%%%%%%%scale1 = 1;scale1 = 1;%%%%%%%%%
Dm=D_month;
for i1=1:12                  %月份1、2.。。12
    k1=i1:12:length(Dm);     %60到19年的 1、2.。12月
    Dn=Dm(k1);               %对应所有年某月的dm
    Row=find(Dn);            %序号
    D2=[Row,Dn];             %序号＋dm
    n1=length(D2);           %个数
    D3=sortrows(D2,2);       %按dm排序
    D4=D3(:,2);              %取排序都的dm
    for n2=1:n1
        F(n2,1)=(n2-0.35)/n1;
    end
    w0=mean(D4);
    w1=mean((1-F).*D4);
    w2=mean((1-F).^2.*D4);
    a=(2*w1-w0)/(6*w1-w0-6*w2);
    b=a*(w0-2*w1);
    c=gamma(1+1/a)*gamma(1-1/a);
    d=b/c;
    r=w0-b;
    for n3=1:n1
        P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
        if P>0.5
            t(n3)=sqrt(-2*log(1-P(n3)));

            SPEIm(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
        else
            t(n3)=sqrt(-2*log(P(n3)));

            SPEIm(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
        end
        SPEI0(:,1)=D3(:,1);
        SPEI0(:,2)=SPEIm;
        mySPEI=sortrows(SPEI0,1);
        SPEI=mySPEI(:,2);
    end
    SPEIn1(:,i1)=SPEI;
end
SPEIn1 = SPEIn1';
SPEI1 = SPEIn1(:);
%%%%%%%%%%%%scale1 = 1;scale1 = 1;%%%%%%%%%

%%%%%%%%%%%%scale3 = 3;scale3 = 3;%%%%%%%%%
scale3 = 3;
L=length(D_month);
D_monthly3=zeros(L-(scale3-1),1);
for ii=1:(L-(scale3-1)) %数据累加成12个月
    D_monthly3(ii)=sum(D_month(ii:ii+scale3-1));
end
Dm3=D_monthly3;
Dspr = Dm3(3:12:length(Dm3));
Dsum = Dm3(6:12:length(Dm3));
Daut = Dm3(9:12:length(Dm3));
Dwin = Dm3(12:12:length(Dm3));
%%%%%%%%%%%%scale3 = 3;scale3 = 3;%%%%%%%%%

%%%%%%%%%%%%春春春春春春春春春春春春春;%%%%%%%%%
Row=find(Dspr); %序号
D2=[Row,Dspr]; %序号+data
n1=length(D2); %个数
D3=sortrows(D2,2);%data排序
D4=D3(:,2);
for n2=1:n1
    F(n2,1)=(n2-0.35)/n1;
end
w0=mean(D4);
w1=mean((1-F).*D4);
w2=mean((1-F).^2.*D4);
a=(2*w1-w0)/(6*w1-w0-6*w2); %白塔
b=a*(w0-2*w1);
c=gamma(1+1/a)*gamma(1-1/a);
d=b/c;                      %阿尔法
r=w0-b;                     %r
for n3=1:n1
    P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
    if P>0.5
        t(n3)=sqrt(-2*log(1-P(n3)));

        SPEIm(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
    else
        t(n3)=sqrt(-2*log(P(n3)));

        SPEIm(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
    end
end
SPEI0(:,1)=D3(:,1);
SPEI0(:,2)=SPEIm;
mySPEI=sortrows(SPEI0,1);
SPEIspr=mySPEI(:,2);
%%%%%%%%%%%%春春春春春春春春春春春春春;%%%%%%%%%

%%%%%%%%%%%%夏夏夏夏夏夏夏夏夏夏夏夏夏夏夏;%%%%%%%%%
Row=find(Dsum); %序号
D2=[Row,Dsum]; %序号+data
n1=length(D2); %个数
D3=sortrows(D2,2);%data排序
D4=D3(:,2);
for n2=1:n1
    F(n2,1)=(n2-0.35)/n1;
end
w0=mean(D4);
w1=mean((1-F).*D4);
w2=mean((1-F).^2.*D4);
a=(2*w1-w0)/(6*w1-w0-6*w2); %白塔
b=a*(w0-2*w1);
c=gamma(1+1/a)*gamma(1-1/a);
d=b/c;                      %阿尔法
r=w0-b;                     %r
for n3=1:n1
    P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
    if P>0.5
        t(n3)=sqrt(-2*log(1-P(n3)));

        SPEIm(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
    else
        t(n3)=sqrt(-2*log(P(n3)));

        SPEIm(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
    end
end
SPEI0(:,1)=D3(:,1);
SPEI0(:,2)=SPEIm;
mySPEI=sortrows(SPEI0,1);
SPEIsum=mySPEI(:,2);
%%%%%%%%%%%%夏夏夏夏夏夏夏夏夏夏夏夏夏夏夏;%%%%%%%%%

%%%%%%%%%%%%秋秋秋秋秋秋秋秋秋秋秋秋;%%%%%%%%%
Row=find(Daut); %序号
D2=[Row,Daut]; %序号+data
n1=length(D2); %个数
D3=sortrows(D2,2);%data排序
D4=D3(:,2);
for n2=1:n1
    F(n2,1)=(n2-0.35)/n1;
end
w0=mean(D4);
w1=mean((1-F).*D4);
w2=mean((1-F).^2.*D4);
a=(2*w1-w0)/(6*w1-w0-6*w2); %白塔
b=a*(w0-2*w1);
c=gamma(1+1/a)*gamma(1-1/a);
d=b/c;                      %阿尔法
r=w0-b;                     %r
for n3=1:n1
    P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
    if P>0.5
        t(n3)=sqrt(-2*log(1-P(n3)));

        SPEIm(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
    else
        t(n3)=sqrt(-2*log(P(n3)));

        SPEIm(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
    end
end
SPEI0(:,1)=D3(:,1);
SPEI0(:,2)=SPEIm;
mySPEI=sortrows(SPEI0,1);
SPEIaut=mySPEI(:,2);
%%%%%%%%%%%%秋秋秋秋秋秋秋秋秋秋秋秋;%%%%%%%%%

%%%%%%%%%%%%冬冬冬冬冬冬冬冬冬冬冬冬冬冬;%%%%%%%%%
Row=find(Dwin); %序号
D2=[Row,Dwin]; %序号+data
n1=length(D2); %个数
SPEI0 = zeros(n1,2);
D3=sortrows(D2,2);%data排序
D4=D3(:,2);
for n2=1:n1
    Fw(n2,1)=(n2-0.35)/n1;
end
w0=mean(D4);
w1=mean((1-Fw).*D4);
w2=mean((1-Fw).^2.*D4);
a=(2*w1-w0)/(6*w1-w0-6*w2); %白塔
b=a*(w0-2*w1);
c=gamma(1+1/a)*gamma(1-1/a);
d=b/c;                      %阿尔法
r=w0-b;                     %r
for n3=1:n1
    P(n3)=1-1/(1+(d/(D4(n3)-r))^a);
    if P>0.5
        t(n3)=sqrt(-2*log(1-P(n3)));

        SPEImw(n3)=-(t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3));
    else
        t(n3)=sqrt(-2*log(P(n3)));

        SPEImw(n3)=t(n3)-(2.515517+0.802853*t(n3)+0.010328*t(n3).^2)/(1+1.432788*t(n3)+0.189269*t(n3).^2+0.001308*t(n3).^3);
    end
end
SPEI0(:,1)=D3(:,1);
SPEI0(:,2)=SPEImw;
mySPEI=sortrows(SPEI0,1);
SPEIwin=mySPEI(:,2);
SPEIwin(length(SPEIaut))=nan;

%%%%%%%%%%%%冬冬冬冬冬冬冬冬冬冬冬冬冬冬;%%%%%%%%%

SP = [SPEI12 SPEIspr SPEIsum SPEIaut SPEIwin];
SP0 = nan(length(SPEI1)-length(SPEI12),5);
SP = [SP;SP0];
SP = [SPEI1 SP];
xlswrite('SPEI-SPI.xlsx',SP,'SPEI年春夏秋冬');

 

HFUT1202
关注
r语言的runmed函数_干旱指数新宠—SPEI的C语言版、R语言版与单独执行软件SPEI Calculator...
weixin_30825953的博客
02-11 1458
Drought indices(SPEI)Drought indices()所属R语言包：SPEICalculation of the Standardized Precipitation-Evapotranspiration Index (SPEI) and the Standardized Precipitation Index (SPI).计算的标准化降水的蒸散指数(SPEI)和标准化降水指...
Python：SPEI指数的计算
houhaha_的博客
08-09 1563
将nc文件修改为想要格式：参考sample_precip.nc和sample_tavg.nc(1）根据ERA5小时数据求月总降水和月平均气温在python中(2) 在nco中修改nc文件的格式ncatted -a units,precip,c,c,millimeter in.nc #给precip变量添加单位属性为“millimeter”ncatted -a units,tavg,c,c,K in.nc #给tavg变量添加单位属性为“K”
R计算SPEI指数
snackpdd的博客
04-18 4492
library(SPEI) # ===========读入数据 # 降水数据 prec = read.table('D:\\station\\prec.txt',header = T) # 温度数据 temp = read.table('D:\\station\\temp.txt',header = T) # 每一个站点对应的纬度数据 lat = read.table('D:\\station\\lat.txt',header = T) # 输出路径 print(ncol(temp)) print(
基于Python的Climate Indices库计算SPEI(标准化降水蒸散发指数)04—不同站点不同时间尺度的SPEI计算
EWBA_GIS_RS_ER的博客
12-21 1221
不同站点不同时间尺度SPEI的计算
SPEI干旱指数和ETo潜在蒸散量计算在python上的实现
04-30 8538
1、计算公式 2、输入数据 3、需求 计算每个站点的spei和ETo以及按年、季节、月统计每个站点的干旱等级、干旱湿润指标统计等 4、结果 代码
日数据 计算SPEI SPEI1 SPEI12
08-09
计算SPEI的过程包括以下几个关键步骤： 1. **数据收集**：首先，需要收集每日的降水量数据，这是SPEI计算的基础。同时，也需要获取相应的潜在蒸发量数据，这通常可以通过气象参数如温度、湿度、风速等估算得到。 2...
spei计算程序附测试数据
05-14
测试数据集通常包括了历史气象观测数据，比如降水量、平均气温、最高气温、最低气温等，这些数据是计算SPEI的必要输入。通过对这些数据的处理，用户可以直观地看到SPEI计算的输出结果，并与已知的干旱事件进行比较，...
用MALTAB计算SPEI指数的代码
01-06
计算标准化降水蒸散指数（Standardized Precipitation Evapotranspiration Index，SPEI）的Matlab代码可以使用一些统计和气象数据处理的工具箱。以下是一个简单的示例代码，其中使用了Matlab的Climate Data Toolbox...
SPEI的计算流程
m0_54800671的博客
03-12 7989
介绍了SPEI的计算流程和PET的三种计算方式
gma 教程 | 气候气象 | 计算标准化降水蒸散指数（SPEI）
weixin_42155937的博客
06-28 1万+
【基于 Excel 降水和蒸散数据计算 SPEI】 【基于 GTiff 栅格降水和蒸散数据计算 SPEI】
r语言计算标准化降水指数
11-05
r语言程序，用于计算标准化降水指数（SPI），用于干旱分级
标准化降水蒸散发指数 SPEI 算法及代码（Java）
热门推荐
qq_40821274的博客
06-10 2万+
SPEI计算及Java算法实现前言关于SPEI 前言 关于SPEI
SPEI指数的原理和计算过程
weixin_44259522的博客
12-22 2777
SPEI（Standardized Precipitation Evapotranspiration Index）的计算过程涉及标准化降水指数（SPI）和蒸散发的影响。在标准化降水指数（SPI）和标准化蒸散发蒸发指数（SPEI）中，"尺度"通常指的是时间尺度，表示考虑降水或蒸散发影响的时间跨度。例如，1月尺度的SPI表示了一个月内的降水情况，而6月尺度的SPI表示了六个月内的降水情况。最终的SPEI值用于表示相对于历史同期的标准化降水和蒸散发的差异，以评估干旱和湿润的情况。是期望的蒸散发量的均值，
gma 教程：计算标准化降水蒸散指数（SPEI）
最新发布
weixin_42155937的博客
03-01 2448
【标准化降水蒸散指数】。基于 Log-Logistic 分布计算标准化降水蒸散指数。PRE:array。降水量（mm）。PET:array。潜在蒸散量（mm）。Axis =int。计算轴。如果不设置（None），多维数据会将所有数据展开到一维计算。Scale =int。时间尺度。默认为 1。例如：1月、3月或其他。int。周期。默认为 12。例如：月数据可以以 12 为周期。str。用于内部拟合/变换计算的分布类型。默认为 LogLogistic。支持的分布类型‘Gamma’：伽马分布；
基于Python的Climate Indices库计算SPEI(标准化降水蒸散发指数)03—单站点不同时间尺度的SPEI计算
EWBA_GIS_RS_ER的博客
09-23 7825
基于Python计算特定站点不同时间尺度的SPEI
matlab利用nc数据计算spei
05-13
下面将介绍如何使用MATLAB利用nc数据计算SPEI。 第一步是读取nc文件中的数据。通过使用“ncinfo”来读取nc文件，确定变量、尺寸和属性等信息。然后使用“ncread”函数将数据读入Matlab工作空间。 接下来，我们需要...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

HFUT1202

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值