-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path.Rhistory
165 lines (165 loc) · 7.15 KB
/
.Rhistory
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
q();
q();
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/', warning=FALSE, message=FALSE)
install.packages('knitr')
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/', warning=FALSE, message=FALSE)
cd
pwd
??circles
exit
q();
exit
devtools::install_github("rstudio/rmarkdown")
install.packages(c("ade4", "adehabitatLT", "adehabitatMA", "backports", "BH", "boot", "caTools", "checkmate", "CircStats", "classInt", "cluster", "curl", "data.table", "deldir", "devtools", "digest", "doParallel", "dtw", "e1071", "evaluate", "fda", "fields", "FNN", "foreach", "forecast", "foreign", "Formula", "gdalUtils", "geoR", "geosphere", "ggplot2", "git2r", "gmodels", "gstat", "gtools", "highr", "Hmisc", "htmlTable", "htmlwidgets", "httpuv", "hydroGOF", "hydroTSM", "igraph", "irlba", "iterators", "knitr", "lattice", "lazyeval", "LearnBayes", "lmtest", "magic", "mapdata", "mapproj", "maps", "maptools", "MASS", "Matrix", "mgcv", "munsell", "mvtnorm", "openssl", "pgirmess", "pkgconfig", "polyclip", "quantmod", "R.oo", "R.utils", "Rcpp", "RcppArmadillo", "reshape", "reshape2", "rgdal", "rgeos", "RgoogleMaps", "rJava", "rjson", "rlang", "rpart", "rprojroot", "sandwich", "scales", "shiny", "sourcetools", "sp", "spacetime", "spam", "spatstat", "spatstat.data", "spatstat.utils", "spBayes", "spData", "spdep", "stringi", "stringr", "survival", "tibble", "timeDate", "tinytex", "tkrplot", "tseries", "TTR", "viridis", "viridisLite", "withr", "xtable", "xts", "yaml", "zoo"))
?rmarkdown::render
pwd
getwd()
ls()
?render
install.packages("rmarkdown")
install.packages("rmarkdown")
install.packages("rmarkdown")
unlink('Git/hub/geog5330-fall2018/week2/Part2-RMarkdown_cache', recursive = TRUE)
library(rmarkdown)
pwd
getwd()
setwd(' /Users/guofeng/Git/hub/geog5330-fall2018/week5')
knitr::opts_chunk$set(fig.width=6, fig.height=4, fig.path='Figs/', warning=FALSE, message=FALSE, results='hide')
rm(list=ls())
library(sp)
library(spatstat)
library(dismo)
# Materm cluster proccess
cluster <- rMatClust(90, 0.05, 10)
plot(cluster)
# log-Gaussian Cox procces
lgcp <- rLGCP("exp", 3, var=0.2, scale=.1 )
plot(lgcp)
# inhomogeneous LGCP with Gaussian covariance function
m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=owin())
plot(m)
X <- rLGCP("gauss", m, var=0.15, scale =0.5)
plot(attr(X, "Lambda"))
points(X)
files <- list.files(path=paste(system.file(package="dismo"), '/ex' , sep= ''), pattern= 'grd', full.names=TRUE )
mask <- raster(files[1])
set.seed(1963)
bg <- randomPoints(mask, 500 )
# set up the plotting area for two maps
par(mfrow=c(1,2))
plot(!is.na(mask), legend=FALSE)
points(bg, cex=0.5)
# now we repeat the sampling, but limit
# the area of sampling using a spatial extent
e <- extent(-80, -53, -39, -22)
bg2 <- randomPoints(mask, 50, ext=e)
plot(!is.na(mask), legend=FALSE)
plot(e, add=TRUE, col= ' red ')
points(bg2, cex=0.5)
file <- paste(system.file(package="dismo"), '/ex/acaule.csv' , sep='')
ac <- read.csv(file)
coordinates(ac) <- ~lon+lat
projection(ac)<-CRS('+proj=longlat +datum=WGS84')
data(wrld_simpl)
plot(mask)
points(ac)
plot(wrld_simpl, add=TRUE, border='dark grey')
# Generate a cicule centered at each acaule
x <- circles(ac, d=50000, lonlat=TRUE)
pol <- polygons(x)
# Now sample 250 random 'pseudo' location around each acaule sample locations
samp1 <- spsample(pol, 250, type= 'random' , iter=25)
points(samp1)
data(wrld_simpl)
library(maptools)
data(wrld_simpl)
?maxent
xm
ge1=evaluate(testpres, testbackg, gm1)
# sample locations
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep= ',')
bradypus <- bradypus[,-1]
# environment conditions
files <- list.files(path=paste(system.file(package="dismo"), '/ex', sep= ''), pattern= 'grd' , full.names=TRUE )
predictors <- stack(files)
# first layer of the RasterStack
plot(predictors, 1)
#plot(predictors)
# note the "add=TRUE" argument with plot
plot(wrld_simpl, add=TRUE)
# with the points function, "add" is implicit
points(bradypus, col= ' blue ')
# 'biome' is a categorical data that needs special treatment in
# regression case, let's simply drop it from the predicator stacks.
pred_nf <- dropLayer(predictors, 'biome')
# Generate background data. To speed up processing, let's restrict the
# predictions to a more restricted area:
ext = extent(-90, -32, -33, 23)
backg <- randomPoints(pred_nf, n=1000, ext=ext)
colnames(backg) = c('lon', 'lat')
backg=as.data.frame(backg)
# Randomly seperate the bradypus observations into 5 groups and use
# four of them as training data and the rest one for validation
group <- kfold(bradypus, 5)
pres_train <- bradypus[group != 1, ]
pres_test <- bradypus[group == 1, ]
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
# Display the training and validation test data
r = raster(pred_nf, 1)
plot(!is.na(r), col=c(' white ' , ' light grey '), legend=FALSE)
plot(ext, add=TRUE, col= ' red ' , lwd=2)
points(backg_train$lon, backg_train$lat, col= ' yellow ')
points(backg_test$lon, backg_test$lat, col= ' red ')
points(pres_train$lon, pres_train$lat, pch= '+' , col= ' green ')
points(pres_test$lon, pres_test$lat, pch= '+' , col= 'blue')
# Combine sample locations and environment conditions into one data frame
train <- rbind(pres_train, backg_train)
pb_train <- c(rep(1, nrow(pres_train)), rep(0, nrow(backg_train)))
envtrain <- extract(predictors, train)
envtrain <- data.frame( cbind(pa=pb_train, envtrain) )
envtrain[,'biome'] = factor(envtrain[,'biome'], levels=1:14)
head(envtrain)
testpres <- data.frame( extract(predictors, pres_test) )
testbackg <- data.frame( extract(predictors, backg_test) )
testpres[ ,'biome'] = factor(testpres[ ,'biome'], levels=1:14)
testbackg[ ,'biome'] = factor(testbackg[ ,'biome'], levels=1:14)
#In R , GLM is implemented in the 'glm' function, and the link function and
#error distribution are specified with the 'family' argument. Examples are:
#family = binomial(link = "logit")
#family = gaussian(link = "identity")
#family = poisson(link = "log")
# logistic regression:
gm1 <- glm(pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
family = binomial(link = "logit"), data=envtrain)
summary(gm1)
gm2 <- glm(pa ~ bio1+bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
family = gaussian(link = "identity"), data=envtrain)
# Model evaluation
# ROC curve (reveiver operating characterstic curve) and related AUC
# (area under the curve) are commonly used metrics for classification
# performance. The ROC curve is based on confusion matrix. X is false
# positive rate (FP/(FP+TN))and Y is true positive rate (TP/(TP+FN)).
# Applying these two rate to different sets of classification
# threshold leads to the ROC curve
ge1=evaluate(testpres, testbackg, gm1)
ge2=evaluate(testpres, testbackg, gm2)
pg <- predict(predictors, gm2, ext=ext)
par(mfrow=c(1,2))
plot(pg, main='GLM/gaussian, raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
tr <- threshold(ge2, 'spec_sens')
plot(pg > tr, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)
ge1
jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep= '')
jar
pwd
getwd()
setwd('~/Git/')
setwd('hub/geog5330-fall2018/week6/')
dir()