# Simple example
x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
a <- x * 0.6
b <- (1 - x) * 0.6
c <- 1 - (a + b)
lons <- seq(0, 359.5, length = 20)
lats <- seq(-89.5, 89.5, length = 10)
if (FALSE) {
PlotMostLikelyQuantileMap(list(a, b, c), lons, lats,
toptitle = 'Most likely tercile map',
bar_titles = paste('% of belonging to', c('a', 'b', 'c')),
brks = 20, width = 10, height = 8)
}
# More complex example
n_lons <- 40
n_lats <- 20
n_timesteps <- 100
n_bins <- 4
# 1. Generation of sample data
lons <- seq(0, 359.5, length = n_lons)
lats <- seq(-89.5, 89.5, length = n_lats)
# This function builds a 3-D gaussian at a specified point in the map.
make_gaussian <- function(lon, sd_lon, lat, sd_lat) {
w <- outer(lons, lats, function(x, y) dnorm(x, lon, sd_lon) * dnorm(y, lat, sd_lat))
min_w <- min(w)
w <- w - min_w
w <- w / max(w)
w <- t(w)
names(dim(w)) <- c('lat', 'lon')
w
}
# This function generates random time series (with values ranging 1 to 5)
# according to 2 input weights.
gen_data <- function(w1, w2, n) {
r <- sample(1:5, n,
prob = c(.05, .9 * w1, .05, .05, .9 * w2),
replace = TRUE)
r <- r + runif(n, -0.5, 0.5)
dim(r) <- c(time = n)
r
}
# We build two 3-D gaussians.
w1 <- make_gaussian(120, 80, 20, 30)
w2 <- make_gaussian(260, 60, -10, 40)
# We generate sample data (with dimensions time, lat, lon) according
# to the generated gaussians
sample_data <- multiApply::Apply(list(w1, w2), NULL,
gen_data, n = n_timesteps)$output1
# 2. Binning sample data
prob_thresholds <- 1:n_bins / n_bins
prob_thresholds <- prob_thresholds[1:(n_bins - 1)]
thresholds <- quantile(sample_data, prob_thresholds)
binning <- function(x, thresholds) {
n_samples <- length(x)
n_bins <- length(thresholds) + 1
thresholds <- c(thresholds, max(x))
result <- 1:n_bins
lower_threshold <- min(x) - 1
for (i in 1:n_bins) {
result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples
lower_threshold <- thresholds[i]
}
dim(result) <- c(bin = n_bins)
result
}
bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1
# 3. Plotting most likely quantile/bin
if (FALSE) {
PlotMostLikelyQuantileMap(bins, lons, lats,
toptitle = 'Most likely quantile map',
bar_titles = paste('% of belonging to', letters[1:n_bins]),
mask = 1 - (w1 + w2 / max(c(w1, w2))),
brks = 20, width = 10, height = 8)
}
Run the code above in your browser using DataLab