# NOT RUN {
# based on examples in the dismo package
# get predictor variables
library(dismo)
predictor.files <- list.files(path=paste(system.file(package="dismo"), '/ex', sep=''),
pattern='grd', full.names=TRUE)
predictors <- stack(predictor.files)
# subset based on Variance Inflation Factors
predictors <- subset(predictors, subset=c("bio5", "bio6",
"bio16", "bio17"))
predictors
predictors@title <- "base"
# presence points
# presence points
presence_file <- paste(system.file(package="dismo"), '/ex/bradypus.csv', sep='')
pres <- read.table(presence_file, header=TRUE, sep=',')[,-1]
# choose background points
background <- randomPoints(predictors, n=1000, extf = 1.00)
# if desired, change working directory where subfolders of "models" and
# "ensembles" will be created
# raster layers will be saved in subfolders of /models and /ensembles:
getwd()
# first calibrate the ensemble
# calibration is done in two steps
# in step 1, a k-fold procedure is used to determine the weights
# in step 2, models are calibrated for all presence and background locations
# factor is not used as it is not certain whether correct levels will be used
# it may therefore be better to use dummy variables
# step 1: determine weights through 4-fold cross-validation
ensemble.calibrate.step1 <- ensemble.calibrate.weights(
x=predictors, p=pres, a=background, k=4,
SINK=TRUE, species.name="Bradypus",
MAXENT=1, MAXLIKE=0, GBM=1, GBMSTEP=0, RF=0, GLM=1, GLMSTEP=0,
GAM=1, GAMSTEP=0, MGCV=1, MGCVFIX=0,
EARTH=1, RPART=1, NNET=1, FDA=1, SVM=1, SVME=1, GLMNET=1,
BIOCLIM.O=0, BIOCLIM=1, DOMAIN=1, MAHAL=0, MAHAL01=1,
ENSEMBLE.tune=TRUE, PROBIT=TRUE,
ENSEMBLE.best=0, ENSEMBLE.exponent=c(1, 2, 3),
ENSEMBLE.min=c(0.65, 0.7),
Yweights="BIOMOD",
PLOTS=FALSE, formulae.defaults=TRUE)
# step 1 generated the weights for each algorithm
model.weights <- ensemble.calibrate.step1$output.weights
x.batch <- ensemble.calibrate.step1$x
p.batch <- ensemble.calibrate.step1$p
a.batch <- ensemble.calibrate.step1$a
MAXENT.a.batch <- ensemble.calibrate.step1$MAXENT.a
factors.batch <- ensemble.calibrate.step1$factors
dummy.vars.batch <- ensemble.calibrate.step1$dummy.vars
# step 2: calibrate models with all available presence locations
# weights determined in step 1 calculate ensemble in step 2
ensemble.calibrate.step2 <- ensemble.calibrate.models(
x=x.batch, p=p.batch, a=a.batch, MAXENT.a=MAXENT.a.batch,
factors=factors.batch, dummy.vars=dummy.vars.batch,
SINK=TRUE, species.name="Bradypus",
models.keep=TRUE,
input.weights=model.weights,
ENSEMBLE.tune=FALSE, PROBIT=TRUE,
Yweights="BIOMOD",
PLOTS=FALSE, formulae.defaults=TRUE)
# step 3: use previously calibrated models to create ensemble raster layers
# re-evaluate the created maps at presence and background locations
# (note that re-evaluation will be different due to truncation of raster layers
# as they wered saved as integer values ranged 0 to 1000)
ensemble.raster.results <- ensemble.raster(xn=predictors,
models.list=ensemble.calibrate.step2$models,
input.weights=model.weights,
SINK=TRUE, evaluate=TRUE,
RASTER.species.name="Bradypus", RASTER.stack.name="base")
# use the base map to check for changes in suitable habitat
# this type of analysis is typically done with different predictor layers
# (for example, predictor layers representing different possible future climates)
# In this example, changes from a previous model (ensemble.raster.results)
# are contrasted with a newly calibrated model (ensemble.raster.results2)
# step 1: 4-fold cross-validation
ensemble.calibrate2.step1 <- ensemble.calibrate.weights(
x=x.batch, p=p.batch, a=a.batch, MAXENT.a=MAXENT.a.batch,
factors=factors.batch, dummy.vars=dummy.vars.batch,
k=4,
SINK=TRUE, species.name="Bradypus",
MAXENT=1, MAXLIKE=0, GBM=1, GBMSTEP=0, RF=0, GLM=1, GLMSTEP=0,
GAM=1, GAMSTEP=0, MGCV=1, MGCVFIX=0,
EARTH=1, RPART=1, NNET=1, FDA=1, SVM=1, SVME=1, GLMNET=1,
BIOCLIM.O=0, BIOCLIM=1, DOMAIN=1, MAHAL=0, MAHAL01=1,
ENSEMBLE.tune=TRUE, PROBIT=TRUE,
ENSEMBLE.best=0, ENSEMBLE.exponent=c(1, 2, 3),
ENSEMBLE.min=c(0.65, 0.7),
Yweights="BIOMOD",
PLOTS=FALSE, formulae.defaults=TRUE)
model.weights2 <- ensemble.calibrate2.step1$output.weights
ensemble.calibrate2.step2 <- ensemble.calibrate.models(
x=x.batch, p=p.batch, a=a.batch, MAXENT.a=MAXENT.a.batch,
factors=factors.batch, dummy.vars=dummy.vars.batch,
SINK=TRUE, species.name="Bradypus",
models.keep=TRUE,
input.weights=model.weights2,
ENSEMBLE.tune=FALSE, PROBIT=TRUE,
Yweights="BIOMOD",
PLOTS=FALSE, formulae.defaults=TRUE)
ensemble.raster.results2 <- ensemble.raster(
xn=predictors,
models.list=ensemble.calibrate2.step2$models,
input.weights=model.weights2,
SINK=TRUE, evaluate=TRUE,
RASTER.species.name="Bradypus", RASTER.stack.name="recalibrated")
base.file <- paste(getwd(), "/ensembles/presence/Bradypus_base.grd", sep="")
other.file <- paste(getwd(), "/ensembles/presence/Bradypus_recalibrated.grd", sep="")
changed.habitat <- ensemble.habitat.change(base.map=base.file,
other.maps=c(other.file),
change.folder="ensembles/change")
change.file <- paste(getwd(), "/ensembles/change/Bradypus_recalibrated_presence.grd", sep="")
par.old <- graphics::par(no.readonly=T)
dev.new()
par(mfrow=c(2,2))
raster::plot(raster(base.file), breaks=c(-1, 0, 1), col=c("grey", "green"),
legend.shrink=0.8, main="base presence")
raster::plot(raster(other.file), breaks=c(-1, 0, 1), col=c("grey", "green"),
legend.shrink=0.8, main="other presence")
raster::plot(raster(change.file), breaks=c(-1, 0, 1, 10, 11),
col=c("grey", "blue", "red", "green"),
legend.shrink=0.8, main="habitat change", sub="11 remaining, 10 lost, 1 new")
graphics::par(par.old)
areas <- ensemble.area(raster(change.file))
areas
# }
Run the code above in your browser using DataLab