Skip to content

Iturbide et al 2016 (submitted)

Jump to bottom
miturbide edited this page Nov 10, 2016 · 21 revisions

In this tutorial we are using the same presence data used in other sections of this wiki, this is distributions of Quercus robur phylogenies H11 and H5 in Europe (object Oak_phylo2). Regarding climatic variables, we are using object biostackENSEMBLES, containing a subset of the data used in the manuscript (E-OBS observational data and the regional climate projections from SMHI_ECHAM5). This worked example is fully reproducible, as the data is available with the mopa package.

data(Oak_phylo2)
presences <- Oak_phylo2
data(biostackENSEMBLES)

Object biostackENSEMBLES is a list with two RasterStacks (see section 2.1. Data preprocessing), the first contains variables for the reference period ($baseline) and the second for the future period 2071-2100 ($future_71_100_SMHI_ECHAM5).

spplot(biostackENSEMBLES$baseline)

The spatial sampling background is given by the environmental data used for modelling. Thus, we can extract a grid (class SpatialPoints) from any of the rasters inside biostackENSEMBLES$baseline with function background.

sp_grid <- background(biostackENSEMBLES$baseline$bio2)

In case we want to delimit the background to a smaller area (i.e. the bounding coordinates of the presence data distribution), function delimit can be applied, which returns polygon/s data of the spatial delimitation for each phylogeny and also the corresponging coordinates in a matrix or matrixes (del$bbs.grid in the following example). Bounding coordinates can be specified by hand, or by using function boundingCoords.

bc <- boundingCoords(presences)
bc
##$H11
##[1] -7.97 13.87 37.89 58.48

##$H5
##[1]  3.42 31.00 37.67 60.99

del <- delimit(bounding.coords = bc, grid = sp_grid)
str(del$bbs.grid)

In this case, we are using the whole study domain for both phylogenies, thus we need to prepare the bounding coordinates the following way:

bc <- boundingCoords(coordinates(sp_grid))
bc
##[1] -23.28691  75.00081  25.49938  71.51928
del <- delimit(bounding.coords = bc, grid = sp_grid)

#TS method

#step1
unsuitable.bg <- OCSVMprofiling(xy = presences,
                               varstack = biostackENSEMBLES$baseline,
                               bbs.grid = del$bbs.grid)

#step2
ext <- bgRadio(xy = presences, 
               background = unsuitable.bg$absence,
               start = 0.249, by = 0.249,
               unit = "decimal degrees")

pa_TS <- list()
for (i in 1:10){
  pa_TS[[i]] <- pseudoAbsences(xy = presences, background = ext,
                               exclusion.buffer = 0.249,
                               prevalence = 0.5, kmeans = FALSE) 
}

presaus.TS <- list()
for(i in 1:10){
presaus.TS[[i]] <- bindPresAbs(presences = presences, absences = pa_TS[[i]])
}
 presaus.RS <- list()
 for (i in 1:10){
   pa_RS <- pseudoAbsences(xy = presences, background = del$bbs.grid,
                               exclusion.buffer = 0.249,
                               prevalence = 0.5, kmeans = FALSE) 
   presaus.RS[[i]] <- bindPresAbs(presences, pa_RS)
 }

###Modelling

dirs <- paste("mydirectory/modTS/round",1:10, sep="")
sapply(1:10, function(x){dir.create(dirs[x])})
for(i in 1:10){
  modmarsTS <- allModeling(data = presaus.TS[[i]], varstack = biostackENSEMBLES$baseline,
                          
                          k = 10, algorithm = "mars", destdir = dirs[i],
                          
                          projection = CRS("+proj=longlat +init=epsg:4326"))
  
  modglmTS <- allModeling(data = presaus.TS[[i]], varstack = biostackENSEMBLES$baseline,
                          
                          k = 10, algorithm = "glm", destdir = dirs[i],
                          
                          projection = CRS("+proj=longlat +init=epsg:4326"))
}
auc_mars_TS <- list()
auc_glm_TS <- list()

