forked from hguerra/biomass
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmaxent_example.R
More file actions
76 lines (55 loc) · 1.9 KB
/
maxent_example.R
File metadata and controls
76 lines (55 loc) · 1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
require("sp")
require("rgdal")
require("dismo")
# only run if the maxent.jar file is available, in the right folder
jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep='')
# checking if maxent can be run (normally not part of your script)
if (file.exists(jar) & require(rJava)) {
# get predictor variables
fnames <- list.files(path=paste(system.file(package="dismo"), '/ex', sep=''), pattern='grd', full.names=TRUE )
predictors <- stack(fnames)
#plot(predictors)
# file with presence points
occurence <- paste(system.file(package="dismo"), '/ex/bradypus.csv', sep='')
occ <- read.table(occurence, header=TRUE, sep=',')[,-1]
# witholding a 20% sample for testing
fold <- kfold(occ, k=5)
occtest <- occ[fold == 1, ]
occtrain <- occ[fold != 1, ]
# fit model, biome is a categorical variable
me <- maxent(predictors, occtrain, factors='biome')
# see the maxent results in a browser:
# me
# use "args"
# me2 <- maxent(predictors, occtrain, factors='biome', args=c("-J", "-P"))
# plot showing importance of each variable
plot(me)
# response curves
# response(me)
# predict to entire dataset
r <- predict(me, predictors)
# with some options:
# r <- predict(me, predictors, args=c("outputformat=raw"), progress='text',
# filename='maxent_prediction.grd')
plot(r)
points(occ)
#testing
# background data
bg <- randomPoints(predictors, 1000)
#simplest way to use 'evaluate'
e1 <- evaluate(me, p=occtest, a=bg, x=predictors)
# alternative 1
# extract values
pvtest <- data.frame(extract(predictors, occtest))
avtest <- data.frame(extract(predictors, bg))
e2 <- evaluate(me, p=pvtest, a=avtest)
# alternative 2
# predict to testing points
testp <- predict(me, pvtest)
head(testp)
testa <- predict(me, avtest)
e3 <- evaluate(p=testp, a=testa)
e3
threshold(e3)
plot(e3, 'ROC')
}