# Copyright 2015 Province of British Columbia
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and limitations under the License.
#' Perform Mann-Kendall trend test on many wells
#'
#' Uses the zyp package to calculate Mann-Kendall trend test on pre-whitened
#' data (to remove autocorrelation) for many wells, using one or both of two
#' pre-whitening methods, see zyp documentation.
#'
#' @param dataframe Dataframe containing: an ID column (specified in \code{by})
#' and a column of values
#' @param wells Vector of well numbers to test. Default NULL does all in
#' dataframe
#' @param byID The name of the ID column
#' @param col The name of the column with the GWL values
#' @param method "both" (default), "yuepilon", or "zhang"
#' @return A dataframe of results for all wells evaluated.
#' @export
gwl_zyp_test <- function(dataframe, wells = NULL,
byID = "Well_Num", col = "mean_GWL", method = "both") {
if (!byID %in% names(dataframe)) stop(byID, " is not a column in the data frame")
if (!col %in% names(dataframe)) stop(col, " is not a column in the data frame")
if (is.null(wells)) {
wells <- unique(dataframe[[byID]])
} else {
wells <- wells
}
# create an empty dataframe to store results
mk.results <- data.frame(ID = character(), test_type = character(),
lbound = numeric(), trend = numeric(), trendp = numeric(),
ubound = numeric(), tau = numeric(), sig = numeric(),
nruns = numeric(), autocor = numeric(),
valid_frac = numeric(), linear = numeric(),
intercept = numeric(), stringsAsFactors = FALSE)
for (well in wells) {
d <- dataframe[[col]][dataframe[[byID]] == well]
if (method == "both" | method == "yuepilon") {
zyp.yuepilon <- zyp::zyp.trend.vector(d, method = "yuepilon", conf.intervals = TRUE)
mk.results[nrow(mk.results) + 1, 1:2] <- c(well, "yuepilon")
mk.results[nrow(mk.results), 3:13] <- zyp.yuepilon
}
if (method == "both" | method == "zhang") {
zyp.zhang <- zyp::zyp.trend.vector(d, method = "zhang", conf.intervals = TRUE)
mk.results[nrow(mk.results) + 1, 1:2] <- c(well, "zhang")
mk.results[nrow(mk.results), 3:13] <- zyp.zhang
}
}
names(mk.results)[1] <- byID
mk.results
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.