for(i in 1:10){
  auc_mars_TS[[i]] <-loadTestValues(models = modmarsTS, test = "auc",sourcedir = dirs[i])
  auc_glm_TS[[i]] <-loadTestValues(models = modglmTS, test = "auc", sourcedir = dirs[i])
}
ind.mars.TS <- list()
ind.glm.TS <- list()

for(i in 1:10){
  ind.mars.TS[[i]] <- indextent(testmat = auc_mars_TS[[i]], diagrams = TRUE)
  ind.glm.TS[[i]] <- indextent(testmat = auc_glm_TS[[i]], diagrams = TRUE)
}
mars.models <- list()
glm.models <- list()

for(i in 1:10){
mars.models[[i]] <- loadDefinitiveModel(ind.mars.TS[[i]], slot = "allmod")
glm.models[[i]] <- loadDefinitiveModel(ind.glm.TS[[i]], slot = "allmod")
}

###RS method

dirsRS <- paste("mydirectory/modRS/round",1:10, sep="")
sapply(1:10, function(x){dir.create(dirs[x])})


#predictions
prob.mars <- list()
prob.glm <- list()
for(i in 1:10){
  prob.mars[[i]] <- predictAllmod(model = mars.models[[i]], environment = biostackENSEMBLES$future_71_100_SMHI_ECHAM5)
  prob.glm[[i]] <- predictAllmod(glm.models[[i]], biostackENSEMBLES$future_71_100_SMHI_ECHAM5)
}
prob.data <- list("mars" = prob.mars, "glm" = prob.glm)

datarange <- list()
for (i in 1:length(prob.data)){
  maps <- list()
  data <- prob.data[[i]]
  H11 <- list()
  H5 <- list()
    for(j in 1:10){
      H11[[j]] <- data[[j]]$H11
      H5[[j]] <- data[[j]]$H5
    }
    names(H11) <- paste("round", 1:10, sep="")
    names(H5) <- names(H11)
    H11 <- calc(stack(H11), max, na.rm=T) - calc(stack(H11), min, na.rm=T)
    H5 <- calc(stack(H5), max, na.rm=T) - calc(stack(H5), min, na.rm=T)
    H11@data@values[which(is.infinite(H11@data@values))] <-NA
    H5@data@values[which(is.infinite(H5@data@values))] <-NA
    st <- stack(H11, H5)
    names(st) <- c("H11", "H5")  
    datarange[[i]] <- st
}
names(datarange) <- names(prob.data)

spplot(datarange$mars, at = seq(0,1,.1), col.regions = rev(gray.colors(10, end = 1, start = 0, gamma = 1.2)), strip = T)


spplot(datarange$glm, at = seq(0,1,.1), col.regions = rev(gray.colors(10, end = 1, start = 0, gamma = 1.2)), strip = T)

par(mfrow = c(1,2))
gridboxes <- length(na.omit(datarange$mars$H11@data@values))

#phylogeny H11
mars <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$mars$H11@data@values)))
glm <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$glm$H11@data@values)))
plot(mars[,2], mars[,1],xlim= c(0,1), ylim = c(0,100), col = "red", typ = "l", ylab = "gridboxes (%)", xlab = "range", main = "H11")
lines(glm[,2], glm[,1], typ = "l")
segments(x0 = seq(0.2,0.8,0.2), y0 = 0, x1 = seq(0.2,0.8,0.2), y1 = 100, col = "grey", lty = 3)
segments(x0 = 0, y0 = seq(20,80,20), x1 = 1, y1 = seq(20,80,20), col = "grey", lty = 3)
mars.5 <- round(100 - mars[which(mars[,2]>0.5),][1,][1], digits = 1)
glm.5 <- round(100 - glm[which(glm[,2]>0.5),][1,][1], digits = 1)
mars.25 <- round(100 - mars[which(mars[,2]>0.25),][1,][1], digits = 1)
glm.25 <- round(100 - glm[which(glm[,2]>0.25),][1,][1], digits = 1)
segments(x0 = 0.5, y0 = 0, x1 = 0.5, y1 = 100, col = "blue")
segments(x0 = 0.25, y0 = 0, x1 = 0.25, y1 = 100, col = "green")
legend(x = 0.51, y = 45, legend = paste(c(mars.5, glm.5),"   ",c(mars.25, glm.25), sep=" "),text.col = "green", lty = 1, col = c("red", "black"), border = F, box.lwd = 0, cex = .8, bty = "n")
legend(x = 0.51, y = 45, legend = c(mars.5, glm.5), text.col = "blue", lty = 1, col = c("red","black"), border = F, box.lwd = 0, cex = .8, bty = "n")

