[matlab project practice] analysis of spatial and temporal characteristics of drought in a region based on SPI index

Data is introduced
2001 year -2020 Annual precipitation data
Data name : China region 1km Resolution monthly precipitation data set （2000-2020 year ）
Data sources ：http://www.geodata.cn/data/datadetails.html?dataguid=2329433&docid=2688
Temporal resolution ：1 month
Time range of original data ：2000-01-2020-12
Time range of proposed data ：2001-01-2020-12
Spatial resolution ：1km
The scope of space ：china
Company :mm
Zoom factor ： no

IPCC The first 4 The evaluation report points out ,1906——2005 The global average temperature rose in 0.74 ℃. In this context , Extreme drought in the world 、 The frequency and intensity of extreme weather and climate events such as high temperature are increasing ; And the changing climate may lead to extreme weather and climate events at the time of occurrence 、 frequency 、 Strength 、 Change in spatial scope and duration . According to the statistics of the World Meteorological Organization ﹐ Meteorological disasters account for about of natural disasters 70%, And drought disasters account for 50%. As the most common natural disaster in the world , Drought has a high frequency 、 It lasts a long time 、 Wide range of influence 、 Great threat to agricultural production 、 Characteristics that have a far-reaching impact on the ecological environment and social economy .
Most of China is located in two subsystems of the East Asian monsoon ——“ East Asian tropical monsoon ( The South China Sea monsoon )” and “ East Asian subtropical monsoon ” Jointly affected areas , It is a country with serious drought . In recent years , Drought events in China show an increasing trend . say concretely ,1951——1990 year 40 There was 8 A major drought event occurred in ﹐1991——2000 year 10 There was 5 A major drought event occurred in ﹐2001—2011 year 11 There was 7 A major drought event occurred in . Statistics from the Ministry of Civil Affairs show ,1999—2005 year , China's drought disaster areas are expanding , The area of crop failure caused by drought has increased sharply [3]. Northwest China, located in arid and semi-arid areas, is the most sensitive zone to respond to global climate change . As the climate warms ﹐ The phenomenon of land desertification in this area tends to be serious , Meteorological disasters have increased , Deterioration of ecological environment [6. A large number of studies have shown that some parts of the region show signs of warmth and humidity , But the overall drought trend is already in 、 Recognized on a small space-time scale . In order to study the occurrence and development of drought , Scientists have done a lot of research on drought indicators . Generally speaking, it is divided into single factor indicators and multi factor indicators . The single factor index highlights the precipitation as the main influencing factor to reflect the drought change , A common single factor indicator is the percentage of precipitation anomaly 、 Days without rain 、Z Index and standardized precipitation index (Standardized precipitation index,SPI) etc. , It is easy to obtain such index data , Simple calculation , And because the index does not involve specific drought mechanism , So it has strong space-time adaptability ; The multi factor indicators are mostly from the perspective of water balance , It mainly emphasizes the mechanism and process of drought formation , Include Palmer Drought index ﹑ Comprehensive meteorological drought index and relative humidity index , The physical mechanism of such indicators is relatively clear , The calculation is complicated , High requirements for data , Some parameters need to be determined by experiments , There are also some parameters that cannot be obtained by experiments , Only through experience , Thus, the calculation accuracy is greatly reduced . Many scholars at home and abroad have compared the effects of monitoring drought events with different indicators , Based on precipitation data SPI The index calculation is simple , Good stability , It has the advantages of multiple time scales and space-time comparability , It has been widely used in drought monitoring at home and abroad .
In this paper SPI The index makes a spatiotemporal analysis of the drought situation in a certain region

SPI The index is McKee When assessing the drought situation in Colorado, USA . For a certain region , The precipitation in a certain period of time generally fluctuates regularly . Based on this fact , It is considered that if the precipitation in a certain region in a certain period of time is higher than the average precipitation for many years in that period ( It can be considered as normal precipitation ) Less to a certain extent , It is thought that drought occurred in this region during this period ; conversely ﹐ If it is too much to a certain extent ﹐ Then flood occurs . It is assumed that the change of precipitation follows gamma Distribution , The cumulative frequency distribution of precipitation is transformed into standard normal distribution by mathematical method ﹐ Finally, we can get SPI.SPI The index is dimensionless ﹑ The characteristics of Standardization ﹐ Be able to compare different regions 、 Precipitation levels in different periods , It can better reflect the intensity and duration of drought .SPI The index has the advantage of multiple time scales , It can not only reflect the change of precipitation in a short time scale , It can also reflect the evolution of water resources in a long time scale . In this paper, we calculate 1 Months 、3 Months 、6 Monthly sum 12 Months SPI value ﹐ Use them separately SPI,SPI3.SPI6 and SPI12 representative .

