干旱指数新宠—SPEI的C语言版、R语言版与单独执行软件SPEI Calculator

hillside

      标题原冠以“国内首发”，是有感于国内论文往往很少明确说明或愿意提供使用的程序与软件，而正这是学习者非常需要的。曾建议有兴趣的网友百度一下或狗狗之。然而，我自己狗狗多页之后，发现国内对于R包已有介绍（见附录），仅此一条消息就将我由“首发”发落至“二手”。至于SPEI Calculator我也不知是首发还是二手。幸好我并不在意首发权，重要的是分享精神。

      以下介绍内容摘自2013年5月刚刚出版的《基于SPEI的云南省夏玉米生长季干旱时空特征分析》（《资源科学》）：
      “目前常用于农业干旱监测、预报的干旱指数有帕尔默干旱强度指数（PDSI）、标准化降水指数（SPI）、水分亏缺指数（CWDI）、地表湿润指数、气象干旱指数（CI）等。PDSI、CWDI和地表湿润指数是基于降水、径流和蒸发量等干旱影响因子提出的，很好地解决了区域干旱监测、预报，但是缺乏时空比较的有效性。SPI能够反映多时间、多尺度标量的降水异常，但忽略了气温上升引起的蒸散量变化对干旱的影响。实证研究表明，气温上升明显影响干旱程度。因此，Vicente-Serrano等提出标准化降水蒸散指数（Standardized PrecipitationEvapotranspiration Index，SPEI），该指数基于降水和蒸散，结合了PDSI对温度的灵敏性和SPI多时空的特性，目前已成为国内外学者研究分析干旱演变趋势新的理想指标。国外研究侧重于气候变化条件下不同时间尺度的干旱监测与农作物生产评估等方面。国内学者也尝试用SPEI分析干旱与气候因子关系、干旱演变趋势、干旱时空特征等。”
SPEI: Calculation of the Standardised Precipitation-Evapotranspiration Index

A set of functions for computing potential evapotranspiration and several widely used drought indices including the Standardized Precipitation-Evapotranspiration Index (SPEI).

Version:	1.4
Depends:	lmomco
Published:	2013-05-09
Author:	Santiago Beguería and Sergio M. Vicente-Serrano.

About the SPEI
Drought monitoring and analysis based on climatic indices

Drought is a major cause of agricultural, economic and environmental damage. Drought effects are apparent after a long period with a shortage of precipitation, making it very difficult to determine their onset, extent and end. Thus, it is hard to objectively quantify the characteristics of drought episodes in terms of their intensity, magnitude, duration and spatial extent. Much effort has been devoted to developing techniques for drought analysis and monitoring. Among these, the definition of quantitative indices is the most widespread approach, but subjectivity in the definition of drought has made it very difficult to establish a unique and universal drought index. Most studies related to drought analysis and monitoring systems have been conducted using either i) the Palmer Drought Severity Index (PDSI) (Palmer, 1965), based on a soil water balance equation, or ii) the Standardised Precipitation Index (SPI; McKee et al., 1993), based on a precipitation probabilistic approach.

注：我将两个程序打包在一个文件夹中。 附件中的“SPEI的R语言版”内置了C语言版源程序，可单独使用。
                                           试用体会：该程序（自带实例数据）操作十分简便

     我已经试用了C语言版与可执行版SPEI Calculator。在windows系统命令行下进入存放目录，直接运行SPEI Calculator的执行文件spei.exe即可生成输出文件。形式如：spei 时间间隔（直接输入数字，如1、3等，含义为月份）源文件名 输入文件名。C语言版在Linux或Cygwin下用GCC编译即可，运行方式与SPEI Calculator相同。在计算尺度上，通常介绍可以计算3、6、9、12、24、48月多种尺度，但据笔者判断，实际上可计算任意月份跨度的SPEI（原始数据为逐月数据）。

       再在命令行窗口键盘敲入：1    spei  1 tampa.txt tampa_spei_1.txt                                       ( 计算1个月尺度)或
2    spei 12 tampa.txt tampa_spei_12.txts                                （计算12个月尺度）或            3    spei  ntampa.txt tampa_spei_12.txts                                 （计算n个月尺度）

具体介绍见以下说明文档：
spei
Create a time series of the Standardised Precipitation-Evapotranspiration Index

