## gridUncertainty.R Subset two or more grids to the common members
##
## Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es)
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
#' @title Calculate the uncertainty of the climate change signal.
#' @description Calculate the uncertainty of the climate change signal.
#' @param historical A grid with annual historical data
#' @param anomaly A grid (with time dimension = 1).
#' @param method "simple" or "advanced".
#' @return A grid with three values (uncertainty categories); 0 = robust signal,
#' 1 = no signal or no change, 2 = conflicting signals).
#' @author M. Iturbide
#' @export
computeUncertainty <- function(historical, anomaly, method = "simple"){
method <- match.arg(method, c("simple", "advanced"))
uncer3 <- aggregateGrid(anomaly, aggr.mem = list(FUN = modelAgreement, th = 80))
if(method == "advanced"){
si <- signal(historical, anomaly)
uncer1 <- aggregateGrid(si, aggr.mem = list(FUN = signalAgreement, th = 66, condition = "GT"))
uncer2 <- aggregateGrid(si, aggr.mem = list(FUN = signalAgreement, th = 66, condition = "LT"))
uncer23 <- binaryGrid(gridArithmetics(uncer2, uncer3, operator = "+"), "GT", 0)
uncer.a.aux1 <- gridArithmetics(uncer1, 1, -1, operator = c("-", "*"))
uncer.a.aux2 <- gridArithmetics(uncer23, 1, -1, operator = c("-", "*"))
eval(parse(text = paste("uncer.a.aux2$Data[which(uncer.a.aux2$Data == 1)] <- 2")))
uncer.a.aux <- gridArithmetics(uncer.a.aux1, uncer.a.aux2, operator = c("+"))
eval(parse(text = paste("uncer.a.aux$Data[which(uncer.a.aux$Data > 1)] <- 2")))
uncer.a.aux
} else {
uncer3
}
}
#' @title Calculate signal agreement.
#' @description Calculate signal agreement.
#' @param x binary numeric
#' @param th threshold
#' @param condition character (choices: "GT", "GE", "LT", "LE")
#' @return binary numeric (0 = no agreement, 1 = agreement).
#' @author M. Iturbide
#' @export
signalAgreement <- function(x, th = 66, condition = "GT") {
condition <- match.arg(condition, choices = c("GT", "GE", "LT", "LE"))
ineq <- switch(condition,
"GT" = ">",
"GE" = ">=",
"LT" = "<",
"LE" = "<=")
eval(parse(text = paste("as.numeric((sum(x)/length(x)*100)", ineq, "th)")))
}
#' @title Calculate model agreement on the sign of change.
#' @description Calculate if signal exists in the anomalies relative to the historical variability.
#' @param x numeric of anomalies
#' @param th percentage value (Default is 80).
#' @return binary numeric (0 = no agreement, 1 = agreement).
#' @author M. Iturbide
#' @export
#'
modelAgreement <- function(x, th = 80) {
mp <- mean(x, na.rm = TRUE)
if (is.na(mp)) {
1
} else {
if (mp > 0) {
as.numeric(sum(as.numeric(x > 0), na.rm = TRUE) > as.integer(length(x) * th / 100))
} else if (mp < 0) {
as.numeric(sum(as.numeric(x < 0), na.rm = TRUE) > as.integer(length(x) * th / 100))
} else if (mp == 0) {
1
}
}
}
#' @title Calculate if signal exists in the anomalies.
#' @description Calculate if signal exists in the anomalies relative to the historical variability.
#' @param historical A grid with historical data to calculate the variability threshols (ideally annual data)
#' @param anomaly A grid (with time dimension = 1).
#' @return A binary grid (0 = no signal, 1 = signal).
#' @author M. Iturbide
#' @export
signal <- function(historical, anomaly) {
vth <- getVariabilityThreshold(historical)
anomaly$Data <- abs(anomaly$Data)
sig <- gridArithmetics(anomaly, vth, operator = "-")
binaryGrid(sig, condition = "GT", threshold = 0, values = c(0, 1))
}
#' @title Calculate the variability threshold.
#' @description Subset two or more grids to the common members.
#' @param grid A grid (ideally with annual data).
#' @return A grid.
#' @author M. Iturbide
#' @export
getVariabilityThreshold <- function(grid) {
climatology(grid, clim.fun = list(FUN = vth.aux))
}
#' @title Compute variability threshold
#' @description Compute variability threshold
#' @param x numeric
#' @param period value of the number of years to consider as a period (default is 20)
#' @return numeric. A single value of the variability threshold as defined in the IPCC AR6 WGI (Atlas chapter)
#' @author M. Iturbide
#' @export
vth.aux <- function(x, period = 20) {
(sqrt(2) * 1.645 * sd(x))/sqrt(period)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.