CDT: Climate Data Tools

# Standard Normal Homogeneity Test (Alexandersson and Moberg, 1997)
homog.SNHT.Test <- function(x, h = NULL){
    nl <- length(x)
    ilen <- 1:nl
    ina <- !is.na(x)
    ilen <- ilen[ina]
    x <- x[ina]
    n <- length(x)
    if(n < 3) return(NULL)
    ll <- 1:n
    sdx <- stats::sd(x)
    if(sdx < 1e-12) return(NULL)
    x <- (x - mean(x)) / stats::sd(x)
    Ts <- numeric(n)

    if(!is.null(h)){
        idh <- which(ll/n >= h & ll/n <= 1 - h)
        Tsh <- Ts
    }

    z1 <- ll * (cumsum(x) / ll)^2
    z2 <- rev(ll * (cumsum(rev(x)) / ll)^2)
    Ts[1:(n - 1)] <- z1[-n] + z2[-1]

    if(!is.null(h)){
        Tsh[idh] <- Ts[idh]
        idm <- which.max(Tsh)
        Tm <- max(Tsh)
    }else{
        idm <- which.max(Ts)
        Tm <- max(Ts)
    }

    Vcritic <- fSNHT.Vc(n)
    test <- Tm >= Vcritic[, 2]
    out <- NULL
    if(sum(test) > 0){
        mx.cl <- max(Vcritic[test, 1])
        out <- list(index.change = ilen[idm] + 1, stat.change = Tm, max.confL = mx.cl)
    }
    out
}

###################################################################
# Pettitt Test (based on the rank)

homog.PettittRank.Test <- function(x, h = NULL){
    nl <- length(x)
    ilen <- 1:nl
    ina <- !is.na(x)
    ilen <- ilen[ina]
    x <- x[ina]
    n <- length(x)
    if(n < 3) return(NULL)
    ll <- 1:n
    Sr <- 2 * cumsum(rank(x)) - ll * (n + 1)
    Sr <- abs(Sr)

    if(!is.null(h)){
        idh <- which(ll/n >= h & ll/n <= 1 - h)
        Srh <- numeric(n)
        Srh[idh] <- Sr[idh]
        idm <- which.max(Srh)
        Sm <- max(Srh)
    }else{
        idm <- which.max(Sr)
        Sm <- max(Sr)       
    }

    Pm <- 1 - exp((-6 * Sm^2) / (n^3 + n^2))
    Vcritic <- c(90.0, 92.0, 95.0, 97.5, 99.0, 99.9)
    test <- Pm * 100 >= Vcritic
    out <- NULL
    if(sum(test) > 0){
        mx.cl <- max(Vcritic[test])
        out <- list(index.change = ilen[idm] + 1, stat.change = Sm, max.confL = mx.cl)
    }
    out
}

###################################################################
# CUSUM-type statistics (Gallagher et al.,2013)
# no cropping: powerful around the center
# cropping bounds: powerful around the boundaries

homog.CUSUM.Simple <- function(x, h = NULL){
    nl <- length(x)
    ilen <- 1:nl
    ina <- !is.na(x)
    ilen <- ilen[ina]
    x <- x[ina]
    n <- length(x)
    if(n < 3) return(NULL)
    ll <- 1:n
    s <- sum(x)
    sig <- stats::sd(x)
    if(sig < 1e-12) return(NULL)

    if(!is.null(h)){
        idh <- which(ll/n >= h & ll/n <= 1 - h)
        Z <- numeric(n)
        Z1 <- cumsum(x[1:max(idh)])
        ivn <- idh / n
        xup <- (Z1[idh] - (s * ivn)) / sqrt(n)
        xlo <- sqrt(n) * sig * sqrt(ivn * (1 - ivn))
        Z[idh] <- xup / xlo
        C <- n * Z^2
        idm <- which.max(C)
        Cm <- max(C)
        Vcritic <- matrix(
            c(NA, 0.010, 0.025, 0.050, 0.100,
            90.0, 9.209, 8.752, 8.312, 7.728,
            95.0, 10.788, 10.338, 9.885, 9.304,
            97.5, 12.331, 11.904, 11.409, 10.827,
            99.0, 14.364, 13.925, 13.421, 12.827,
            99.9, 19.278, 18.890, 18.377, 17.741),
            nrow = 6, byrow = TRUE)

        test <- Cm >= Vcritic[-1, which(Vcritic[1, ] == h)]
    }else{
        cus <- (cumsum(x) - (ll * s / n)) / sqrt(n)
        C <- abs(cus) / sig
        idm <- which.max(C)
        Cm <- max(C)
        Vcritic <- matrix(
            c(90.0, 1.224,
            95.0, 1.358,
            97.5, 1.480,
            99.0, 1.628,
            99.9, 1.949),
            ncol = 2, byrow = TRUE)

        test <- Cm >= Vcritic[, 2]
    }

    out <- NULL
    if(sum(test) > 0){
        mx.cl <- max(Vcritic[c(FALSE, test), 1])
        out <- list(index.change = ilen[idm]+1, stat.change = Cm, max.confL = mx.cl)
    }
    out
}