#phylogeny H5  
mars <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$mars$H5@data@values)))
glm <- cbind(1:gridboxes *100/gridboxes , sort(na.omit(datarange$glm$H5@data@values)))
plot(mars[,2], mars[,1], xlim= c(0,1), ylim = c(0,100), col = "red", typ = "l", ylab = "gridboxes (%)", xlab = "range", main = "H5")
lines(glm[,2], glm[,1], typ = "l")
segments(x0 = seq(0.2,0.8,0.2), y0 = 0, x1 = seq(0.2,0.8,0.2), y1 = 100, col = "grey", lty = 3)
segments(x0 = 0, y0 = seq(20,80,20), x1 = 1, y1 = seq(20,80,20), col = "grey", lty = 3)
mars.5 <- round(100 - mars[which(mars[,2]>0.5),][1,][1], digits = 1)
glm.5 <- round(100 - glm[which(glm[,2]>0.5),][1,][1], digits = 1)
mars.25 <- round(100 - mars[which(mars[,2]>0.25),][1,][1], digits = 1)
glm.25 <- round(100 - glm[which(glm[,2]>0.25),][1,][1], digits = 1)
segments(x0 = 0.5, y0 = 0, x1 = 0.5, y1 = 100, col = "blue")
segments(x0 = 0.25, y0 = 0, x1 = 0.25, y1 = 100, col = "green")
legend(x = 0.51, y = 45, legend = paste(c(mars.5, glm.5),"   ",c(mars.25, glm.25), sep=" "),text.col = "green", lty = 1, col = c("red", "black"), border = F, box.lwd = 0, cex = .8, bty = "n")
legend(x = 0.51, y = 45, legend = c(mars.5, glm.5), text.col = "blue", lty = 1, col = c("red","black"), border = F, box.lwd = 0, cex = .8, bty = "n")

legend(x = 0.5, y = 95, legend = c("mars", "glm"), lty = 1, col = c("red","black"), bty = "n",
       border = F, box.lwd = 0, cex = 1)


mars.H11 <- list()
glm.H11 <- list()
for(i in 1:10){
  mars.H11[[i]] <- prob.mars[[i]]$H11
  glm.H11[[i]] <- prob.glm[[i]]$H11
}
names(mars.H11) <- as.character(1:10)
names(glm.H11) <- names(mars.H11)
mars.H11 <- stack(mars.H11)
glm.H11 <- stack(glm.H11)

ov.mars.H11 <- nicheOver(stack = mars.H11, metric = "D")
library(lattice)
levelplot(ov.mars.H11)

ind <- which(ov.mars.H11 == min(ov.mars.H11, na.rm = T), arr.ind = T)
extremes <- stack(list(prob.mars[[ind[1]]]$H11, prob.mars[[ind[2]]]$H11))
spplot(extremes, strip = F)

center

ov.glm.H11 <- nicheOver(stack = glm.H11, metric = "D")
ind <- which(ov.glm.H11 == min(ov.glm.H11, na.rm = T), arr.ind = T)
extremes <- stack(list(prob.glm[[ind[1]]]$H11, prob.glm[[ind[2]]]$H11))
spplot(extremes, strip = F)

mopa - Santander MetGroup (Univ. Cantabria - CSIC)

mopa v0.4.0

  1. Data pre-processing
  2. Pseudo absence data generation and modelling

mopa v1.0.0

  1. Data formats and pre-processing
  2. Pseudo-absence data generation
  3. Model fitting and prediction
  4. Further options to assess model performace
  5. Analysis of SDM projections
Clone this wiki locally