Nothing
enirg.predict <-
function(enirg.results, qtegv.maps = NULL, qlegv.maps = NULL, load.map = FALSE,
method = "normal", prediction.name = "predicted") {
if (class(enirg.results) != "enirg")
stop("The function predict.enirg needs an object of class 'enirg'!")
egv.map.names <- c(qlegv.maps, qtegv.maps)
li.map.names <- paste(enirg.results$species, "_li_", colnames(enirg.results$co), sep="")
pred.map.name <- paste(enirg.results$species, "_", prediction.name, "_pred", sep = "")
hsm.map.name <- paste(enirg.results$species, "_", prediction.name, "_hsm", sep = "")
HSM <- list()
if (!is.null(egv.map.names)) {
if (length(egv.map.names) != length(enirg.results$egvs))
stop("Number of prediction maps does not match number of maps used to estimate the niche!")
cat("\n\nCalculating li prediction maps from new set of EGVs ...\n\n")
li.map.names <- paste(enirg.results$species, "_", prediction.name,"_li_Mar", sep="")
for(i in 1:enirg.results$nf) {
li.map.names <- c(li.map.names, paste(enirg.results$species, "_", prediction.name,"_li_Spec", i, sep=""))
}
li.calc <- c()
for (j in 1:ncol(enirg.results$co)) {
pre1 <- paste(egv.map.names, enirg.results$co[, j], sep = " * ")
calc <- paste(li.map.names[j], "= 0", sep = "")
for (i in 1:length(pre1)) calc <- paste(calc, pre1[i], sep = " + ")
li.calc <- c(li.calc, calc)
}
if (method == "normal") {
for(i in 1:length(li.calc)) execGRASS("r.mapcalc", expression = li.calc[i], flags = "overwrite")
}
if (method == "large") {
for(i in 1:length(li.calc)) system(paste("r.mapcalc --overwrite", li.calc[i], sep = ""))
}
}
HSM[["EGV_map_names"]] <- egv.map.names
HSM[["li_map_names"]] <- li.map.names
HSM[["HSM_map_name"]] <- hsm.map.name
cat("\n\nCalculating Mahalanobis distances ...\n\n")
tmp <- as.matrix(enirg.results$obs.li[, c(3:(enirg.results$nf + 3), 2)])
mat.presences <- matrix(nrow = sum(tmp[, enirg.results$nf + 2]), ncol = ncol(enirg.results$co))
colnames(mat.presences) <- colnames(enirg.results$co)
samples <- 1:nrow(tmp)
position <- 1
for (i in samples) {
variables <- c(tmp[i, 1:(enirg.results$nf + 1)])
rep <- tmp[i, enirg.results$nf + 2]
ocupa <- (position):(position + rep - 1)
for (j in ocupa) {
mat.presences[j, 1:(enirg.results$nf + 1)] <- variables
}
position <- position + rep
}
mean.centre <- apply(mat.presences, 2, median)
cov <- t(as.matrix(mat.presences)) %*% as.matrix(mat.presences)/nrow(mat.presences)
cova <- solve(cov)
cat(paste("\n\nCalculating the Habitat Suitability Map ", enirg.results$species, "_", prediction.name,
"_hsm ...\n\n", sep = ""))
f3 <- function(i) enirg.results$co[, i]/sqrt(crossprod(enirg.results$co[, i])/length(enirg.results$co[,
i]))
c1 <- matrix(unlist(lapply(1:(enirg.results$nf + 1), f3)), length(enirg.results$egvs))
li_number <- length(li.map.names)
c1 <- paste("(", li.map.names, "-", mean.centre, ")", sep = "")
c <- (paste(c1, "*", cova, sep = ""))
c3 <- NULL
for (i in 1:(enirg.results$nf + 1)) {
c3[i] <- ""
for (j in 1:enirg.results$nf) c3[i] <- paste(c3[i], c[j + (i - 1) * (enirg.results$nf +
1)], "+", sep = "")
j <- j + 1
c3[i] <- paste(c3[i], c[j + (i - 1) * (enirg.results$nf + 1)], sep = "")
}
calc <- paste("(", c3[1], ")", "*", c1[1], sep = "")
i <- 1
if (length(c3) > 2)
for (i in 2:(length(c3) - 1)) calc <- paste(calc, "+", "(", c3[i], ")",
"*", c1[i], sep = "")
calc <- paste(calc, "+", "(", c3[length(c3)], ")", "*", c1[i + 1], sep = "")
calc.pred <- paste(pred.map.name, " = ", calc, sep = "")
if (method == "normal") {
execGRASS("r.mapcalc", expression = calc.pred, flags = "overwrite", legacyExec = TRUE)
statistic <- execGRASS("r.univar", map = pred.map.name, flags = "g", intern = TRUE, legacyExec = TRUE)
map.max <- as.numeric(strsplit(statistic[agrep(statistic, pattern = "max")],
"=")[[1]][[2]])
calc.hsm <- paste(hsm.map.name, "=1-(", pred.map.name, "/", map.max, ")", sep = "")
execGRASS("r.mapcalc", expression = calc.hsm, flags = "overwrite", legacyExec = TRUE)
if (load.map == TRUE) {
HSM[[hsm.map.name]] <- raster(readRAST(hsm.map.name))
}
}
if (method == "large") {
calc.pred <- paste("r.mapcalc --overwrite '", calc.pred, "'", sep = "")
system(calc.pred)
statistic <- execGRASS("r.univar", map = pred.map.name, flags = "g", intern = TRUE, legacyExec = TRUE)
map.max <- as.numeric(strsplit(statistic[agrep(statistic, pattern = "max")],
"=")[[1]][[2]])
calc.hsm <- paste("r.mapcalc --overwrite '", hsm.map.name, " = 1-(", pred.map.name, "/", map.max + 0.01, ")'", sep = "")
system(calc.hsm)
}
cat(paste("'", hsm.map.name, "' HSM map was sucessfully created!\n\n", sep = ""))
cat("You can find HSM in your mapset in GRASS. You can see it through QGIS or using 'raster' library instead.\n\n")
presences_map <- SpatialPointsDataFrame(enirg.results$presences, data = data.frame(presences = enirg.results$presences[, 3]))
pred.points <- paste(enirg.results$species, "_", prediction.name, sep = "")
if(pred.points %in% list.maps()$vector) {
cat(paste("'", pred.points, "' map is already present in your mapset!\n"))
cat(paste("Removing '", pred.points, "' map ...\n\n"))
execGRASS("g.remove", type = "vector",
name = pred.points, flags = "f")
}
writeVECT(SDF = presences_map, vname = pred.points,
v.in.ogr_flags = c("overwrite", "o"))
execGRASS("v.db.addcolumn", map = pred.points, columns = "pred double precision")
execGRASS("v.what.rast", map = pred.points, raster = hsm.map.name, column = "pred")
vect <- readVECT(pred.points)
vect <- na.exclude(data.frame(cbind(vect@data$presences, vect@data$pred)))
names(vect) <- c("observed", "predicted")
HSM[["predictions"]] <- vect
hsm.report <- strsplit(execGRASS("r.report", map = hsm.map.name,
nsteps = 10, units = "c", intern = T)[16:25], split = " ")
hist.hsm <- c()
for (i in 1:10) {
hist.hsm <- c(hist.hsm, as.numeric(gsub("([0-9]+).*$", "\\1", hsm.report[[i]][length(hsm.report[[i]])])))
}
intervals <- seq(0, 1, by = 0.1)
hist.predicted <- hist(vect$predicted, breaks = intervals, plot = FALSE)$counts
HSM[["validation"]] <- cbind(intervals[-11], intervals[-1], hist.hsm, hist.predicted)
colnames(HSM[["validation"]])[1:2] <- c("int.inf", "int.sup")
class(HSM) <- c("hsm")
return(HSM)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.