## ---- include=FALSE----------------------------------------------------------- library(knitr) library(rspatial) opts_chunk$set( fig.width = 5, fig.height = 5, fig.cap = '', collapse = TRUE ) ## ----getData, echo=TRUE------------------------------------------------------- if (!require("rspatial")) remotes::install_github('rspatial/rspatial') library(rspatial) bf <- sp_data('bigfoot') dim(bf) head(bf) ## ---- sup1-------------------------------------------------------------------- plot(bf[,1:2], cex=0.5, col='red') library(maptools) data(wrld_simpl) plot(wrld_simpl, add=TRUE) ## ---- sup3-------------------------------------------------------------------- library(raster) wc <- raster::getData('worldclim', res=10, var='bio') plot(wc[[c(1, 12)]], nr=2) ## ---- sup5-------------------------------------------------------------------- bfc <- extract(wc, bf[,1:2]) head(bfc) # Any missing values? i <- which(is.na(bfc[,1])) i plot(bf[,1:2], cex=0.5, col='blue') plot(wrld_simpl, add=TRUE) points(bf[i, ], pch=20, cex=3, col='red') ## ---- sup7-------------------------------------------------------------------- plot(bfc[ ,'bio1'] / 10, bfc[, 'bio12'], xlab='Annual mean temperature (C)', ylab='Annual precipitation (mm)') ## ----------------------------------------------------------------------------- library(dismo) # extent of all points e <- extent(SpatialPoints(bf[, 1:2])) e # 5000 random samples (excluding NA cells) from extent e set.seed(0) bg <- sampleRandom(wc, 5000, ext=e) dim(bg) head(bg) ## ----------------------------------------------------------------------------- d <- rbind(cbind(pa=1, bfc), cbind(pa=0, bg)) d <- data.frame(d) dim(d) ## ----------------------------------------------------------------------------- de <- d[bf[,1] > -102, ] dw <- d[bf[,1] <= -102, ] ## ---- sup10a------------------------------------------------------------------ library(rpart) cart <- rpart(pa~., data=dw) printcp(cart) plotcp(cart) ## ---- sup10b------------------------------------------------------------------ plot(cart, uniform=TRUE, main="Regression Tree") # text(cart, use.n=TRUE, all=TRUE, cex=.8) text(cart, cex=.8, digits=1) ## ---- sup11------------------------------------------------------------------- library(randomForest) # create a factor to indicated that we want classification fpa <- as.factor(dw[, 'pa']) ## ---- sup12a------------------------------------------------------------------ crf <- randomForest(dw[, 2:ncol(dw)], fpa) crf plot(crf) ## ---- sup12b------------------------------------------------------------------ varImpPlot(crf) ## ---- sup14a------------------------------------------------------------------ trf <- tuneRF(dw[, 2:ncol(dw)], dw[, 'pa']) trf mt <- trf[which.min(trf[,2]), 1] mt ## ---- sup14b------------------------------------------------------------------ rrf <- randomForest(dw[, 2:ncol(d)], dw[, 'pa'], mtry=mt) rrf plot(rrf) varImpPlot(rrf) ## ----------------------------------------------------------------------------- # Extent of the western points ew <- extent(SpatialPoints(bf[bf[,1] <= -102, 1:2])) ew ## ---- sup17a------------------------------------------------------------------ rp <- predict(wc, rrf, ext=ew) plot(rp) ## ---- sup17b------------------------------------------------------------------ eva <- evaluate(dw[dw$pa==1, ], dw[dw$pa==0, ], rrf) eva ## ---- sup18------------------------------------------------------------------- plot(eva, 'ROC') ## ---- sup19------------------------------------------------------------------- tr <- threshold(eva) tr plot(rp > tr[1, 'spec_sens']) ## ---- sup20a------------------------------------------------------------------ rc <- predict(wc, crf, ext=ew) plot(rc) ## ---- sup20b------------------------------------------------------------------ rc2 <- predict(wc, crf, ext=ew, type='prob', index=2) plot(rc2) ## ---- sup22------------------------------------------------------------------- de <- na.omit(de) eva2 <- evaluate(de[de$pa==1, ], de[de$pa==0, ], rrf) eva2 plot(eva2, 'ROC') ## ---- sup24------------------------------------------------------------------- eus <- extent(SpatialPoints(bf[, 1:2])) eus rcusa <- predict(wc, rrf, ext=eus) plot(rcusa) points(bf[,1:2], cex=.25) ## ---- sup26------------------------------------------------------------------- mex <- getData('GADM', country='MEX', level=1) pm <- predict(wc, rrf, ext=mex) pm <- mask(pm, mex) plot(pm) ## ---- sup28------------------------------------------------------------------- fut <- getData('CMIP5', res=10, var='bio', rcp=85, model='AC', year=70) names(fut) names(wc) names(fut) <- names(wc) futusa <- predict(fut, rrf, ext=eus, progress='window') plot(futusa)