标准降水指数（SPI）、标准化降水蒸发蒸腾指数（SPEI）和帕尔默干旱指数（PDSI）的下载、读取和可视化

news2024/10/7 6:46:05

数据下载网址：WWDT Data (dri.edu)https://wrcc.dri.edu/wwdt/data/PRISM/

以SPI为例说明，

标准化降水指数(Standardized Precipitation Index, SPI)是由Mckee et al(.1993)分析美国科罗拉多干旱时,发现降水服从偏态分布,基于此提出了标准化降水指数。标准化降水指数通过计算给定时间内降雨量的累计概率,比较客观地表述了多时间尺度下的降水概率,反映降水因素对干旱的影响。Mckee et al.(1993)认为降水量的减少是导致干旱的主要因素,其他气候因子一般变化较小,对干旱产生的影响不大,因此SPI使用降水量作为唯一参数。SPI能较好表达因降水量的大小反映干旱状况,操作简单,所需资料仅为降水量,同时该干旱指数在各个地区和各个时段都具有良好的计算稳定性,是继PDSI之后另一种被广泛应用的干旱指数。
标准化降水指数介绍_标准化降水指数spipython-CSDN博客

下载数据时间跨度为从2002年4月至2014年6月，共147个月。提前准备好如下形式的控制文件：

读取SPI数据，并绘制加州中央山谷区域的SPI-6、SPI-12和SPI-24指标时间序列图

clearvars -except Time02_14
info6 = ncinfo('E:\RSE\data\PRISM\SPI\6\spi6_2003_1_PRISM.nc');
fid6=fopen('E:\RSE\data\PRISM\SPI\6\Contr_spi6.txt','r');
fid12=fopen('E:\RSE\data\PRISM\SPI\12\Contr_spi12.txt','r');
fid24=fopen('E:\RSE\data\PRISM\SPI\24\Contr_spi24.txt','r');
load E:\RSE\Result1\Time02_14.mat;
num_file= fscanf(fid6,'%d',1); fig=1;
dir_in6=fscanf(fid6,'%s',1);dir_in12=fscanf(fid12,'%s',1);dir_in24=fscanf(fid24,'%s',1);
for i = 1:num_file 
    file_name{i,1} = fscanf(fid6,'%s',1);
    file_name{i,2} = fscanf(fid12,'%s',1);
    file_name{i,3} = fscanf(fid24,'%s',1);
end
for k=1:num_file
    data6(:,:,k)=rot90(ncread(strcat(dir_in6,file_name{k,1}),'data')); 
    data12(:,:,k)=rot90(ncread(strcat(dir_in12,file_name{k,2}),'data'));
    data24(:,:,k)=rot90(ncread(strcat(dir_in24,file_name{k,3}),'data'));
end

lon = 360-125:0.041666:360-66.5;  lat =50:-0.041666:24.13;[LON,LAT]=meshgrid(lon,lat);
mask = 'E:\RSE\data\Basin\studyregion.vec';rows=3932;
spi6_serie=gmt_grid2regionserie(data6,mask,rows,lon,lat);
spi12_serie=gmt_grid2regionserie(data12,mask,rows,lon,lat);
spi24_serie=gmt_grid2regionserie(data24,mask,rows,lon,lat);
% save SPI6_24.mat spi6_serie spi12_serie spi24_serie

figure(fig)
plot(Time02_14,spi6_serie,'color',[0 0 0],'Linewidth',1.5);hold on;
plot(Time02_14,spi12_serie,'color',[0 1 0],'Linewidth',1.5);
plot(Time02_14,spi24_serie,'--','color',[1 0 1],'Linewidth',1.5);
title('','FontName','Time New Roman','Fontsize',30);
yline(0,'LineStyle','-','Color','red','LabelVerticalAlignment','middle','LabelHorizontalAlignment','center','LineWidth',1.5);
xlim([2002 2015]);
xticks([2002 2004 2006 2008 2010 2012 2014]);
xlabel('Time (year)','Fontname','Time New Roman',"FontSize",10);
ylabel('Index','Fontname','Time New Roman',"FontSize",10);
l1=legend('SPI6','SPI12','SPI24','location','northeast');
set(l1,'box','off',"FontSize",10,'Fontname','Time New Roman');
box on;grid on;fig=fig+1;

值得说明滴是，“Time02_14”为以年为单位的时间变量；函数“gmt_grid2regionserie”如下：

function [plot_region]=gmt_grid2regionserie(grid,dir_msk,rows,lon,lat)
fid         = fopen(dir_msk,'r');
num_points  = rows;
        river_name  = fscanf(fid,'%s',1);
        xv          = zeros(num_points,1);
        yv          = zeros(num_points,1);
        for i=1:num_points
            a       = fscanf(fid,'%f %f',[1 2]);
            xv(i)   = a(1);
            yv(i)   = a(2);
        end
        fclose(fid);
lon_min=min(xv);lon_max=max(xv);
lat_min=min(yv);lat_max=max(yv);
ii=size(grid,2);jj=size(grid,1);num_file=size(grid,3);
if ii==1405 && jj==621
    load E:\RSE\result\lg_msk.mat;%%加州区域的掩膜，区域内为1区域外为0
else
        lg_msk = zeros(jj,ii);
        for j=1:jj
            for i=1:ii
                if lon(i)>=lon_min && lon(i)<=lon_max && lat(j)>=lat_min && lat(j)<=lat_max
                    in =inpolygon(lon(i),lat(j),xv,yv);
                    if in==1
                        lg_msk(j,i) = 1;
                    end
                end
            end
        end
end
grid_weight=zeros(jj,ii);
for j=1:jj
    for i=1:ii
        grid_weight(j,i)=cos(lat(j)*pi/180);
    end
end
grid_weight=grid_weight.*lg_msk;
tmmp=zeros(jj,ii);
% grid(logical(isnan(grid)))=0; % set possible NaNs to zeros
for k=1:num_file
    tmmp(:,:)=grid(:,:,k);
    grid_weight_new=grid_weight;
    tmmppp=tmmp.*grid_weight_new;
    tmmppp(logical(isnan(tmmppp)))=0;% in case, there are NaNs
    grid_sum=sum(sum(tmmppp));
    grid_weight_new(logical(isnan(tmmp)))=0;
    weight_sum=sum(sum(grid_weight_new));
    plot_region(k,1)=grid_sum/weight_sum;
end
end

结果：