###################################################################
# CUSUM-type statistics with trends (Gallagher et al.,2013)

homog.CUSUM.Trend <- function(x, h = NULL){
    nl <- length(x)
    ilen <- 1:nl
    ina <- !is.na(x)
    ilen <- ilen[ina]
    x <- x[ina]
    n <- length(x)
    if(n < 3) return(NULL)
    ll <- 1:n
    epst <- as.vector(stats::residuals(stats::lm(x ~ ilen)))
    sig <- sqrt(sum(epst^2) / (n - 2))
    if(sig < 1e-12) return(NULL)

    if(!is.null(h)){
        idh <- which(ll/n >= h & ll/n <= 1 - h)
        Ts <- numeric(n)
        T1 <- cumsum(epst[1:max(idh)])
        ivn <- idh / n
        xlo <- sig * sqrt(n) * sqrt(ivn * (1 - ivn) * (1 - 3 * ivn * (1 - ivn)))
        Ts[idh] <- (-1) * T1[idh] / xlo
        Tc <- Ts^2
        idm <- which.max(Tc)
        Tm <- max(Tc)
        Vcritic <- matrix(
            c(NA, 0.010, 0.025, 0.050, 0.100,
            90.0, 10.341, 10.061, 9.790, 9.459,
            95.0, 11.956, 11.684, 11.415, 11.077,
            97.5, 13.532, 13.271, 12.989, 12.662,
            99.0, 15.540, 15.301, 15.011, 14.692,
            99.9, 20.600, 20.218, 20.114, 19.758),
            nrow = 6, byrow = TRUE)
        test <- Tm >= Vcritic[-1, which(Vcritic[1, ] == h)]
    }else{
        D <- cumsum(epst) / (sig * sqrt(n))
        Dc <- abs(D) / sig
        idm <- which.max(Dc)
        Dm <- max(Dc)
        Vcritic <- matrix(
            c(90.0, 0.836,
            95.0, 0.906,
            97.5, 0.970,
            99.0, 1.047,
            99.9, 1.222),
            ncol = 2, byrow = TRUE)
        test <- Dm >= Vcritic[, 2]
    }

    out <- NULL
    if(sum(test) > 0){
        mx.cl <- max(Vcritic[c(FALSE, test), 1])
        out <- list(index.change = ilen[idm] + 1, stat.change = Tm, max.confL = mx.cl)
    }
    out
}

###################################################################
# SNHT Critical Values

