SPI干旱指数[关闭]

Question

我有以下数据：

2001-01,2020-12的pcp#降水月度数据

我使用以下方法来计算 SPI 干旱指数

library(SPEI)

spi1 <- spi(pcp,1,kernel = list(type = "rectangular", shift = 0),distribution = "Gamma", fit = "ub-pwm", na.rm = FALSE,ref.start=NULL, ref.end=NULL, x=FALSE, params=NULL)

第一个问题： 当我 plot(spi1) 我在 Y 轴上得到我不想要的 SPEI，我想要的是 SPI， 第二： 如何分别绘制每个月，例如当你调用 spi1 时，它会给你每个月的索引值，我想为每个月绘制它

Answer 1

对于第一个答案，你可以重写函数plot.spei

plot.spei <- 
function (x, ...) 
{
    ## label <- ifelse(as.character(x$call)[1] == "spei", "SPEI", 
    ##     "SPI")

    ser <- ts(as.matrix(x$fitted[-c(1:x$scale), ]), end = end(x$fitted), 
        frequency = frequency(x$fitted))
    ser[is.nan(ser - ser)] <- 0
    se <- ifelse(ser == 0, ser, NA)
    tit <- dimnames(x$coefficients)[2][[1]]
    if (start(ser)[2] == 1) {
        ns <- c(start(ser)[1] - 1, 12)
    }
    else {
        ns <- c(start(ser)[1], start(ser)[2] - 1)
    }
    if (end(ser)[2] == 12) {
        ne <- c(end(ser)[1] + 1, 1)
    }
    else {
        ne <- c(end(ser)[1], end(ser)[2] + 1)
    }
    n <- ncol(ser)
    if (is.null(n)) 
        n <- 1
    par(mar = c(4, 4, 2, 1) + 0.1)
    if (n > 1 & n < 5) 
        par(mfrow = c(n, 1))
    if (n > 1 & n >= 5) 
        par(mfrow = c({
            n + 1
        }%/%2, 2))
    for (i in 1:n) {
        datt <- ts(c(0, ser[, i], 0), frequency = frequency(ser), 
            start = ns, end = ne)
        datt.pos <- ifelse(datt > 0, datt, 0)
        datt.neg <- ifelse(datt <= 0, datt, 0)
        plot(datt, type = "n", xlab = "", main = tit[i], ...)
        if (!is.null(x$ref.period)) {
            k <- ts(5, start = x$ref.period[1, ], end = x$ref.period[2, 
                ], frequency = 12)
            k[1] <- k[length(k)] <- -5
            polygon(k, col = "light grey", border = NA, density = 20)
            abline(v = x$ref.period[1, 1] + (x$ref.period[1, 
                2] - 1)/12, col = "grey")
            abline(v = x$ref.period[2, 1] + (x$ref.period[2, 
                2] - 1)/12, col = "grey")
        }
        grid(col = "black")
        polygon(datt.pos, col = "blue", border = NA)
        polygon(datt.neg, col = "red", border = NA)
        lines(datt, col = "dark grey")
        abline(h = 0)
        points(se, pch = 21, col = "white", bg = "black")
    }
}

然后使用ylab参数

plot(spi1, ylab = "SPI")

如果你想单独绘制它，你可以提取类ts的拟合值并在R中为时间序列对象应用基本绘图。

par(mfrow = c(3, 4))
listofmonths <- split(fitted(spi1), cycle(fitted(spi1)))
names(listofmonths) <- month.abb

require(plyr)
l_ply(seq_along(listofmonths), function(x) {
       plot(x = seq_along(listofmonths[[x]]), y = listofmonths[[x]],
            type = "l", xlab = "", ylab = "SPI")
       title(names(listofmonths)[x])
   })

你也可以尝试这些类型的情节

monthplot(fitted(spi1), labels = month.abb, cex.axis = 0.8)
boxplot(fitted(spi1) ~ cycle(fitted(spi1)), names = month.abb, cex.axis = 0.8)