## ----pkgs, echo=TRUE----------------------------------------------------------
if (!require("rspat")) remotes::install_github("rspatial/rspat")
if (!require("predicts")) install.packages("predicts")
if (!require("geodata")) install.packages("geodata")


## ----getData, echo=TRUE-------------------------------------------------------
library(terra)
library(rspat)
bf <- spat_data("bigfoot")
dim(bf)
head(bf)


## ----sup1---------------------------------------------------------------------
plot(bf[,1:2], cex=0.5, col="red")

library(geodata)
wrld <- geodata::world(path=".")
bnds <- wrld[wrld$NAME_0 %in% c("Canada", "Mexico", "United States"), ]
lines(bnds)


## ----sup3---------------------------------------------------------------------
wc <- geodata::worldclim_global("bio", res=10, ".")
plot(wc[[c(1, 12)]], nr=2)


## ----sup5x--------------------------------------------------------------------
bfc <- extract(wc, bf[,1:2])
head(bfc, 3)


## ----sup5---------------------------------------------------------------------
bfc <- bfc[,-1]


## ----sup7---------------------------------------------------------------------
plot(bfc[ ,"wc2.1_10m_bio_1"], bfc[, "wc2.1_10m_bio_12"], col="red",
        xlab="Annual mean temperature (°C)", ylab="Annual precipitation (mm)")


## -----------------------------------------------------------------------------
ext_bf <- ext(vect(bf[, 1:2])) + 1
ext_bf


## -----------------------------------------------------------------------------
set.seed(0)
window(wc) <- ext_bf
bg <- spatSample(wc, 5000, "random", na.rm=TRUE, xy=TRUE)
head(bg)


## -----------------------------------------------------------------------------
plot(bg[, c("x", "y")])


## -----------------------------------------------------------------------------
bg <- bg[, -c(1:2)]


## -----------------------------------------------------------------------------
plot(bg[,1], bg[,12], xlab="Annual mean temperature (°C)", 
  	      ylab="Annual precipitation (mm)", cex=.8)

points(bfc[,1], bfc[,12], col="red", cex=.6, pch="+")
legend("topleft", c("observed", "background"), col=c("red", "black"), pch=c("+", "o"), pt.cex=c(.6, .8))   


## -----------------------------------------------------------------------------
#eastern points
bfe <- bfc[bf[,1] > -102, ]
#western points
bfw <- bfc[bf[,1] <= -102, ]


## -----------------------------------------------------------------------------
dw <- rbind(cbind(pa=1, bfw), cbind(pa=0, bg))
de <- rbind(cbind(pa=1, bfe), cbind(pa=0, bg))

dw <- data.frame(dw)
de <- data.frame(na.omit(de))

dim(dw)
dim(de)


## ----sup10a-------------------------------------------------------------------
library(rpart)
cart <- rpart(pa~., data=dw)


## ----sup10ab------------------------------------------------------------------
printcp(cart)
plotcp(cart)


## ----sup10abc-----------------------------------------------------------------
cart <- rpart(pa~., data=dw, cp=0.02)


## ----sup10b-------------------------------------------------------------------
library(rpart.plot)
rpart.plot(cart, uniform=TRUE, main="Regression Tree")


## ----sup10bb, fig.width=6, fig.height=4---------------------------------------
x <- predict(wc, cart)
x <- mask(x, wc[[1]])
x <- round(x, 2)
plot(x, type="class", plg=list(x="bottomleft"))


## -----------------------------------------------------------------------------
set.seed(123)
i <- sample(nrow(dw), 0.2 * nrow(dw))
test <- dw[i,]
train <- dw[-i,]


## ----sup11--------------------------------------------------------------------
fpa <- as.factor(train[, 'pa'])


## ----sup12a-------------------------------------------------------------------
library(randomForest)
crf <- randomForest(train[, 2:ncol(train)], fpa)
crf


## ----sup12b-------------------------------------------------------------------
varImpPlot(crf)


## ----sup14a-------------------------------------------------------------------
trf <- tuneRF(train[, 2:ncol(train)], train[, "pa"])
trf
mt <- trf[which.min(trf[,2]), 1]
mt


## ----sup14b-------------------------------------------------------------------
rrf <- randomForest(train[, 2:ncol(train)], train[, "pa"], mtry=mt, ntree=250)
rrf
plot(rrf)


## ----sup17a, fig.width=6, fig.height=4----------------------------------------
rp <- predict(wc, rrf, na.rm=TRUE)
plot(rp)


## ----sup17b-------------------------------------------------------------------
library(predicts)
eva <- pa_evaluate(predict(rrf, test[test$pa==1, ]), predict(rrf, test[test$pa==0, ]))
eva


## ----sup18--------------------------------------------------------------------
plot(eva, "ROC")


## ----sup18b-------------------------------------------------------------------
par(mfrow=c(1,2))
plot(eva, "boxplot")
plot(eva, "density")


## ----sup19--------------------------------------------------------------------
tr <- eva@thresholds
tr
plot(rp > tr$max_spec_sens)


## ----sup20a, fig.width=6, fig.height=4----------------------------------------
rc <- predict(wc, crf, na.rm=TRUE)
plot(rc)


## ----sup20b, fig.width=6, fig.height=4----------------------------------------
rc2 <- predict(wc, crf, type="prob", na.rm=TRUE)
plot(rc2, 2)


## ----sup22--------------------------------------------------------------------
eva2 <- pa_evaluate(predict(rrf, de[de$pa==1, ]), predict(rrf, de[de$pa==0, ]))
eva2
par(mfrow=c(1,2))
plot(eva2, "ROC")
plot(eva2, "boxplot")


## ----sup24, fig.width=6, fig.height=4-----------------------------------------
plot(rc)
points(bf[,1:2], cex=.25)


## ----sup26, fig.width=6, fig.height=4-----------------------------------------
window(wc) <- NULL
pm <- predict(wc, rrf, na.rm=TRUE)
plot(pm)
lines(wrld)


## ----sup28, fig.width=6, fig.height=4-----------------------------------------
fut <- cmip6_world("CNRM-CM6-1", "585", "2061-2080", var="bio", res=10, path=".")
names(fut)
names(wc)
names(fut) <- names(wc)
window(fut) <- ext_bf
pfut <- predict(fut, rrf, na.rm=TRUE)
plot(pfut)