-4

次のデータがあります。

pcp # 2001-01,2020-12 の降水量の月間データ

以下を使用して、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)

最初の質問: プロット (spi1) すると、Y 軸に SPEI が表示されますが、これは望ましくありません。必要なのは SPI です。2 番目: 毎月個別にプロットする方法 (たとえば、spi1 を呼び出すと、インデックス値が得られます)月ごとに、月ごとにプロットしたい

4

1 に答える 1

2

最初の答えについては、関数を書き直すことができます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)
于 2013-07-20T06:22:52.253 に答える