Description
The program calculates a time series of the Standardised Precipitation-
Evapotransporation Index (SPEI) at a given time interval from an input data file
containing monthly time series of precipitation and mean temperature, plus the
geographic coordinates of the observatory.
Usage

spei [timeInterval] [inputFile] [outputFile]

Arguments
timeInterval A time interval, in months.
inputFile The name of an input file, with extension.
outputFile The name of the output file to be produced, with extension.
Details
The SPEI index is a standardised monthly climatic balance computed as the difference
between the cummulative precipitation and the potential evapotranspiration. Details on
the index calculation and applications can be found in the references below.
The SPEI can be calculated at the monthly scale with timeInterval = 1, or accumulated
at more than one month with timeInterval > 1. Typical values are 1, 3, 6, 12 and 24
months. If the accumulated index is calculated, the starting date of the resulting SPEI
series will be lagged a number of months equal to timeInterval – 1.
The input file (inputFile) can have any extension, but must be a plain text file (ASCII).
The file structure is as follows:
tampa
27.96
1900;01
12
110.70;14.30
...
The first line contains the name of the observatory, and is only used for identification
purposes. The second line is the latitude of the observatory, in degrees. The third line
contains the year and month of the first record in the time series, separated by a semicolon
(;). The fourth line contains the seasonality of the time series, and must be set to 12.
Finally, from the fifth line the data series of monthly precipitation and mean temperature,
separated by a semi-colon (;). The series must be continuous. Gaps and missing-values
are not allowed.
The output file (outputFile) can have any extension. It will be a plain text (ASCII) file,
with the following structure:
tampa
27.959999
1900;12
6
1.456516
...
The first three lines contain the name of the station, latitude and initial date of the SPEI
series. The fourth line contains a value indicating the cummulative parameter used,
timeInterval (six months in the example). The SPEI time series is given from the fifth
line on.
The program is run from the Windows console. The easiest way is to locate the program
and the input file(s) in the same directory. If you need to run the program from a different
location, it might be necessary to modify the path system variable to include the path to
the directory where the program was installed.
It is easy to create a batch script for automating the calculation of the SPEI over a large
number of observatories or for several accumulated periods.
A hint on the usage of the program is obtained if the spei is invoked with no arguments or
with a wrong number of arguments.
See Also
spi program.
Examples
spi 1 tampa.txt tampa_spei_1.txt
spi 12 tampa.txt tampa_spei_12.txt
The above lines calculate the monthly SPEI and the 12-months cummulative SPEI time
series for Tampa (Florida).


附录1：http://www.biostatistic.net/thread-85413-1-1.html相关网页介绍