The principle is as follows ：

according to SPI Normal distribution curve of index , According to the classification standard of the National Climate Center, the drought and flood grades , And determine the corresponding SPI Exponential limit ( surface 1).

The experimental results are as follows

This area is not likely to be dry It's all wet ha-ha

Sort out SPI The trend chart is as follows , The selected time scale is 3 Months ：

Import ArcMap Changeable color

Code ：

clear all;
[a,R]=geotiffread('C:\Users\yuanyuan\Desktop\ raster data SPI Calculation \ Precipitation data of a certain area \2001_01.tif'); % First import projection information 
info=geotiffinfo('C:\Users\yuanyuan\Desktop\ raster data SPI Calculation \ Precipitation data of a certain area \2001_01.tif');% First import projection information 
[m,n]=size(a);

datasum=[];   % Generate a number of pixels * Matrix of months 
for j=2001:2020
    for i=1:12      % Start and end month 
        if i<10
        filename=['C:\Users\yuanyuan\Desktop\ raster data SPI Calculation \ Precipitation data of a certain area \',int2str(j),'_0',int2str(i),'.tif'];     % Read in file name    
        data=importdata(filename);  % Import data 
        data=reshape(data,m*n,1);   %reshape  Change the matrix form to m*n That's ok 、1 Column 
        datasum=[datasum  data];          % Put the monthly data into datasum Each column of 
    
        else 
        filename=['C:\Users\yuanyuan\Desktop\ raster data SPI Calculation \ Precipitation data of a certain area \',int2str(j),'_',int2str(i),'.tif'];     % Read in file name    
        data=importdata(filename);  % Import data 
        data=reshape(data,m*n,1);   %reshape  Change the matrix form to m*n That's ok 、1 Column 
        datasum=[datasum  data];          % Put the monthly data into datasum Each column of 
    
        end
    end    
end
%%  Data processing ： non-negative 
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 Calculation 
% scale  Time scale 
% nseas  Number of seasons 
scale=[1,3,6,12];        
nseas=12;
for k=1:240
%  One month scale SPI, Other time scales are similar 
    Z(k) = mean(SPI(datasum(:,k),scale(1),nseas)); %  Take the average 
end
figure(2)
subplot(221)
plot(Z)
ylabel('SPI Index ')
title(' The time scale is 1 Months ：SPI_1')

for k=1:240
% 3 Month scale SPI, Other time scales are similar 
    Z(k) = mean(SPI(datasum(:,k),scale(2),nseas)); %  Take the average 
end
figure(2)
subplot(222)
plot(Z)
ylabel('SPI Index ')
title(' The time scale is 3 Months ：SPI_3')

for k=1:240
% 6 Month scale SPI, Other time scales are similar 
    Z(k) = mean(SPI(datasum(:,k),scale(3),nseas)); %  Take the average 
end
figure(2)
subplot(223)
plot(Z)
ylabel('SPI Index ')
title(' The time scale is 6 Months ：SPI_6')

for k=1:240
% 12 Month scale SPI, Other time scales are similar 
    Z(k) = mean(SPI(datasum(:,k),scale(4),nseas)); %  Take the average 
end
figure(2)
subplot(224)
plot(Z)
ylabel('SPI Index ')
title(' The time scale is 12 Months ：SPI_12')
%%
for k=1:240
% 3 Month scale SPI, Other time scales are similar 
    Z1{
    k} = SPI(datasum(:,k),scale(1),nseas); %  Take the average 
end
Z2=cell2mat(Z1');
%  Seek trend chart 
Z2M=mean(Z2);
Trend=reshape(Z2M,[m,n]);
geotiffwrite(strcat(' result /','Trend1.tif'),Trend,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);% Pay attention to modifying the path 

Code download link ：
https://download.csdn.net/download/qq_45047246/86248543