Review ENMeval
require(dismo)
## Loading required package: dismo
## Loading required package: raster
## Loading required package: sp
require(ENMeval)
## Loading required package: ENMeval
## Loading required package: rJava
## Loading required package: parallel
require(phyloclim)
## Loading required package: phyloclim
## Loading required package: ape
##
## Attaching package: 'ape'
## The following objects are masked from 'package:raster':
##
## rotate, zoom
require(sp)
require(rgdal)
## Loading required package: rgdal
## Warning: package 'rgdal' was built under R version 3.4.3
## rgdal: version: 1.2-20, (SVN revision 725)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
## Path to GDAL shared files: /usr/local/anaconda3/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-7
require(rgeos)
## Loading required package: rgeos
## Warning: package 'rgeos' was built under R version 3.4.3
## rgeos version: 0.3-26, (SVN revision 560)
## GEOS runtime version: 3.6.2-CAPI-1.10.2 0
## Linking to sp version: 1.2-7
## Polygon checking: TRUE
require(ENMTools)
## Loading required package: ENMTools
require(spocc)
## Loading required package: spocc
taxon = 'Quercus virginiana'
#Build a best model
#get occ
occ = occ(taxon, from = 'gbif', limit =900)
## Warning in C_valid_tz(tzone): System timezone name is unknown. Please set
## environment variable TZ.
occdf = occ2df(occ)
summary(occ)
## <source> gbif
## <time> 2018-05-17 16:28:16
## <found> 1353
## <returned> 900
## <query options>
## scientificName: Quercus virginiana
## limit: 900
##
## <source> bison
## <time>
## <found>
## <returned>
## <query options>
##
## <source> inat
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ebird
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ecoengine
## <time>
## <found>
## <returned>
## <query options>
##
## <source> antweb
## <time>
## <found>
## <returned>
## <query options>
##
## <source> vertnet
## <time>
## <found>
## <returned>
## <query options>
##
## <source> idigbio
## <time>
## <found>
## <returned>
## <query options>
##
## <source> obis
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ala
## <time>
## <found>
## <returned>
## <query options>
##
#get climate
# Env = stack("/data/spbio/climgrids/bio.gri")
Env = stack("~/bio.gri")
ext = extent(c(-125, -40, 25, 60))
Env = crop(Env, ext)
#clean
loc = cbind(occdf$longitude, occdf$latitude)
loc = loc[loc[,1]<= -40,]
loc = na.omit(loc)
df = data.frame(occdf)
extr = extract(Env, cbind(df$longitude, df$latitude))
df = df[!is.na(extr[,1]),]
#thin
require(spThin)
## Loading required package: spThin
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.1-4 (2018-04-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: fields
## Loading required package: maps
## See www.image.ucar.edu/~nychka/Fields for
## a vignette and other supplements.
## Loading required package: knitr
df = data.frame(occdf)
extr = extract(Env, cbind(df$longitude, df$latitude))
df = df[!is.na(extr[,1]),]
thin<-thin(loc.data = df,
lat.col = "latitude",
long.col = "longitude",
spec.col = "name",
thin.par = 10,
reps = 10,
locs.thinned.list.return = T,
write.files = T,
max.files = 2,
out.dir = paste(taxon, "/", sep=''),
write.log.file = T)
## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Thu May 17 16:28:35 2018
## lat.long.thin.count
## 232
## 10
## [1] "Maximum number of records after thinning: 232"
## [1] "Number of data.frames with max records: 10"
## [1] "Writing new *.csv files"
## Warning in thin(loc.data = df, lat.col = "latitude", long.col =
## "longitude", : Created new output directory: Quercus virginiana/
## [1] "Writing file: Quercus virginiana/thinned_data_thin1.csv"
## [1] "Writing file: Quercus virginiana/thinned_data_thin2.csv"
newthin = read.csv(paste(taxon, "/thinned_data_thin1.csv", sep=''))
ENMeval model testing
#ENMeval
thinres = ENMevaluate(occ=newthin[,2:3],
env = Env,
method='block',
parallel=T,
numCores=4,
bg.coords = bg2,
fc=c("L", "LQ", "H"),
RMvalues=seq(1,4,1),
rasterPreds=F)
## Doing evaluations using spatial blocks...
## Of 4 total cores using 4
## Running in parallel...
## Warning in calc.aicc(nparm, occ, predictive.maps): Can't calculate AICc
## when rasterPreds=F - returning NA's.
## ENMeval completed in 2 minutes 50.6 seconds.
print(thinres@results)
## settings features rm full.AUC Mean.AUC Var.AUC Mean.AUC.DIFF
## 1 L_1 L 1 0.9278 0.9203615 0.017103171 0.03477701
## 2 LQ_1 LQ 1 0.9435 0.9215963 0.008023755 0.03435251
## 3 H_1 H 1 0.9534 0.9212871 0.009490735 0.03925093
## 4 L_2 L 2 0.9252 0.9178986 0.018479871 0.03613958
## 5 LQ_2 LQ 2 0.9407 0.9220281 0.008606305 0.03346379
## 6 H_2 H 2 0.9481 0.9178433 0.010990158 0.04004436
## 7 L_3 L 3 0.9244 0.9183450 0.017373956 0.03483506
## 8 LQ_3 LQ 3 0.9386 0.9249125 0.008348839 0.02943635
## 9 H_3 H 3 0.9460 0.9152286 0.013418348 0.04062877
## 10 L_4 L 4 0.9239 0.9200470 0.015813067 0.03318606
## 11 LQ_4 LQ 4 0.9378 0.9278261 0.008344391 0.02850927
## 12 H_4 H 4 0.9444 0.9113981 0.016795819 0.04310065
## Var.AUC.DIFF Mean.OR10 Var.OR10 Mean.ORmin Var.ORmin AICc
## 1 0.014450729 0.1336207 0.05992370 0.02155172 0.0018579073 NA
## 2 0.009293970 0.1637931 0.05836306 0.01293103 0.0006688466 NA
## 3 0.010057706 0.1767241 0.04882580 0.03017241 0.0024524376 NA
## 4 0.015605278 0.1465517 0.06747919 0.02155172 0.0018579073 NA
## 5 0.009871021 0.1637931 0.05281411 0.01293103 0.0006688466 NA
## 6 0.011829264 0.1422414 0.05239298 0.01724138 0.0011890606 NA
## 7 0.014499008 0.1293103 0.05083234 0.02155172 0.0018579073 NA
## 8 0.009995038 0.1551724 0.05885850 0.01293103 0.0006688466 NA
## 9 0.014426900 0.1379310 0.05390408 0.03017241 0.0036414982 NA
## 10 0.013158813 0.1422414 0.05794193 0.02155172 0.0018579073 NA
## 11 0.009711300 0.1594828 0.08152497 0.00862069 0.0002972652 NA
## 12 0.017818580 0.1379310 0.05390408 0.03879310 0.0060196195 NA
## delta.AICc w.AIC nparam
## 1 NA NA 14
## 2 NA NA 19
## 3 NA NA 62
## 4 NA NA 9
## 5 NA NA 13
## 6 NA NA 33
## 7 NA NA 7
## 8 NA NA 12
## 9 NA NA 26
## 10 NA NA 9
## 11 NA NA 11
## 12 NA NA 24
eval.plot(thinres@results, "Mean.AUC", )
#predictbest
setsort = thinres@results[order(thinres@results[,'Mean.ORmin']),]
setsort2 = setsort[order(setsort[,'Mean.AUC'], decreasing=TRUE),]
top = setsort2[1,]
best.thin = which(as.character(thinres@results[,1]) == as.character(setsort2[1,1]))
pred.thin = predict(Env, thinres@models[[best.thin]])
plot(pred.thin, col=viridis::viridis(99))
#thresholding
ev.set <- evaluate(newthin[,2:3], thinres@bg.pts, thinres@models[[best.thin]], Env)
th1 = threshold(ev.set)
p1.nomit = pred.thin >= th1$no_omission
p1.equal = pred.thin >= th1$equal_sens_spec ##probably this one is better
Paleoclimate
lgm = stack('/data/spbio/climgrids/miroc/lgm_envirem.gri')
names(lgm) = names(Env) #Make names match
lgm = crop(lgm, extent(Env))
lgm.pred.thin = predict(lgm, thinres@models[[best.thin]])
lgm.binary = lgm.pred.thin>= th1$equal_sens_spec
plot(lgm.binary+p1.equal, col = c('black','red', 'green', 'blue'))
