[Biomod-commits] biomod2 / Maxent

[Sami Domisch]

Dear Damien and Wilfried,

first of all many thanks for providing the new biomod2 package! I played 
a bit with the code in the multi species manual, and tried to include 
MAXENT within the sfLapply-function in snowfall. However, snowfall 
returns an error, stating that it "cannot open the connection". After 
turning MAXENT off, the models run fine.

 > Error in checkForRemoteErrors(val) :
 >  9 nodes produced errors; first error: cannot open the connection

I guess that this error derives from the temporary MAXENT-folder created 
in setwd()? The single nodes try to access that folder which is deleted 
after each run?

Is it possible to place the MaxentTmpData-folder within each species' 
folder during the model run? And related to that, is it possible to save 
the modelling results (save.image) also in each species' folder? In 
other words, how can I access the single species folder within the 
lapply-wrapper?

Again, thanks a lot for this nice tool!
Best, Sami


#----------------------------------
library(biomod2)
### Load data
mySpeciesOcc <- read.table("mySpeciesOcc.txt", h=T)
coord = read.table("coordinates.txt", h=T)
mySpeciesOcc = cbind(coord, mySpeciesOcc)

### Fill 'NAs'
mySpeciesOcc[mySpeciesOcc==0] <- NA

require(raster)

myExpl = stack("present/bio1",
                "present/bio12",
                 "present/bio15")

sp.names <- names(mySpeciesOcc[3:12])

### Create function to run models within 'lapply'

MyBiomodSF <- function(sp.n){
   cat('/n',sp.n,'modeling...')
   ### definition of data for this run
   ## i.e keep only the column of our species
   myResp <- as.numeric(mySpeciesOcc[,sp.n])
   #myResp <- as.numeric(mySpeciesOcc[,sp.n])
   # get NAs id
   na.id <- which(is.na(myResp))
   # remove NAs to enforce PA sampling to be done on explanatory rasters
   myResp <- myResp[-na.id]
  # myRespCoord =mySpeciesOcc[1:2]
   myRespCoord = mySpeciesOcc[-na.id,c('long','lat')]
   myRespName = sp.n
   ### Initialisation
   myBiomodData <- BIOMOD_FormatingData(resp.var = myResp,
                                        expl.var = myExpl,
                                        resp.xy = myRespCoord,
                                        resp.name = myRespName,
                                        PA.nb.rep = 1,
                                        PA.nb.absences = 10000,
                                        PA.strategy = 'random')

   ### Options definition
   myBiomodOption <- BIOMOD_ModelingOptions()
   ### Modelling
   myBiomodModelOut <- BIOMOD_Modeling(
     myBiomodData,
     #models = 
c('SRE','CTA','RF','MARS','FDA','MAXENT','GLM','GAM','GBM','ANN'),
     models = c('SRE','RF','MAXENT'),
     models.options = myBiomodOption,
     NbRunEval=3,
     DataSplit=70,
     Yweights=NULL,
     VarImport=3,
     models.eval.meth = c('TSS','ROC'),
     SaveObj = TRUE,
     rescal.all.models = TRUE)

   #   ### Building ensemble-models

     myBiomodEM <- BIOMOD_EnsembleModeling(
       modeling.output = myBiomodModelOut,
       chosen.models = 'all',
       eval.metric = c('TSS'),
       eval.metric.quality.threshold = c(0.4), # e.g. 0.4
       prob.mean = T,
       prob.cv = T,
       prob.ci = T,
       prob.ci.alpha = 0.05,
       prob.median = T,
       committee.averaging = T,
       prob.mean.weight = T,
       prob.mean.weight.decay = 'proportional' )

   ### Do projections on current variable
   myBiomomodProj <- BIOMOD_Projection(
     modeling.output = myBiomodModelOut,
     new.env = myExpl,
     proj.name = 'current',
     selected.models = 'all',
     binary.meth= 'ROC',
     compress = 'xz',
     clamping.mask = F)

   ### Do ensemble-models projections on current varaiable
   myBiomodEF <- BIOMOD_EnsembleForecasting(
     projection.output = myBiomomodProj,
     EM.output = myBiomodEM,
     binary.meth = 'TSS',
     total.consensus = TRUE)

}

#----------------------------------


### Init snowfall
library(snowfall)
sfInit(parallel=TRUE, cpus=8, type="SOCK")

### Export packages
sfLibrary('biomod2', character.only=TRUE)
sfLibrary('rgdal', character.only=TRUE)

### Export data to nodes
sfExportAll()

### Check the items each node has loaded so far:
sfClusterEval(ls())

### Do the run
start <- Sys.time()
result <- sfClusterApplyLB(sp.names, MyBiomodSF)
Sys.time()-start

### Remove all variables from nodes
# sfRemoveAll(except=c("sp.names", "MyBiomodSF"))

### Stop snowfall
sfStop()
#----------------------------------