【MATLAB项目实战】基于SPI指数的某地区地区干旱时空特征分析

数据介绍
2001年-2020年降水数据
数据名称:中国区域1km分辨率逐月降水量数据集（2000-2020年）
数据来源：http://www.geodata.cn/data/datadetails.html?dataguid=2329433&docid=2688
时间分辨率：1 month
原数据时间范围：2000-01-2020-12
拟用数据时间范围：2001-01-2020-12
空间分辨率：1km
空间范围：china
单位:mm
缩放因子：暂无

IPCC第4 次评估报告指出，1906——2005年全球平均气温升高了0.74 ℃。在这个背景下,全球范围内特大干旱、高温等极端天气气候事件发生频率和强度呈增加趋势;而不断变化的气候又可能导致极端天气气候事件在发生时间、频率、强度、空间范围和持续时间上产生变化。据世界气象组织的统计数据﹐气象灾害约占自然灾害的70%,而干旱灾害又占气象灾害的50%。作为全球最为常见的自然灾害,干旱具有发生频率高、持续时间长、影响范围广、对农业生产威胁大、对生态环境和社会经济影响深远的特点。
中国大部分地区位于东亚季风的两类子系统——“东亚热带季风(南海季风)”和“东亚副热带季风”共同影响的地区,是干旱严重的国家。近年来,中国干旱事件呈现增长的趋势。具体来说，1951——1990年40年间有8年出现重大干旱事件﹐1991——2000 年10年间有5年出现重大干旱事件﹐2001—2011年11年间有7年出现重大干旱事件。民政部的统计资料显示,1999—2005年,中国干旱灾害发生区域不断扩大,因旱灾造成作物绝收面积急剧增加[3]。位于干旱半干旱区的中国西北部是全球气候变化响应最敏感的地带。随着气候变暖﹐该地区的土地荒漠化现象趋于严重,气象灾害增多,生态环境恶化[6。大量研究表明该地区部分出现暖湿迹象,但总体干旱化趋势已经在中、小时空尺度上得到识别。为了研究干旱的发生和发展规律,科学家们对干旱指标做了大量的研究。总的来说分为单因子指标和多因子指标两类。单因子指标突出降水量这一主要影响因子来反映干旱变化,比较常见的单因子指标有降水距平百分率、无雨日数、Z指数和标准化降水指数(Standardized precipitation index,SPI)等,这类指标资料获取比较容易,计算简单,而且由于指标不涉及具体的干旱机理,所以具有较强的时空适应性;而多因子指标多从水分平衡角度出发,主要强调干旱形成的机理和过程,包括Palmer干旱指标﹑综合气象干旱指数和相对湿润度指数等,这类指标的物理机制较为明确,计算复杂,对资料要求高,有些参数需要实验来确定,还有一些参数不能依靠实验获得,只有通过经验获取，从而大大降低了计算精度。国内外很多学者对不同指标监测干旱事件的效果做了比较,发现基于降水资料的SPI指数计算简单,稳定性好,具有多时间尺度和时空可比性等优势,在国内外的干旱监测中得到广泛应用。
本文采用SPI指数对某地区地区干旱情况进行时空分析

SPI指数是McKee等在评估美国科罗拉多州干旱状况时提出的。对于某一地区而言,确定时间段内的降水量一般呈规律性的波动变化。基于这个事实,认为某一地区确定时间段内的降水量若比该时间段内多年平均降水量(可以认为是正常降水量)偏少到一定程度,就认为该地区该时段发生了干旱;反之﹐若偏多到一定程度﹐则发生洪涝。假定降水量的变化服从gamma分布,运用数学方法将降水量的累计频率分布转化为标准正态分布﹐最后求得SPI。SPI指数具有无量纲﹑标准化的特点﹐能够比较不同地区、不同时段的降水水平,可较好地反映干旱强度和持续时间。SPI指数具有多时间尺度的优势,不仅可以反映短时间尺度内降水量的变化，也可以反映长时间尺度的水资源演变情况。本文计算了1个月、3个月、6个月和12个月的SPI值﹐分别用SPI,SPI3.SPI6和 SPI12代表。

具体原理如下：

根据SPI指数的正态分布曲线,参照国家气候中心的划分标准划分旱涝等级，并确定相应的SPI指数界限值(表1)。

实验结果如下

这个地区不太会干旱呢 都还挺湿润呢 哈哈

整理出SPI趋势图如下，选择的时间尺度为3个月：

导入ArcMap可更改颜色

代码：

clear all;
[a,R]=geotiffread('C:\Users\yuanyuan\Desktop\栅格数据SPI计算\某地区降水数据\2001_01.tif'); %先导入投影信息
info=geotiffinfo('C:\Users\yuanyuan\Desktop\栅格数据SPI计算\某地区降水数据\2001_01.tif');%先导入投影信息
[m,n]=size(a);

datasum=[];   %生成一个像素个数*月数的矩阵
for j=2001:2020
    for i=1:12      %起止月份
        if i<10
        filename=['C:\Users\yuanyuan\Desktop\栅格数据SPI计算\某地区降水数据\',int2str(j),'_0',int2str(i),'.tif'];     %读入文件名   
        data=importdata(filename);  %导入数据
        data=reshape(data,m*n,1);   %reshape 改变矩阵形式为m*n行、1列
        datasum=[datasum  data];          %把每月的数据依次放到datasum的每一列
    
        else 
        filename=['C:\Users\yuanyuan\Desktop\栅格数据SPI计算\某地区降水数据\',int2str(j),'_',int2str(i),'.tif'];     %读入文件名   
        data=importdata(filename);  %导入数据
        data=reshape(data,m*n,1);   %reshape 改变矩阵形式为m*n行、1列
        datasum=[datasum  data];          %把每月的数据依次放到datasum的每一列
    
        end
    end    
end
%% 数据处理：非负
for i=1:size(datasum,1)
    for j=1:size(datasum,2)
        if datasum(i,j)<0
            datasum(i,j)=0;
        end
    end
end
%% SPI计算
% scale 时间尺度
% nseas 季节数量
scale=[1,3,6,12];        
nseas=12;
for k=1:240
% 一个月尺度SPI，其他时间尺度与此类似
    Z(k) = mean(SPI(datasum(:,k),scale(1),nseas)); % 取一下均值
end
figure(2)
subplot(221)
plot(Z)
ylabel('SPI指数')
title('时间尺度为1个月：SPI_1')

for k=1:240
% 3个月尺度SPI，其他时间尺度与此类似
    Z(k) = mean(SPI(datasum(:,k),scale(2),nseas)); % 取一下均值
end
figure(2)
subplot(222)
plot(Z)
ylabel('SPI指数')
title('时间尺度为3个月：SPI_3')

for k=1:240
% 6个月尺度SPI，其他时间尺度与此类似
    Z(k) = mean(SPI(datasum(:,k),scale(3),nseas)); % 取一下均值
end
figure(2)
subplot(223)
plot(Z)
ylabel('SPI指数')
title('时间尺度为6个月：SPI_6')

for k=1:240
% 12个月尺度SPI，其他时间尺度与此类似
    Z(k) = mean(SPI(datasum(:,k),scale(4),nseas)); % 取一下均值
end
figure(2)
subplot(224)
plot(Z)
ylabel('SPI指数')
title('时间尺度为12个月：SPI_12')
%%
for k=1:240
% 3个月尺度SPI，其他时间尺度与此类似
    Z1{k} = SPI(datasum(:,k),scale(1),nseas); % 取一下均值
end
Z2=cell2mat(Z1');
% 求趋势图
Z2M=mean(Z2);
Trend=reshape(Z2M,[m,n]);
geotiffwrite(strcat('结果/','Trend1.tif'),Trend,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);%注意修改路径

代码下载链接：
https://download.csdn.net/download/qq_45047246/86248543