dim_spi_n
Calculates the standardized precipitation index (SPI) by fitting a gamma or a Pearson Type III distribution to monthly precipitation values.
Available in version 6.1.0 and later.
Prototype
function dim_spi_n ( x : numeric, ; float, double nrun [1] : integer, opt : logical, dims [*] : integer ) return_val : float or double
Arguments
xMonthly precipitation of type 'float' or 'double' and any dimensionality. The size of the specified dims must be divisible by 12. Since a distribution is being fit, there should be a 'reasonably' large sample size. At least 30 years of monthly data (360=12*30) is recommended.
nrunA scalar that specifies the number of months over which the standardized precipitation index is to be calculated. Common values are 3, 6, 12, 24,36.
optOptions parameter.
As of NCL version 6.3.0, if opt=True, you can set opt@spi_type = 3 to have this function calculate the standardized precipitation index (SPI) using the Pearson type III distribution.
dimsThe dimension(s) of x to be used to estimate the SPI. Usually, this is the record ('time') dimension.
Return value
The returned SPI will be the same shape, size and type as x.
Description
This function calculates the Standardized Precipitation Index (SPI), a probability (ie., statistical) index that gives a representation of abnormal wetness and dryness. It is an alternative to the (more complicated) physically based Palmer Severe Drought Index (PSDI) which uses a simple water balance model.
In 2009, the World Meteorological Organization (WMO) approved the Lincoln Declaration on Drought Indices (LDDI). The LDDI recommends that "the Standardized Precipitation Index (SPI) be used to characterize the meteorological droughts around the world", in addition to other drought indices that were in use in their service. In support of this recommendation, it was suggested that a "comprehensive user manual" describing the SPI should be developed. The manual provides a description of the index, the computation methods, specific examples of where it is currently being used, the strengths and limitations and mapping capabilities.
Some advantages of the SPI:
- It requires only monthly precipitation.
- It can be compared across regions with markedly different climates.
- The standardization of the SPI allows the index to determine the rarity of a current drought.
- It can be created for differing periods of 1-to-36 months.
A shortcoming of the SPI, as noted by Trenberth et al (2014):
"the SPI are based on precipitation alone and provide a measure only for water supply.
They are very useful as a measure of precipitation deficits or meteorological drought
but are limited because they do not deal with the ET [evapotranspiration] side of the issue."
As of NCL version 6.3.0, this function will calculate the SPI using two possible methods:
- Applying a 2-parameter gamma distribution fit (default).
- Applying a Pearson type III distribution (opt@spi_type = 3).
The default implementation of dim_spi_n uses a 2-parameter gamma distribution fit (dim_gamfit_n) where the shape and scale parameters are maximum likelihood estimates as described in
A Note on the Gamma Distribution
Thom (1958): Monthly Weather Review, pp 117-122.
specifically: eqn 22 for gamma; just above eqn 21
However, there is some variation in the methods used to derive the
SPI. Guttman (1998, 1999) recommends that the Pearson III
distribution be used, which is available in NCL V6.3.0 by
setting opt@spi_type=3. This option uses code made available
at the National Climate Data Center
(NCDC). See spi.f
Note: In 2015, the NCDC was merged with the National Geophysical Data Center (NGDC) and the National Oceanic Data Center (NODC) into the National Centers for Environmental Information (NCEI).
Generally, the Pearson III distribution is likely to give essentially equivalent results to the 2-parameter gamma distribution fit. In some instances, where monthly and seasonal precipitation of zero is common, it will give slightly better results.
Generally, monthly precipitation is not normally distributed so a transformation is performed such that the derived SPI values follow a normal distribution. The SPI is the number of standard deviations that the observed value would deviate from the long-term mean, for a normally distributed random variable. One interpretation of the resultant values is:
[+,-]2.00 and above/below: exceptionally [wet,dry]
[+,-]1.60 to 1.99: extremely [wet,dry]
[+,-]1.30 to 1.59: severely [wet,dry]
[+,-]0.80 to 1.29: moderately [wet,dry]
[+,-]0.51 to 0.79: abnormally [wet,dry]
[+,-]0.50: near normal
An
explanation of the SPI at different lengths and
sample spatial pattterns
over the USA at different run times are available.
More information can be obtained at the ClimateDataGuide.
References:
McKee, T.B., N.J. Doesken, and J. Kleist, 1993.
The relationship of drought frequency and duration ot time scales.
Eighth Conference on Applied Climatology, American Meteorological Society
Jan 17-23, 1993, Anaheim CA, pp. 179-186.
McKee, T.B., N.J. Doesken, and J. Kleist, 1995.
Drought monitoring with multiple time scales.
Ninth Conference on Applied Climatology, American Meteorological Society
Jan 15-20, 1995, Dallas TX, pp. 233-236.
Guttman, N.B., 1998.
Comparing the Palmer Drought Index and the Standardized Precipitation Index.
Journal of the American Water Resources Association, 34(1), 113-121.
Guttman, N.B., 1999. Accepting the Standardized Precipitation Index: A calculation algorithm.
Journal of the American Water Resources Association, 35(2), 311-322.
Trenberth et al (2014)
Global warming and changes in drought
Nature Climate Change 4, 17-22; doi:10.1038/nclimate2067
WMO:
Handbook of Drought Indicators
Lincoln Declaration on Drought Indices
Standardized Precipitation Index User Guide
See Also
Examples
Several applications of the SPI with figures are HERE.
Example 1
Read an ASCII file containing 117 years of monthly precipitation at Boulder, Colorado (source: http://www.esrl.noaa.gov/psd/boulder/). Compute 12-month estimates of the SPI.
diri = "./" fili = "Boulder.precip.1894-2010.txt" ncol = 14 nrow = numAsciiRow(diri+fili) ; contributed.ncl data = asciiread(diri+fili,(/nrow,ncol/), "float") prc = ndtooned(data(:,1:ncol-2)) ; make one dimensional prc@units = "inches" prc@long_name = "Boulder Precipitation" printVarSummary(prc) print("min="+min(prc)+" max="+max(prc)) nprc = dimsizes(prc) ; check size if ((nprc%12).ne.0) then print("prc size must be divisible by 12") ; full 12-month years only exit end if spi = dim_spi_n(prc, 12, False, 0) ; create a yyyymm array for printing purposes year = toint(data(:,0)) nyear = dimsizes(year) yrStrt = year(0) yrLast = year(nyear-1) yyyymm = yyyymm_time(yrStrt, yrLast, "integer") ; contributed.ncl print(yyyymm+sprintf("%8.2f", prc)+sprintf("%8.2f", spi))The output has _FillValue at the beginning (nrun-1) temporal locations.
yyyymm 12
(0) 189401 -999.00
(1) 189402 -999.00
(2) 189403 -999.00
(3) 189404 -999.00
(4) 189405 -999.00
(5) 189406 -999.00
(6) 189407 -999.00
(7) 189408 -999.00
(8) 189409 -999.00
(9) 189410 -999.00
(10) 189411 -999.00
(11) 189412 -0.39
(12) 189501 -0.34
(13) 189502 -0.39
(14) 189503 -0.24
(15) 189504 -0.33
(16) 189505 -0.37
(17) 189506 0.30
[SNIP]
(1391) 200912 0.80
(1392) 201001 0.72
(1393) 201002 0.96
(1394) 201003 1.21
(1395) 201004 0.75
(1396) 201005 0.67
(1397) 201006 0.83
(1398) 201007 1.04
(1399) 201008 1.23
(1400) 201009 1.20
(1401) 201010 0.67
(1402) 201011 0.59
(1403) 201012 0.39
Example 2Same as Example 1 but (i) use a netCDF version of the file and (ii) compute SPI for 3, 6, 12, 24, 36 and 48 lengths.
diri = "./" fili = "Boulder.precip.1894-2010.nc" f = addfile(diri+fili, "r") prc = f->PRECIP ; PRECIP(time) nprc = dimsizes(prc) ; # monthly precipitation values if ((nprc%12).ne.0) then ; % is NCL syntax for the 'mod' function print("prc size must be divisible by 12") ; full 12-month years only exit end if len = (/3, 6, 12, 24, 36, 48 /) klen = dimsizes(len) spi = new( (/klen, nprc/) , typeof(prc), prc@_FillValue) do k=0,klen-1 spi(k,:) = dim_spi_n(prc, len(k), False, 0) ; spi(nlen,nprc) end do print(yyyymm+sprintf("%8.2f", prc) \ +sprintf("%8.2f", spi(0,:))+sprintf("%8.2f", spi(1,:)) \ +sprintf("%8.2f", spi(2,:))+sprintf("%8.2f", spi(3,:)) \ +sprintf("%8.2f", spi(4,:))+sprintf("%8.2f", spi(5,:)) )The output has _FillValue at the beginning (nrun-1) temporal locations.
yyyymm prc 3 6 12 24 36 48
(0) 189401 0.16 -999.00 -999.00 -999.00 -999.00 -999.00 -999.00
(1) 189402 0.82 -999.00 -999.00 -999.00 -999.00 -999.00 -999.00
(2) 189403 1.40 -0.45 -999.00 -999.00 -999.00 -999.00 -999.00
(3) 189404 2.30 -0.15 -999.00 -999.00 -999.00 -999.00 -999.00
(4) 189405 4.50 0.41 -999.00 -999.00 -999.00 -999.00 -999.00
(5) 189406 0.80 0.13 -0.14 -999.00 -999.00 -999.00 -999.00
(6) 189407 3.08 0.71 0.38 -999.00 -999.00 -999.00 -999.00
[SNIP]
(1398) 201007 2.31 0.71 1.39 1.04 1.42 0.98 0.92
(1399) 201008 1.07 0.67 1.03 1.23 1.16 0.88 0.88
(1400) 201009 0.25 -0.54 0.30 1.20 0.95 0.68 0.78
(1401) 201010 0.95 -1.17 -0.13 0.67 0.89 0.60 0.46
(1402) 201011 0.61 -1.29 -0.14 0.59 0.94 0.60 0.44
(1403) 201012 0.48 -0.76 -0.87 0.39 0.81 0.40 0.14
Example 3
Consider prc(time,lat,lon) and PRC(lat,lon,time), calculate a 12-month SPI. Note that the dims changes depending upon order of the variable dimensions.
prc = f->prc ; prc(time,lat,lon) => (0,1,2) run = 12 spi = dim_spi_n(prc, run, False, 0) => (ntim,nlat,mlon) PRC = f->PRC ; PRC(lat,lon,time) => (0,1,2) SPI = dim_spi_n(PRC, run, False, 2) => (nlat,mlon,ntim) ; optional ... add meta data spi@long_name = "spi: run="+run copy_VarCoords(prc, spi) ; contributed.ncl SPI@long_name = "SPI: run="+run copy_VarCoords(PRC, SPI) ; contributed.ncl
Example 4
An example of the SPI derived from Global Precipitation Climatology Project (GPCP) data spanning 1979-2010 is here.
Example 5
Consider a file which contains time spanning more than full years. The last full year with data is 2013. Use the ind function to directly read in only the year-months specified.
diri = "./" fili = "precip.mon.mean.v22.197901-201405.nc" ; 5 additional months f = addfile(diri+fili,"r") YYYYMM = cd_calendar(f->time, -1) ; ALL dates on file ymStrt = 197901 ymLast = 201312 ; last full year ntStrt = ind(YYYYMM.eq.ymStrt) ; index of ymStrt ntLast = ind(YYYYMM.eq.ymLast) ; index of ymLast delete(YYYYMM) ; no longer needed ; read full year subset using appropriate index values prc = f->precip(ntStrt:ntLast,:,:) ; PREC(time,lat,lon) printVarSummary(prc) dimp = dimsizes(prc) ; # monthly precipitation values ntim = dimp(0) ; 420 nlat = dimp(1) mlon = dimp(2) len = (/3, 6, 12, 24, 36, 48 /) klen = dimsizes(len) spi = new((/klen,ntim,nlat,mlon/) ,float,prc@_FillValue) do k=0,klen-1 spi(k,:,:,:) = dim_spi_n(prc, len(k), False, 0) end do copy_VarCoords(prc,spi(0,:,:,:)) spi@long_name = "SPI" spi!0 = "len" spi&len = len printVarSummary(spi) ; [len | 6] x [time | 420] x [lat | 72] x [lon | 144] ; print values at one location; Change _FillValue from -9.96921e+36f to -999.9 for nicer print LAT = 45.5 LON = 292.5 prc@_FillValue = -999.0 spi@_FillValue = -999.0 yyyymm = cd_calendar(f->time, -1) ; time@units = "days since ..." print(yyyymm+sprintf("%8.2f", prc(:,{LAT},{LON})) \ +sprintf("%8.2f", spi(0,:,{LAT},{LON})) \ +sprintf("%8.2f", spi(1,:,{LAT},{LON})) \ +sprintf("%8.2f", spi(2,:,{LAT},{LON})) \ +sprintf("%8.2f", spi(3,:,{LAT},{LON})) \ +sprintf("%8.2f", spi(4,:,{LAT},{LON})) \ +sprintf("%8.2f", spi(5,:,{LAT},{LON})) )The output:
Variable: spi
Type: float
Total Size: 104509440 bytes
26127360 values
Number of Dimensions: 4
Dimensions and sizes: [len | 6] x [time | 420] x [lat | 72] x [lon | 144]
Coordinates:
len: [3..48]
time: [65378..78131]
lat: [88.75..-88.75]
lon: [1.25..358.75]
Number Of Attributes: 2
long_name : SPI
_FillValue : -9.96921e+36
(0) 197901 6.87 -999.00 -999.00 -999.00 -999.00 -999.00 -999.00
(1) 197902 2.69 -999.00 -999.00 -999.00 -999.00 -999.00 -999.00
(2) 197903 4.46 1.97 -999.00 -999.00 -999.00 -999.00 -999.00
(3) 197904 4.16 1.10 -999.00 -999.00 -999.00 -999.00 -999.00
(4) 197905 4.80 1.66 -999.00 -999.00 -999.00 -999.00 -999.00
(5) 197906 2.43 0.76 1.81 -999.00 -999.00 -999.00 -999.00
(6) 197907 2.63 -0.17 0.67 -999.00 -999.00 -999.00 -999.00
(7) 197908 4.36 -0.57 1.09 -999.00 -999.00 -999.00 -999.00
(8) 197909 3.57 0.22 0.70 -999.00 -999.00 -999.00 -999.00
(9) 197910 3.47 0.49 0.26 -999.00 -999.00 -999.00 -999.00
(10) 197911 3.38 -0.02 -0.37 -999.00 -999.00 -999.00 -999.00
(11) 197912 2.98 -0.23 -0.13 1.03 -999.00 -999.00 -999.00
(12) 198001 1.77 -0.85 -0.21 0.18 -999.00 -999.00 -999.00
(13) 198002 1.17 -1.91 -0.98 -0.10 -999.00 -999.00 -999.00
(14) 198003 3.44 -1.25 -0.92 -0.29 -999.00 -999.00 -999.00
(15) 198004 3.08 -0.53 -0.90 -0.51 -999.00 -999.00 -999.00
(16) 198005 1.43 -0.51 -1.52 -1.14 -999.00 -999.00 -999.00
(17) 198006 2.76 -1.24 -1.68 -1.03 -999.00 -999.00 -999.00
(18) 198007 4.78 -0.66 -0.85 -0.60 -999.00 -999.00 -999.00
(19) 198008 2.62 -0.11 -0.51 -0.90 -999.00 -999.00 -999.00
(20) 198009 4.24 0.80 -0.41 -0.86 -999.00 -999.00 -999.00
(21) 198010 3.24 0.02 -0.45 -0.92 -999.00 -999.00 -999.00
(22) 198011 3.64 0.24 0.10 -0.97 -999.00 -999.00 -999.00
(23) 198012 3.12 -0.18 0.19 -0.84 0.17 -999.00 -999.00
(24) 198101 2.06 -0.55 -0.36 -0.78 -0.43 -999.00 -999.00
(25) 198102 3.72 -0.25 -0.00 -0.29 -0.30 -999.00 -999.00
(26) 198103 2.58 -0.23 -0.32 -0.46 -0.49 -999.00 -999.00
(27) 198104 2.77 0.14 -0.31 -0.53 -0.66 -999.00 -999.00
(28) 198105 3.49 -0.10 -0.23 -0.10 -0.80 -999.00 -999.00
(29) 198106 3.95 0.24 -0.04 0.11 -0.63 -999.00 -999.00
(30) 198107 3.37 0.32 0.26 -0.14 -0.52 -999.00 -999.00
(31) 198108 5.45 1.37 0.70 0.33 -0.42 -999.00 -999.00
(32) 198109 5.52 2.13 1.57 0.57 -0.20 -999.00 -999.00
(33) 198110 5.91 2.16 2.11 1.06 0.09 -999.00 -999.00
(34) 198111 2.80 1.29 1.79 1.00 0.02 -999.00 -999.00
(35) 198112 4.63 0.82 1.57 1.10 0.22 0.69 -999.00
(36) 198201 3.64 0.31 1.58 1.38 0.43 0.38 -999.00
(37) 198202 3.36 0.96 1.32 1.35 0.68 0.44 -999.00
(38) 198203 2.70 0.35 0.82 1.42 0.59 0.28 -999.00
(39) 198204 3.13 0.20 0.31 1.55 0.61 0.20 -999.00
(40) 198205 0.75 -1.21 -0.09 1.00 0.52 -0.16 -999.00
(41) 198206 3.79 -1.02 -0.41 0.92 0.60 -0.05 -999.00
(42) 198207 2.46 -1.88 -1.03 0.74 0.35 -0.06 -999.00
(43) 198208 3.63 -0.28 -1.14 0.44 0.48 -0.13 -999.00
(44) 198209 3.41 -0.40 -1.09 0.09 0.39 -0.14 -999.00
(45) 198210 1.48 -0.61 -1.74 -0.74 0.19 -0.33 -999.00
(46) 198211 5.10 -0.19 -0.36 -0.34 0.39 -0.17 -999.00
(47) 198212 2.72 -0.41 -0.55 -0.67 0.33 -0.21 0.28
(48) 198301 3.49 0.40 -0.14 -0.69 0.49 -0.04 0.03
(49) 198302 3.15 -0.02 -0.17 -0.74 0.41 0.14 0.07
[SNIP]
(413) 201306 4.98 0.66 -0.52 -0.98 -0.59 -0.08 -0.20
(414) 201307 3.82 1.50 0.05 -0.45 -0.57 0.01 -0.28
(415) 201308 4.25 1.52 0.58 -0.20 -0.79 0.24 -0.16
(416) 201309 4.34 1.20 1.24 -0.16 -0.71 0.15 0.01
(417) 201310 1.65 0.07 0.97 -0.82 -0.91 -0.04 -0.20
(418) 201311 3.52 -0.39 0.50 -0.36 -0.73 -0.11 -0.19
(419) 201312 2.88 -0.88 -0.12 -0.42 -0.74 -0.41 -0.27
Example 6
Same as Example 1 but specify a Pearson-III vi the opt@spi_type option
fili = "Boulder.precip.1894-2010.txt" ncol = 14 nrow = numAsciiRow(diri+fili) ; contributed.ncl data = asciiread(diri+fili,(/nrow,ncol/), "float") prc = ndtooned(data(:,1:ncol-2)) ; make one dimensional prc@units = "inches" prc@long_name = "Boulder Precipitation" printVarSummary(prc) print("min="+min(prc)+" max="+max(prc)) nprc = dimsizes(prc) ; check size if ((nprc%12).ne.0) then print("prc size must be divisible by 12") ; full 12-month years only exit end if opt = True opt@spi_type = 3 ; calculate using Pearson III distribution spi = dim_spi_n(prc, 12, opt, 0) ; create a yyyymm array for printing purposes year = toint(data(:,0)) nyear = dimsizes(year) yrStrt = year(0) yrLast = year(nyear-1) yyyymm = yyyymm_time(yrStrt, yrLast, "integer") ; contributed.ncl print(yyyymm+sprintf("%8.2f", prc)+sprintf("%8.2f", spi))The output has _FillValue at the beginning (nrun-1) temporal locations.
yyyymm 12
(0) 189401 -999
(1) 189402 -999
(2) 189403 -999
(3) 189404 -999
(4) 189405 -999
(5) 189406 -999
(6) 189407 -999
(7) 189408 -999
(8) 189409 -999
(9) 189410 -999
(10) 189411 -999
(11) 189412 -0.35
(12) 189501 -0.30
(13) 189502 -0.37
(14) 189503 -0.19
(15) 189504 -0.36
(16) 189505 -0.45
(17) 189506 0.24
[SNIP]
(1391) 200912 0.80
(1392) 201001 0.72
(1393) 201002 0.94
(1394) 201003 1.19
(1395) 201004 0.72
(1396) 201005 0.63
(1397) 201006 0.81
(1398) 201007 1.05
(1399) 201008 1.25
(1400) 201009 1.20
(1401) 201010 0.68
(1402) 201011 0.60
(1403) 201012 0.41