帕默尔干旱指数通常指帕默尔干旱强度指数,即PDSI,是一种在气象气候农业水文水资源等领域广泛应用的干旱指数,由美国气象学家韦恩·帕默尔(Wayne Palmer)于1965年提出。

帕默尔干旱指数(PDSI)是一个基于水量供需关系的干旱指数,在当地水分供不应求时即为干旱,否则为湿润。水分供给量相对容易求得,通常可以降水量代替;而水分需求量的推算则较为复杂,因其涉及到受到气温、土壤性质、土地利用等因素影响的蒸散发及土壤水分变化等方面。对于这一问题,帕默尔提出了“当前情况下的气候适宜”(climatically appropriate for existing conditions,CAFEC)的概念,定义了“气候适宜降水量”来作为水分需求量,并以实际降水量与其差值来界定水分盈亏状况。PDSI不仅能考虑当前的水分供需状况,还能考虑前期干湿状况及其持续时间来对当前干旱状况的影响,物理意义明确,是一个能较为客观合理地定量描述干旱的干旱指数。

PDSI值通常介于-4到4之间,当其值大于0时则为湿润,反之为干旱,不同大小的值反映不同的干旱或湿润等级。

下面是计算PDSI干旱指数的R代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#' Calculation of (sc)PDSI
#'
#' Calculates monthly scPDSI series from temperature and precipitation data.
#' Uses binaries compiled from official University of Nebraska C++ Code.
#' @details This function transforms the input climate data into the necessary
#' files for the PDSI binary executable and saves them to a temp directory.
#' The binary is called and the resulting files of interest for PDSI and
#' scPDSI are read in again and returned.
#'
#' For details on the algorithm, see the comments in the C++ source file. For
#' reference, the original source code (scpdsi-orig.cpp) is put in the
#' toplevel installation directory of the package, you may find it using
#' \code{system.file(package = "pdsi")}. For building the Linux binaries, a
#' modified version of the original source code (scpdsi.cpp, string constants
#' cast to \code{(char *)} for compatibility with modern g++, Windows-specific
#' code commented out) is also distributed with the package in the same
#' directory. To build the binary on Linux, use
#' \code{pdsi::build_linux_binary()}.
#' @param awc Available soil water capacity (in mm)
#' @param lat Latitude of the site (in decimal degrees)
#' @param climate \code{data.frame} with monthly climate data consisting of 4
#' columns for year, month, temperature (deg C), and precipitation (mm)
#' @param start Start year for PDSI calculation
#' @param end End year for PDSI calculation
#' @param mode one of c("both", "pdsi", "scpdsi")
#' @param verbose should stderr/stdout be shown?
#' @return For mode "both", a \code{list} of two \code{data.frames},
#' one holding the standard PDSI, one holding the scPDSI. For modes
#' "pdsi" or "scpdsi" only the respective \code{data.frame}.
#' @references Methodology based on Research Paper No. 45; Meteorological
#' Drought; by Wayne C. Palmer for the U.S. Weather Bureau, February 1965.
#' @keywords utils
#' @examples
#' library(bootRes)
#' data(muc.clim)
#' pdsi(120, 50, muc.clim, 1960, 2000)
#' @importFrom bootRes pmat
#' @import digest
#' @export
pdsi <- function(awc, lat, climate, start, end, mode = "both",
verbose = TRUE) {

## check the system we are on
the_system <- Sys.info()["sysname"]

## create temp dir
tdir <- paste(getwd(), "/", digest(Sys.time()), sep = "")
dir.create(tdir)

## truncate and reformat climate data
climate_start <- which(climate[,1] == start-1)[1]
climate_end <- which(climate[,1] == end)[12]
climate <- climate[climate_start:climate_end,]
climate_reform <- pmat(climate, start = 1, end = 12)

## split in temp and prec
pmat_temp <- climate_reform[,1:12]
pmat_prec <- climate_reform[,13:24]

## write to files
temp_path <- file.path(tdir, "monthly_T")
prec_path <- file.path(tdir, "monthly_P")
write.table(pmat_temp, temp_path, col.names = F, quote = F)
write.table(pmat_prec, prec_path, col.names = F, quote = F)

## calculate mean values and write to files
normal_temp <- round(t(as.vector(colMeans(pmat_temp))), 3)
normal_prec <- round(t(as.vector(colMeans(pmat_prec))), 3)
normal_temp_path <- file.path(tdir, "mon_T_normal")
normal_prec_path <- file.path(tdir, "mon_P_normal")
write.table(normal_temp, normal_temp_path, col.names = F, quote = F,
row.names = F)
write.table(normal_prec, normal_prec_path, col.names = F, quote = F,
row.names = F)

## write parameter files to tempdir
params <- t(c(awc, lat))
param_path <- file.path(tdir, "parameter")
write.table(params, param_path, col.names = F, quote = F,
row.names = F)

## run executable (depending on platform)
if (the_system == "Windows") {
exec_path <- file.path(system.file(package = "pdsi"), "exec", "sc-pdsi.exe")
} else {
if (the_system == "Linux") {
exec_path <- file.path(system.file(package = "pdsi"), "scpdsi.o")
if (!file.exists(exec_path))
stop("You need to build the binary first. On a Linux machine with recent g++ installed, call function `pdsi::build_linux_binary()`.")
} else {
if (the_system == "Darwin") {
exec_path <- file.path(system.file(package = "pdsi"), "exec", "pdsi")
} else {
stop("Unsupported OS.")
}
}
}

oldwd <- getwd()
setwd(tdir)

cmd <- paste(exec_path, " -m -i", shQuote(tdir), start, end)
system(cmd, ignore.stdout = !verbose, ignore.stderr = !verbose)

setwd(oldwd)

## read (sc)PDSI in again and return it
scpdsi_path <- file.path(tdir, "monthly", "self_cal", "PDSI.tbl")
pdsi_path <- file.path(tdir, "monthly", "original", "PDSI.tbl")

if (any(c("scpdsi", "both") == mode)) {
scPDSI <- read.fwf(scpdsi_path, c(5, rep(7, 12)))
colnames(scPDSI) <- c("YEAR", toupper(month.abb))
}

if (any(c("pdsi", "both") == mode)) {
PDSI <- read.fwf(pdsi_path, c(5, rep(7, 12)))
colnames(PDSI) <- c("YEAR", toupper(month.abb))
}

unlink(tdir, recursive = TRUE)

if (mode == "both") {
list(PDSI, scPDSI)
} else {
if (mode == "pdsi") {
PDSI
} else {
if (mode == "scpdsi") {
scPDSI
} else {
stop("`mode` has to be one of 'pdsi', 'scpdsi', or 'both'.")
}
}
}
}