NCL Home > Documentation > Functions > Drought

# 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

x

Monthly 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.

nrun

A 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.

opt

Options 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.

dims

The 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:

1. Applying a 2-parameter gamma distribution fit (default).

2. 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

```

## 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 2

Same 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

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

```