Skip to content

Commit

Permalink
Corrected ExDet when n(analogue) < 2
Browse files Browse the repository at this point in the history
  • Loading branch information
@pjbouchet
pjbouchet committed Jun 1, 2021
1 parent 9745ac9 commit 35a6306
Show file tree
Hide file tree
Showing 6 changed files with 59 additions and 39 deletions.
47 changes: 29 additions & 18 deletions R/ExDet.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,11 +45,6 @@ ExDet <- function(ref, tg, xp){
rep(0, nrow(tg) * ncol(tg))),
dim = c(dim(tg), 3)),
c(1, 2), min)/(maxref - minref))

#---------------------------------------------
# Rename columns
#---------------------------------------------

names(nt1.df) <- xp

#---------------------------------------------
Expand All @@ -69,7 +64,6 @@ ExDet <- function(ref, tg, xp){
#---------------------------------------------

# Zero indicates within univariate range

univ.rge <- which(nt1 == 0)
mic_nt1[univ.rge] <- NA

Expand Down Expand Up @@ -127,7 +121,7 @@ ExDet <- function(ref, tg, xp){
cov.aa <- cov.combs %>% purrr::map(., ~apply(as.matrix(ref[,.]), 2, mean))
cov.bb <- cov.combs %>% purrr::map(., ~var(as.matrix(ref[,.])))

}else{
} else {

cov.aa <- cov.combs %>% purrr::map(., ~apply(as.matrix(ref[,.]), 2, mean, na.rm = TRUE))
cov.bb <- cov.combs %>% purrr::map(., ~var(as.matrix(ref[,.]), na.rm = TRUE))
Expand All @@ -138,18 +132,35 @@ ExDet <- function(ref, tg, xp){
# Calculate Mahalanobis distance for each combination of covariates
#---------------------------------------------

if(length(xp) == 1){
mah_nt2 <- purrr::pmap(list(cov.combs, cov.aa, cov.bb),
function(a, b, c) stats::mahalanobis(x = as.matrix(tg.univ[,a]),
center = b,
cov = c))

}else{mah_nt2 <- purrr::pmap(list(cov.combs, cov.aa, cov.bb),
function(a, b, c) stats::mahalanobis(x = as.matrix(tg.univ[,a]),
center = b,
cov = c))
# Need min of two analogue observations to calculate Mahalanobis distance
if(nrow(tg.univ) < 2){

warning("Only one prediction point within analogue conditions. Mahalanobis distances cannot be calculated.")
mah_nt2 <- vector(mode = "list", length = length(cov.combs))

} else {

mah_nt2 <- purrr::pmap(.l = list(cov.combs, cov.aa, cov.bb),
.f = function(a, b, c) stats::mahalanobis(x = as.matrix(tg.univ[,a]),
center = b,
cov = c))
}

# if(length(xp) == 1){
#
# mah_nt2 <- purrr::pmap(list(cov.combs, cov.aa, cov.bb),
# function(a, b, c) stats::mahalanobis(x = as.matrix(tg.univ[,a]),
# center = b,
# cov = c))
#
# } else {
#
# mah_nt2 <- purrr::pmap(list(cov.combs, cov.aa, cov.bb),
# function(a, b, c) stats::mahalanobis(x = as.matrix(tg.univ[,a]),
# center = b,
# cov = c))
# }


#---------------------------------------------
# Retrieve missing covariate names
Expand Down Expand Up @@ -180,7 +191,7 @@ ExDet <- function(ref, tg, xp){
#---------------------------------------------

results <- tibble::tibble(ExDet = nt1, mic_univariate = mic_nt1, mic_combinatorial = NA)
results$mic_combinatorial[univ.rge] <- mic_nt2
if(nrow(tg.univ) > 1) results$mic_combinatorial[univ.rge] <- mic_nt2

# Analog conditions have no MIC

Expand Down
4 changes: 2 additions & 2 deletions R/compute_nearby.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,10 +148,10 @@ compute_nearby <- function (samples,
counterfact <- whatif.opt(formula = NULL,
data = make_X(calibration_data = samples,
test_data = samples,
covariate.names),
var_name = covariate.names),
cfact = make_X(calibration_data = samples,
test_data = prediction.grid,
covariate.names),
var_name = covariate.names),
nearby = nearby,
no.partitions = no.partitions,
verbose = verbose)
Expand Down
26 changes: 16 additions & 10 deletions R/map_extrapolation.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,11 +138,11 @@ map_extrapolation <- function(map.type = NULL,
# Define coordinate systems
#---------------------------------------------

crs.ll <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
suppressWarnings(crs.ll <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

# EPSG:3857 (Web Mercator) is needed by leaflet for plotting

crs.webmerc <- sp::CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
suppressWarnings(crs.webmerc <- sp::CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"))

#---------------------------------------------
# Check which type of extrapolation occurred
Expand Down Expand Up @@ -254,7 +254,7 @@ map_extrapolation <- function(map.type = NULL,

pal.univariate <- leaflet::colorNumeric(c("#7f2704", "#fd8d3c", "#fee6ce"),
raster::getValues(projr$univariate),na.color = "transparent")

suppressWarnings(
exleaf <- exleaf %>%
leaflet::addRasterImage(map = .,
x = projr$univariate,
Expand All @@ -273,7 +273,8 @@ map_extrapolation <- function(map.type = NULL,
r.values = raster::getValues(projr$univariate),
decreasing = TRUE,
group = "Univariate",
title = "Univariate")}
title = "Univariate")
)}

if("combinatorial"%in%names(projr)){

Expand All @@ -285,6 +286,7 @@ map_extrapolation <- function(map.type = NULL,
pal.combinatorial <- leaflet::colorNumeric(c("#deebf7", "#6baed6", "#08519c"),
raster::getValues(projr$combinatorial), na.color = "transparent")}

suppressWarnings(
exleaf <- exleaf %>%
leaflet::addRasterImage(map = .,
x = projr$combinatorial,
Expand All @@ -303,14 +305,16 @@ map_extrapolation <- function(map.type = NULL,
decreasing = TRUE,
group="Combinatorial",
r.values = raster::getValues(projr$combinatorial),
title = "Combinatorial")}
title = "Combinatorial")
)}

if("analogue"%in%names(projr)){

pal.analogue <- leaflet::colorNumeric(c("#e5f5e0", "#74c476", "#00441b"),
raster::getValues(projr$analogue),
na.color = "transparent")

suppressWarnings(
exleaf <- exleaf %>%
leaflet::addRasterImage(map = .,
x = projr$analogue,
Expand All @@ -329,7 +333,8 @@ map_extrapolation <- function(map.type = NULL,
decreasing = TRUE,
group="Analogue",
r.values = raster::getValues(projr$analogue),
title = "Analogue")}
title = "Analogue")
)}

}else if(map.type == "mic"){

Expand All @@ -348,6 +353,7 @@ map_extrapolation <- function(map.type = NULL,

mapvars <- dplyr::left_join(x = micvars, y = allvars, by = c("mic" = "ID"))

suppressWarnings(
exleaf <- exleaf %>% leaflet::addRasterImage(map = .,
x = raster::as.factor(extrapolation.values$rasters$mic$all),
colors = pal,
Expand All @@ -367,7 +373,7 @@ map_extrapolation <- function(map.type = NULL,
decreasing = TRUE,
group = "MIC",
r.values = raster::getValues(extrapolation.values$rasters$mic$all),
title = "MIC")
title = "MIC"))


}else if(map.type == "nearby"){
Expand All @@ -379,7 +385,7 @@ map_extrapolation <- function(map.type = NULL,
pal.near <- leaflet::colorNumeric(pals::parula(100),
raster::getValues(gower.values),
na.color = "transparent")

suppressWarnings(
exleaf <- exleaf %>%
leaflet::addRasterImage(map = .,
colors = pal.near,
Expand All @@ -397,7 +403,7 @@ map_extrapolation <- function(map.type = NULL,
decreasing = TRUE,
group ="% nearby",
r.values = raster::getValues(gower.values),
title = "% nearby")
title = "% nearby"))


}
Expand Down Expand Up @@ -487,7 +493,7 @@ map_extrapolation <- function(map.type = NULL,
overlayGroups = lyr.controls,
options = layersControlOptions(collapsed = FALSE))

if(verbose) warning('map_extrapolation relies on the leaflet package, which is built around a Web Mercator projection (EPSG:3857), and therefore requires rasters to be reprojected for plotting. As a result, minor discrepancies may occur between the interactive maps shown in the viewer, and the underlying raw data. The latter can be accessed directly from the extrapolation.values or gower.values objects and visualised using alternative packages such as ggplot2.')
if(verbose) warning('map_extrapolation relies on the leaflet package, which is built around a Web Mercator projection (EPSG:3857), and therefore requires rasters to be reprojected for plotting. As a result, minor discrepancies may occur between the interactive maps shown in the viewer, and the underlying raw data. The latter can be accessed directly from extrapolation object returned by <compute_extrapolation> and visualised using alternative packages such as ggplot2.')

return(exleaf)
}
2 changes: 1 addition & 1 deletion R/print.extrapolation_results_summary.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#' @return invisibly returns the summary part of the object only, printing the results
#' @keywords internal

print.extrapolation_results_summary <- function(x, digits=2, ...){
print.extrapolation_results_summary <- function(x, digits = 2, ...){

class(x) <- class(x)[-1]

Expand Down
1 change: 0 additions & 1 deletion R/summarise_extrapolation.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,6 @@ summarise_extrapolation <- function(extrapolation.object,

res <- res[zeroes.names]


#---------------------------------------------
# Most influential covariates - by extrapolation type
#---------------------------------------------
Expand Down
18 changes: 11 additions & 7 deletions R/utils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,8 +175,8 @@ rescale_cov <- function (ynew, y) { return ((ynew - mean(y, na.rm = TRUE)) / (st
make_X <- function (calibration_data,
test_data,
var_name){
X <- sapply(var_name, function (k) { rescale_cov (ynew = test_data[, k],
y = calibration_data[, k])})
# Changed from: rescale_cov(ynew = test_data[, k], y = calibration_data[, k]) -- June 1st, 2021
X <- sapply(var_name, function (k) {rescale_cov(ynew = test_data[[k]], y = calibration_data[[k]])})
X <- as.data.frame (X)
names (X) <- var_name
return (X)}
Expand All @@ -187,7 +187,7 @@ make_X <- function (calibration_data,

proj_rasters <- function(ll, coordinate.system){

crs.webmerc <- sp::CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
suppressWarnings(crs.webmerc <- sp::CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"))

llr <- ll # Copy list

Expand Down Expand Up @@ -215,22 +215,26 @@ proj_rasters <- function(ll, coordinate.system){
llr$all <- NULL
llr <- purrr::discard(.x = llr, is.null)

suppressWarnings(
llr <- purrr::map(.x = llr, # Same extent as the full raster, allows correct alignment
.f = ~raster::projectRaster(from = .x,
to = ll$all,
method = 'ngb')) %>%
purrr::map(.x = ., # CRS used by leaflet
.f = ~raster::projectRaster(from = .,
crs = crs.webmerc,
method = 'ngb'))
method = 'ngb')))

llr.univariate.ind <- raster::cellFromXY(object = llr$univariate,
xy = sp::coordinates(univariate.xy))

llr$univariate[llr.univariate.ind[which(is.na(llr$univariate[llr.univariate.ind]))]] <- univariate.values[which(is.na(llr$univariate[llr.univariate.ind]))]

duplicate.cells <- rbind(raster::as.data.frame(llr$univariate, xy = TRUE),
raster::as.data.frame(llr$analogue, xy = TRUE)) %>%
r1 <- raster::as.data.frame(llr$univariate, xy = TRUE)
r2 <- raster::as.data.frame(llr$analogue, xy = TRUE)
names(r1) <- names(r2) <- c("x", "y", "ExDet")

duplicate.cells <- rbind(r1, r2) %>%
stats::na.omit(.) %>%
dplyr::select(., x, y) %>%
.[duplicated(.),]
Expand All @@ -256,7 +260,7 @@ if(!class(coordinate.system)=="CRS"){
error = function(e) return(NA))

if(is.na(coord.err)){stop('Unrecognised coordinate system')
}else{coordinate.system <- sp::CRS(coordinate.system)}
}else{supressWarnings(coordinate.system <- sp::CRS(coordinate.system))}
}

return(coordinate.system)
Expand Down

0 comments on commit 35a6306

Please sign in to comment.