Drought indices(SPEI) Drought indices()所属R语言包：SPEI                                         Calculation of the Standardized Precipitation-Evapotranspiration Index (SPEI) and the Standardized Precipitation Index (SPI).                                          计算的标准化降水的蒸散指数（SPEI）和标准化降水指数（SPI）。                                          译者：生物统计家园网 机器人LoveR 描述----------Description---------- Given a time series of the climatic water balance (precipitation minus potential evapotranspiration), gives a time series of the Standardized Precipitation-Evapotranspiration Index (SPEI). 给定一个时间序列的气候水分平衡的（降水量减去潜在蒸散量），给出了一个时间序列的标准化降水的蒸散指数（SPEI）。 用法----------Usage---------- spei(data, scale, kernel = list(type = 'rectangular', shift = 0),         distribution = 'log-Logistic', fit = 'ub-pwm', na.rm = FALSE,         ref.start=NULL, ref.end=NULL, x=FALSE, ...) spi(data, scale, kernel = list(type = 'rectangular', shift = 0),         distribution = 'Gamma', fit = 'ub-pwm', na.rm = FALSE,         ref.start=NULL, ref.end=NULL, x=FALSE, ...) 参数----------Arguments---------- 参数：data a vector, matrix or data frame with time ordered values of precipitation (for the SPI) or of the climatic balance precipitation minus potential evapotranspiration (for the SPEI). 一个向量，矩阵或数据框随时间沉淀（为SPI），或气候的平衡降水量减去潜在蒸散量（SPEI）有序值。 参数：scale an integer, representing the time scale at which the SPEI / SPI will be computed. 的整数，相当于在哪些SPEI / SPI将被计算的时间尺度。 参数：kernel optional, a list defining the type of kernel used for computing the SPEI / SPI at scales higher than one. Defaults to unshifted rectangular kernel. 可选的，一个列表定义类型的内核，用于计算SPEI / SPI高于一个尺度。默认非移动的矩形内核。 参数：distribution optional, name of the distribution function to be used for computing the SPEI / SPI (one of 'log-Logistic', 'Gamma' and 'PearsonIII'). Defaults to 'log-Logistic'. 可选的，要用于计算SPEI / SPI（之一log - MF的“，”伽马“和”PearsonIII）的分布函数的名称。默认的log - MF。 参数：fit optional, name of the method used for computing the distribution function parameters (one of 'ub-pwm', 'pp-pwm' and 'max-lik'). Defaults to 'ub-pwm'. 可选的，用于计算的分布函数的参数（UB-pwm的，页pwm的和最大力）的方法名称。默认为“UB-PWM。 参数：na.rm optional, a logical value indicating whether NA values should be stripped from the computations. Defaults to FALSE, i.e. no NA are allowed in the data. 可选的，一个逻辑值，该值指示NA值是否应该从计算中被剥离。默认为false，即不NA允许的数据。 参数：ref.start optional, starting point of the reference period used for computing the index. Defaults to NULL, indicating that the first value in data will be used as starting point. 可选的，用于计算索引的基准周期的起点。默认为空，表示数据中的第一个值将被用来作为出发点。 参数：ref.end optional, ending point of the reference period used for computing the index. Defaults to NULL, indicating that the last value in data will be used as ending point. 可选的，用于计算索引的基准周期的结束点。默认为空，表示的终点，数据将被用作最后的值。 参数：x optional, a logical value indicating wether the data used for fitting the model should be kept. Defaults to FALSE. 可选的，一个逻辑值，表示天气的数据拟合模型，应予保留。默认为false。 参数：... other parameters to be passed. 其他参数来传递。 Details 详细信息----------Details---------- The spei and spi functions allow computing the SPEI and the SPI indices. These are widely used for drought identification and monitoring. Both functions are identical (in fact, spi is just a wrapper for spei), but they are kept separated for clarity. Basically, the functions standardize a variable following a log-Logistic distribution function (i.e., they transform it to a standard Gaussian variate with zero mean and standard deviation of one). The input variable is a time ordered series of precipitation values for spi, and a series of the climate water balance (precipitation minus potential evapotranspiration) for spei. When used with the default options, they yield values of both indices as defined in the references given below. spei和spi函数允许SPEI和SPI指数计算。这些被广泛用于干旱的识别和监控。这两个函数是相同的（其实，spi是只是一个包装spei），但他们分开为清晰。基本上，功能规范变量的log - MF分布函数（即，他们把它转换成一个标准的高斯变量的均值为零，标准偏差为1）。输入变量是一个时间顺序系列的降水值spi和一系列的气候水分平衡的（降水量减去潜在蒸散量）为spei。当使用默认的选项，他们得到的数值，这两个指数的定义在下面给出的参考。 The SPEI and the SPI were defined for monthly data. Since the PDFs of the data are not homogenous from month to month, the data is split into twelve series (one for each month) and independent PDFs are fit to each series. If data is a vector or a matrix it will be treated as a sequence of monthly values starting in January. If it is a (univariate or multivariate) time series then the function cycle will be used to determine the position of each observation within the year (month), allowing the data to start in a month other than January. SPEI和SPI被定义为月度数据。由于PDF文件中的数据是不是同质的每个月，数据分割成十二个系列（每月）和独立的PDF文件，适合各系列。如果data是序列年1月开始每月的值将被视为一个向量或矩阵。如果它是一个（单变量或多变量）的时间序列，然后功能cycle将被用于确定每个观测位置年（或月）内，允许数据在一个月以外一月开始。 An important advantage of the SPEI and the SPI is that they can be computed at different time scales. This means that it is possible to incorporate the influence of the past values of the variable in its computation, enabling the index to adapt to the memory of the system under study. The magnitude of this memory is controlled by parameter scale: for example, a value of six would imply that data from the current month and of the past five months will be used for computing the SPEI or SPI value for a given month. By default all past data will have the same weight in computing the index, as it was originally proposed by the authors. Other kernels, however, are available through parameter kernel. The parameter kernel is a list defining the shape of the kernel and a time shift. For more details refer to  the function kern. 磺化聚醚酰亚胺和SPI的一个重要的优点是，它们可以在不同的时间尺度计算。这意味着，它可以将过去的值的影响，在其计算的变量，使索引，以适应所研究的系统的存储器中。该存储器的大小被控制的参数scale：例如，值6将意味着，从当前的月份和过去5个月的数据将被用于计算SPEI或SPI值，对于一个给定的月份。默认情况下，所有的过去的数据将具有相同的权重计算指数，因为它最初是由作者。然而，其他的内核，都可以通过参数kernel。参数kernel是一个列表，限定的内核和一个时间偏移的形状。有关详细信息，请参阅的功能kern。 Following the original definitions spei uses a log-Logistic distribution by default, and spi uses a Gamma distribution. This behaviour can be modified, however, through parameter distribution. 的原始定义spei使用默认的log - MF配送，和spi采用了伽玛分布。然而，这种行为可以被修改，通过参数distribution。 The default method for parameter fitting is based on unbiased Probability Weighted Moments ('ub-pwm'), but other methods can be used through parameter fit. Valid aleternatives are the plotting-position PWM ('pp-pwm') and the maximum likelihood method ('max-lik'). 默认参数拟合方法的基础上公正概率加权矩（UB-PWM），但其他的方法可以用来通过参数fit。有效aleternatives是绘图位置PWM（页pwm的“）和最大似然法（”最大力）。 The default behaviour of the functions is using all the values provided in data for parameter fitting. However, this can be modified with help of parameters ref.start and ref.end. These parameters allow defining a subset of values that will be used for parameter fitting, i.e. a reference period. The functions, however, will compute the values of the indices for the whole data set. For these options to work it is necessary that data will be a time series object. The starting and ending points of the reference period will then be defined as pairs of year and month values, e.g. c(1900,1). 功能的默认行为是使用提供的值在data参数拟合。然而，这可以修改与的参数ref.start和ref.end的帮助。这些参数允许定义的值，将用于为参数拟合，即一个特定期间内的一个子集。的功能，但是，将整个数据集的指数计算的值。这些工作是必要的，data将是一个时间序列对象。基准期的开始和结束点，然后将被定义为对年份和月份的值，例如C（1900,1）。 值----------Value---------- Functions spei and spi return an object of class spei. The generic functions print and summary can be used to obtain summaries of the results. The generic accessor functions coefficients and fitted extract useful features of the object. 功能spei和spi类spei返回一个对象。通用功能中的print和summary可以用来获得汇总结果。一般的访问功能coefficients和fitted的对象中提取有用的功能。 An object of class spei is a list containing at least the following components: 一种的类spei的对象是一个列表，包含至少以下组件： call: the call to spei or spi used to generate the object. 检测：呼叫spei或spi用于生成的对象。 fitted: time series with the values of the Standardized Precipitation-Evapotranspiration Index (SPEI) or the Standardized Precipitation Index (SPI). If data consists of several columns the function will treat each column as independent data, and the result will be a multivariate time series. The names of the columns in data will be used as column names in fitted. 装：时间序列的标准化降水的蒸散指数（SPEI）或标准化降水指数（SPI）的值。如果数据由多个列，该函数将治疗作为独立的数据的每一列，和的结果将是一个多变量时间序列。 data中的列的名称将用作列名在装。 coefficients: an array with the values of the coefficients of the distribution function fitted to the data. The first dimension of the array contains the three (or two) coefficients, the second dimension will typically consist of twelve values corresponding to each month, and the third dimension will be equal to the number of columns (series) in data. If a time scale greater than one has been used then the first elements will have NA value since the kernel can not be applied. The first element with valid data will be the one corresponding to the time scale chosen. 系数：一个数组与嵌合的数据的分布函数的值的系数。阵列的第一维中包含的三个（或两个）系数，第二个维度将通常由12值对应每月，和第三个维度将是相等的列数（系列）在data 。如果大于一的时间尺度已被使用，则在第一个元素将有因为内核不能应用于NA值。有效的数据的第一个元素将是一个对应于所选择的时间尺度。 scale: the scale parameter used to generate the object. 规模：scale参数用于生成对象。 kernel: the parameters and values of the kernel used to generate the object. 内核：用于生成对象的内核的参数和值。 distribution: the distribution function used to generate the object. 分布：用于生成对象的分布函数。 fit: the fitting method used to generate the object. 适合：用于生成对象的拟合方法。 na.action: the value of the na.action parameter used. na.action：的na.action使用的参数的值。 data: if requested, the input data used. 数据：如果被请求的，使用所输入的数据。 注意----------Note---------- Dependencies: the spei function depends on the library lmomco. 依存关系的SpeI功能依赖于库中lmomco。 （作者）----------Author(s)---------- Santiago Beguer铆a and Sergio M. Vicente-Serrano. Maintainer: Santiago Beguer铆a 参考文献----------References---------- 参见----------See Also---------- kern for different kernel functions available. thornthwaite, hargreaves and penman for ways of calculating potential evapotranspiration. summary.spei and print.spei for summaries of spei objects. plot.spei for plotting spei objects. kern不同的内核函数。 thornthwaite，hargreaves和penman潜在蒸散量的计算方法。 summary.spei和print.speispei对象的摘要。 plot.spei：策划spei对象的。 实例----------Examples---------- # Load data[数据加载] data(wichita) # One and tvelwe-months SPEI[和tvelwe个月SPEI] wichita$PET <- thornthwaite(wichita$TMED,37.6475) spei1 <- spei(wichita$PRCP-wichita$PET,1) spei12 <- spei(wichita$PRCP-wichita$PET,12) # not executed[不执行] #spei1[spei1] #summary(spei1)[摘要（spei1）] #coefficients(spei1)[系数（spei1）] par(mfrow=c(2,1)) plot(spei1) plot(spei12,'Wichita') # One and tvelwe-months SPI[一，tvelwe个月的SPI] spi_1 <- spi(wichita$PRCP,1) spi_12 <- spi(wichita$PRCP,12) par(mfrow=c(2,1)) plot(spi_1) plot(spi_12) # Define the properties of the time series with ts()[与日ts（）定义的属性的时间序列] plot(spei(ts(wichita$PRCP-wichita$PET,freq=12,start=c(1980,1)),12)) # Time series not starting in January[时间序列年1月开始] plot(spei(ts(wichita$PRCP-wichita$PET,freq=12,start=c(1980,6)),12)) # Using a particular reference period[使用一个特定参考期] plot(spei(ts(wichita$PRCP-wichita$PET,freq=12,start=c(1980,1)),12,         ref.start=c(1980,1), ref.end=c(2000,1))) # Different kernels[不同的内核] spei24 <- spei(wichita$PRCP-wichita$PET,24) spei24_gau <- spei(wichita$PRCP-wichita$PET,24,kernel=list(type='gaussian',shift=0)) par(mfrow=c(2,1)) plot(spei24,'Rectangular kernel') plot(spei24_gau,'Gaussian kernel') # Several time series at a time[在时间的几个时间序列] data(balance) names(balance) bal_spei12 <- spei(balance,12) plot(bal_spei12)       注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。 注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。 注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。                                                                             转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。

