Nothing
thin_max <- function(x, cols, npoints){
inds <- vector(mode="numeric")
this.dist <- as.matrix(dist(x[,cols], upper=TRUE))
inds <- c(inds, as.integer(runif(1, 1, length(this.dist[,1]))))
inds <- c(inds, which.max(this.dist[,inds]))
while(length(inds) < npoints){
min.dists <- apply(this.dist[,inds], 1, min)
this.ind <- which.max(min.dists)
if(length(this.ind) > 1){
print("Breaking tie...")
this.ind <- this.ind[1]
}
inds <- c(inds, this.ind)
}
return(x[inds,])
}
#myRF <- function(x, p, q) {
# tf <- tuneRF(x[, 2:ncol(x)], x[, "pa"], plot=FALSE)
# mt <- tf[which.min(tf[,2]), 1]
# rf1 <- randomForest(x[, 2:ncol(x)], x[, "pa"], mtry=mt, ntree=250,
# na.action = na.omit)
# rp <- terra::predict(p, rf1, na.rm=TRUE)
# erf <- pa_evaluate(predict(rf1, q[q$pa==1, ]), predict(rf1, q[q$pa==0, ]))
# tr <- erf@thresholds
# rf_fnl <- rp > tr$max_spec_sens
# return(RF = rf_fnl)
#}
myGLM <- function(f, x, p, q) {
gm <- glm(f, family = gaussian(link = "identity"), data = x)
pg <- terra::predict(p, gm, type="response")
ge <- pa_evaluate(predict(gm, q[q$pa==1, ]), predict(gm, q[q$pa==0, ]))
gtr <- ge@thresholds
glm_fnl <- pg > gtr$max_spec_sens
return(GLM = glm_fnl)
}
myMAXENT <- function(p, ox, ot, xy) {
maxent_available <- FALSE
if (MaxEnt()) {
maxent_available <- TRUE
xm <- MaxEnt(p, ox)
mx <- terra::predict(xm, p)
xe <- pa_evaluate(xm, p = ot, a = xy, x = p)
mr <- xe@thresholds
MX <- mx > mr$max_spec_sens
return(MAXENT=MX)
}
}
#myCTA <- function(f, x, p, q) {
# cm <- rpart(f, data = x)
# pc <- terra::predict(p, cm)
# ce <- pa_evaluate(predict(cm, q[q$pa==1, ]), predict(cm, q[q$pa==0, ]))
# ctr <- ce@thresholds
# cta_fnl <- pc > ctr$max_spec_sens
# return(CTA = cta_fnl)
#}
#myGBM <- function(f, x, p, q) {
# gbm_mod <- gbm::gbm(f, data = x, interaction.depth = 3,
# n.minobsinnode = 1, cv.folds = 5)
#
# pred_gb <- terra::predict(p, gbm_mod, na.rm=TRUE)
# egb <- predicts::pa_evaluate(predict(gbm_mod, q[q$pa==1, ]),
# predict(gbm_mod, q[q$pa==0, ]))
# gb <- egb@thresholds
# gbm_fnl <- pred_gb > gb$max_spec_sens
# return(GBM = gbm_fnl)
#}
sampleBuffer <- function(p, size, width) {
if (max(distance(p)) > 7000000) {
h <- buffer(p, width = width)
} else {
h <- convHull(p)
}
spatSample(h, size = size, method = "random")
}
#' Species distribution models
#'
#' This function computes species distribution models using
#' two modelling algorithms: generalized linear models,
#' and maximum entropy (only if \code{rJava} is available).
#' Note: this is an experimental function, and may change in the future.
#'
#' @param x A dataframe containing the species occurrences
#' and geographic coordinates. Column 1 labeled as "species", column 2 "lon",
#' column 3 "lat".
#' @param predictors A \code{SpatRaster} to extract values from the
#' locations in x on which the models will be projected.
#' @param pol A vector polygon specifying the boundary to restrict the
#' prediction. If \code{NULL}, the extent of input points is used.
#' @param algorithm Character. The choice of algorithm to run the species
#' distribution model. Available algorithms include:
#' \itemize{
#' \item \dQuote{all}: Calls all available algorithms: GLM, and MAXENT.
#' \item \dQuote{GLM}: Calls only Generalized linear model.
#' \item \dQuote{MAXENT}: Calls only Maximum entropy.
#' }
#' @param size Minimum number of points required to successfully run
#' a species distribution model especially for species with few occurrences.
#' @param thin Whether to thin occurrences
#' @param thin.size The size of the thin occurrences.
#' @param width Width of buffer in meter if x is in longitude/latitude CRS.
#' @param mask logical. Should y be used to mask? Only used if pol is a SpatVector
#' @rdname sdm
#' @importFrom terra distance convHull spatSample vect ext window<- rast nlyr
#' @importFrom terra geom resample crop median deepcopy as.polygons predict
#' @importFrom terra extract
#' @importFrom predicts folds MaxEnt pa_evaluate
#' @importFrom stats glm median formula gaussian dist runif
#' @importFrom smoothr smooth
#' @return A list with the following objects:
#' \itemize{
#' \item \code{ensemble_raster} The ensembled raster that predicts
#' the potential species distribution based on the algorithms selected.
#' \item \code{data} The dataframe of occurrences used to implement the model.
#' \item \code{polygon} Map polygons of the predicted distributions
#' analogous to extent-of-occurrence range polygon.
#' \item \code{indiv_models} Raster layers for the separate models that
#' predict the potential species distribution.
#' }
#' @references
#' Zurell, D., Franklin, J., König, C., Bouchet, P.J., Dormann, C.F., Elith, J.,
#' Fandos, G., Feng, X., Guillera‐Arroita, G., Guisan, A., Lahoz‐Monfort, J.J.,
#' Leitão, P.J., Park, D.S., Peterson, A.T., Rapacciuolo, G., Schmatz, D.R.,
#' Schröder, B., Serra‐Diaz, J.M., Thuiller, W., Yates, K.L., Zimmermann, N.E.
#' and Merow, C. (2020), A standard protocol for reporting species distribution
#' models. \emph{Ecography}, \strong{43}: 1261-1277.
#' @examples
#' \donttest{
#' # get predictor variables
#' library(predicts)
#' f <- system.file("ex/bio.tif", package="predicts")
#' preds <- rast(f)
#' #plot(preds)
#'
#' # get species occurrences
#' b <- file.path(system.file(package="predicts"), "ex/bradypus.csv")
#' d <- read.csv(b)
#'
#' # fit ensemble model for four algorithms
#' m <- sdm(d, predictors = preds, algorithm = "all")
#' # plot(m$ensemble_raster)
#' # plot(m$polygon, add=TRUE)
#' }
#' @export
sdm <- function(x, predictors = NULL, pol = NULL, thin = TRUE, thin.size = 500,
algorithm = "all", size = 50, width = 50000, mask = FALSE) {
x <- .matchnames(x)
name.sp <- unique(x$species)
if (is.null(predictors)) {
stop("you need to specify RasterStack of environmental data")
}
x <- x[, -1]
x <- unique(na.omit(x))
if(thin == TRUE) {
x <- thin_max(x=x, cols = c("lat", "lon"), npoints = thin.size)
}
x$source <- "rw"
x <- vect(x, crs="+proj=longlat +datum=WGS84")
if (length(x) < size) {
h <- geom(x, df = TRUE)
h <- h[, 3:4]
h$source <- "rw"
v <- sampleBuffer(x, size = size-length(x), width = width)
v <- geom(v, df = TRUE)
v <- v[, 3:4]
v$source <- "rn"
x <- rbind(h, v)
} else {
x <- geom(x, df = TRUE)
x <- x[, 3:4]
x$source <- "rw"
}
occ_csv <- x
x <- x[, -3]
b <- extract(predictors, x[, 1:2])
b <- b[, -1]
if(!is.null(pol)) {
pol <- buffer(pol, width=width)
} else {
pol <- ext(vect(x, geom = c("x", "y"))) + 1
}
predictors <- terra::crop(predictors, pol, mask = mask)
bg <- spatSample(predictors, size = nrow(x)*100, method = "random",
na.rm = TRUE, xy = TRUE, warn = FALSE, ext=ext(pol))
xy <- bg[, 1:2]
bg <- bg[, -c(1:2)]
j <- data.frame(rbind(cbind(pa = 1, b), cbind(pa = 0, bg)))
j <- na.omit(j)
# for random forest
i <- sample(nrow(j), 0.25 * nrow(j))
test <- j[i,]
train <- j[-i,]
# for maxent
fold <- predicts::folds(x, k=5)
occtest <- x[fold == 1, ]
occtrain <- x[fold != 1, ]
# for glm
Preds <- names(train)[-1]
y <- names(train)[1]
Formula <- formula(paste(y," ~ ", paste(Preds, collapse=" + ")))
if (algorithm == "all") {
models <- c()
# 1. Random Forest
#tryCatch({
# RF1 <- FALSE
# rf_fnl <- myRF(x = train, p = predictors, q = test)
# RF1 <- TRUE
# if (RF1) {
# names(rf_fnl) <- "RF"
# models <- c(models, rf_fnl)
# }
#}, error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
# 2. GLM
tryCatch({
GLM1 <- FALSE
glm_fnl <- myGLM(f = Formula, x = train, p = predictors, q = test)
GLM1 <- TRUE
if (GLM1) {
names(glm_fnl) <- "GLM"
models <- c(models, glm_fnl)
}
}, error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
# 3. Maxent
tryCatch({
MX1 <- FALSE
MX <- myMAXENT(p=predictors, ox=occtrain, xy=xy, ot=occtest)
MX1 <- TRUE
if (MX1) {
names(MX) <- "MAXENT"
models <- c(models, MX)
}
}, error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
# 4. gbm
#tryCatch({
# GBM1 <- FALSE
# gbm_fnl <- myGBM(f = Formula, x = train, p = predictors, q = test)
# GBM1 <- TRUE
# if (GBM1) {
# names(gbm_fnl) <- "GBM"
# models <- c(models, gbm_fnl)
# }
#}, error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
# 5. CTA
#tryCatch({
# CTA1 <- FALSE
# cta_fnl <- myCTA(f = Formula, x = train, p = predictors, q = test)
# CTA1 <- TRUE
# if (CTA1) {
# names(cta_fnl) <- "CTA"
# models <- c(models, cta_fnl)
# }
#}, error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
models <- rast(models)
} else {
models <- switch(algorithm,
#RF = myRF(x = train, p = predictors, q = test),
GLM = myGLM(f=Formula, x=train, p=predictors, q=test),
MAXENT = myMAXENT(p=predictors, ox=occtrain, xy=xy,
ot=occtest),
#GBM = myGBM(f=Formula, x=train, p=predictors, q=test),
#CTA = myCTA(f=Formula, x=train, p=predictors, q=test)
)
names(models) <- algorithm
}
m <- models
if (nlyr(m) > 1) {
m <- terra::median(m)
}
spo <- m > 0.5
pol <- smooth(as.polygons(spo), method = "ksmooth")
names(pol)[1] <- "median"
pol <- pol[pol$median==1,]
pol$species <- name.sp
occ_csv$species <- name.sp
return(list(ensemble_raster = spo, data = occ_csv,
polygon = pol, indiv_models = models))
}
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.