fSNHT.Vc <- function(n){
dat <- utils::read.table(text = 
"NA 90 92 94 95 97.5 99
10 4.964 5.197 5.473 5.637 6.188 6.769
12 5.288 5.554 5.876 6.068 6.729 7.459
14 5.54 5.831 6.187 6.402 7.152 8.001
16 5.749 6.059 6.441 6.674 7.492 8.44
18 5.922 6.248 6.652 6.899 7.775 8.807
20 6.07 6.41 6.83 7.089 8.013 9.113
22 6.2 6.551 6.988 7.257 8.22 9.38
24 6.315 6.675 7.123 7.4 8.4 9.609
26 6.417 6.785 7.246 7.529 8.558 9.812
28 6.509 6.884 7.353 7.643 8.697 9.993
30 6.592 6.973 7.451 7.747 8.825 10.153
32 6.669 7.056 7.541 7.841 8.941 10.3
34 6.741 7.132 7.625 7.93 9.05 10.434
36 6.803 7.201 7.699 8.009 9.143 10.552
38 6.864 7.263 7.768 8.081 9.23 10.663
40 6.921 7.324 7.835 8.151 9.317 10.771
42 6.972 7.38 7.894 8.214 9.39 10.865
44 7.022 7.433 7.951 8.273 9.463 10.957
46 7.071 7.484 8.007 8.331 9.53 11.04
48 7.112 7.529 8.054 8.382 9.592 11.116
50 7.154 7.573 8.103 8.432 9.653 11.193
52 7.194 7.616 8.149 8.48 9.711 11.259
54 7.229 7.654 8.19 8.524 9.76 11.324
56 7.264 7.69 8.23 8.566 9.81 11.382
58 7.299 7.727 8.268 8.606 9.859 11.446
60 7.333 7.764 8.308 8.647 9.906 11.498
62 7.363 7.796 8.343 8.683 9.947 11.548
64 7.392 7.827 8.375 8.717 9.985 11.599
66 7.421 7.857 8.408 8.752 10.026 11.648
68 7.449 7.886 8.439 8.784 10.067 11.692
70 7.475 7.913 8.467 8.814 10.099 11.737
72 7.499 7.938 8.496 8.844 10.134 11.776
74 7.525 7.965 8.523 8.873 10.171 11.822
76 7.547 7.989 8.548 8.898 10.2 11.858
78 7.57 8.013 8.575 8.926 10.23 11.895
80 7.591 8.035 8.599 8.951 10.259 11.928
82 7.613 8.059 8.623 8.976 10.29 11.966
84 7.634 8.079 8.647 9.001 10.315 11.995
86 7.655 8.102 8.67 9.026 10.347 12.033
88 7.673 8.121 8.691 9.047 10.37 12.059
90 7.692 8.14 8.71 9.067 10.394 12.089
92 7.711 8.16 8.732 9.09 10.417 12.12
94 7.73 8.181 8.752 9.11 10.447 12.153
96 7.745 8.196 8.77 9.127 10.465 12.175
98 7.762 8.214 8.788 9.147 10.484 12.196
100 7.778 8.231 8.807 9.167 10.507 12.228
105 7.819 8.273 8.851 9.214 10.562 12.291
110 7.856 8.312 8.892 9.255 10.608 12.343
115 7.891 8.35 8.931 9.296 10.656 12.401
120 7.921 8.38 8.963 9.33 10.694 12.446
125 7.952 8.413 8.999 9.365 10.735 12.488
130 7.983 8.446 9.032 9.4 10.772 12.538
135 8.01 8.474 9.063 9.431 10.808 12.579
140 8.038 8.501 9.092 9.462 10.845 12.621
145 8.063 8.529 9.12 9.49 10.877 12.66
150 8.086 8.554 9.147 9.519 10.906 12.694
155 8.111 8.578 9.172 9.543 10.933 12.725
160 8.133 8.601 9.195 9.569 10.966 12.759
165 8.155 8.625 9.222 9.596 10.992 12.793
170 8.174 8.643 9.241 9.615 11.016 12.82
175 8.195 8.666 9.265 9.641 11.046 12.851
180 8.214 8.685 9.283 9.658 11.062 12.872
185 8.233 8.706 9.307 9.683 11.089 12.904
190 8.252 8.725 9.325 9.701 11.11 12.93
195 8.268 8.741 9.343 9.72 11.132 12.956
200 8.286 8.761 9.364 9.741 11.156 12.982
225 8.361 8.838 9.446 9.826 11.247 13.083
250 8.429 8.908 9.516 9.898 11.329 13.175
275 8.489 8.97 9.581 9.966 11.399 13.248
300 8.54 9.022 9.635 10.02 11.46 13.326
325 8.587 9.07 9.685 10.071 11.517 13.389
350 8.633 9.117 9.732 10.118 11.565 13.44
375 8.67 9.157 9.775 10.161 11.613 13.494
400 8.706 9.193 9.814 10.202 11.654 13.542
425 8.738 9.224 9.844 10.234 11.692 13.58
450 8.771 9.26 9.882 10.272 11.73 13.623
475 8.798 9.288 9.912 10.302 11.761 13.655
500 8.828 9.317 9.939 10.33 11.795 13.69
525 8.854 9.344 9.967 10.36 11.827 13.73
550 8.878 9.369 9.995 10.386 11.854 13.751
575 8.901 9.391 10.016 10.408 11.878 13.782
600 8.923 9.414 10.04 10.431 11.904 13.813
650 8.963 9.455 10.083 10.476 11.949 13.856
700 9.001 9.493 10.119 10.511 11.986 13.904
750 9.033 9.524 10.152 10.547 12.026 13.947
800 9.063 9.557 10.187 10.58 12.059 13.975
850 9.093 9.587 10.216 10.612 12.096 14.023
900 9.119 9.614 10.244 10.64 12.12 14.041
950 9.143 9.638 10.269 10.665 12.149 14.07
1000 9.168 9.664 10.295 10.692 12.176 14.105
1100 9.211 9.708 10.339 10.736 12.22 14.15
1200 9.246 9.745 10.377 10.775 12.263 14.197
1300 9.283 9.781 10.415 10.812 12.304 14.235
1400 9.313 9.812 10.446 10.845 12.34 14.271
1500 9.347 9.846 10.481 10.88 12.374 14.312
1600 9.372 9.871 10.506 10.904 12.396 14.339
2000 9.464 9.965 10.603 11.002 12.5 14.443
2500 9.551 10.052 10.69 11.089 12.591 14.54
3000 9.618 10.121 10.76 11.161 12.664 14.619
3500 9.675 10.178 10.818 11.219 12.727 14.683
4000 9.727 10.229 10.869 11.271 12.779 14.734
4500 9.766 10.269 10.911 11.313 12.82 14.777
5000 9.803 10.307 10.948 11.349 12.859 14.817
7500 9.938 10.442 11.085 11.487 12.997 14.959
10000 10.031 10.537 11.18 11.584 13.095 15.063
15000 10.152 10.658 11.302 11.707 13.221 15.186
20000 10.236 10.743 11.388 11.791 13.305 15.271
50000 10.48 10.988 11.634 12.039 13.556 15.523",
sep = " ")

    len <- as.numeric(dat[-1, 1])
    vCI <- as.numeric(dat[1, -1])
    Vcritic <- as.matrix(dat[-1, -1])
    val <- sapply(1:6, function(j){
        f <- stats::splinefun(len, Vcritic[, j])
        f(n)
    })

    return(cbind(vCI, val))
}