附录2：气象家园相关帖子：
SPE和SPEIStandardised Precipitation-Evapotransporation Index ）的计算程序

龙行天下

在cmd下计算spi，当时间间隔设置为n，输出的数据比原始数据将会少n-1个，请问楼主这个减少代表什么？另外把时间间隔分别设置成1个月，6个月，12个月，输出的结果差别在哪里？

凯旋

SPEI计算出的都是月时间尺度的值，如何转化成季节、年尺度的值，求教！！！

青火流音雪

楼主求回复。我现在做五十年的某站spei指数，一共600月。做spei-1、spei-6、spei-12分别生成了600、598和589个数，如果要评估干旱，直接用生成的这些数画图就可以吗？

善人/jw

感谢LZ分享~~
但不知道LZ有没有研究过SPEI指数，它的适用性如何呢？是对全国还是对某个地区的干旱表征比较好呢？

hillside

善人/jw 发表于 2013-6-4 20:24

                               
登录/注册后可看大图

感谢LZ分享~~
但不知道LZ有没有研究过SPEI指数，它的适用性如何呢？是对全国还是对某个地区的干旱表征比较 ...

还没有用过，只是见过文献都称好。时间略宽裕时打算进行一些分析。

kongfeng0824

做干旱的有福啦。spi和pdsi以及z指数

qixiangren34

还没有用过，学无止境啊。

Jasonluo

有没有人对同一地区用不同干旱指数分析过呢？

nwsxer

刚刚学习那个PDSI指数，你这又来一个SPEI指数，貌似比较厉害的东西，学习了。话说R听容易的啊，基本都可以看明白

zihan6688

不知道好不好学

柠檬安琪拉

这个指数的时间尺度可不可以做到旬的呢？

黎明的夕阳

{:eb502:}