-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmaxent.USA.R
102 lines (73 loc) · 4.96 KB
/
maxent.USA.R
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
# This code generates a host map prediction for Tree of Heaven across
# the lower 48 United States for use in PoPS modeling
require(BIEN)
require(dismo)
require(sf)
require(rJava)
require(rgdal)
require(parallel)
require(ENMeval)
require(sdmvspecies)
options(java.parameters = "-Xmx8000m")
usa <- readOGR('H:\\Shared drives\\APHIS Projects\\shared resources\\data\\usa_boundaries\\us_lower_48_states.shp')
#usa <- readOGR('C:\\Users\\bjselige\\Downloads\\us_lower_48_states.shp')
usa <- spTransform(usa, CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'))
toh.BIEN <- read.csv('H:\\Shared drives\\APHIS Projects\\PoPS\\Case Studies\\spotted_latternfly\\Host_mapping\\Ailanthus.BIEN.csv')
toh.EDDall <- read.csv('H:\\Shared drives\\APHIS Projects\\PoPS\\Case Studies\\spotted_latternfly\\Host_mapping\\Ailanthus.eddmaps.csv')
#toh.BIEN <- read.csv('C:\\Users\\bjselige\\Documents\\Tree_of_Heaven\\Ailanthus.BIEN.csv')
#toh.EDDall <- read.csv('C:\\Users\\bjselige\\Documents\\Tree_of_Heaven\\Ailanthus.eddmaps.csv')
toh.EDDall <- toh.EDDall[-which(is.na(toh.EDDall$Longitude)),]
toh.EDDpos <- toh.EDDall[which(toh.EDDall$OccStatus=="Positive"), c('Longitude', 'Latitude')]
toh.BIEN <- toh.BIEN[, c('Longitude', 'Latitude')]
toh.coords <- unique(rbind(toh.BIEN, toh.EDDpos))
toh.xy <- SpatialPoints(toh.coords, proj4string = CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 '))
toh.xy <- crop(toh.xy, usa)
biodir <- 'H:\\Shared drives\\APHIS Projects\\shared resources\\data\\worldclim1k\\US\\'
biovars <- stack(lapply(X=list.files(biodir), FUN=function(X){raster(paste(biodir, X, sep=''))}))
rails.d <- raster('H:\\Shared drives\\APHIS Projects\\shared resources\\data\\Rails_Roads\\Products_generated_from_Rails_Roads\\rails.distance.1km.tif')
rails.d <- resample(rails.d, biovars[[1]], method='bilinear'); rails.d <- rails.d + .5
roads.d <- raster('H:\\Shared drives\\APHIS Projects\\shared resources\\data\\Rails_Roads\\Products_generated_from_Rails_Roads\\roads.distance.tif')
roads.d <- resample(roads.d, biovars[[1]], method='bilinear')
#ronra.d <- overlay(roads.d, rails.d, fun=min)
set.seed(1991)
folds <- kfold(toh.xy, k=5)
train <- toh.xy[folds!=5,]
test <- toh.xy[folds==5,]
backg <- spsample(usa, n=length(test), type='regular')
backg.xy <- coordinates(backg); colnames(backg.xy) <- c('Longitude', 'Latitude')
# m.rora <- maxent(p=train, x=stack(roads.d, rails.d))
m.all <- maxent(p=train, x=biovars)
# m.all.ro <- maxent(p=train, x=stack(biovars, roads.d))
# m.all.rnr <- maxent(p=train, x=stack(biovars, ronra.d))
# m.all.rora <- maxent(p=train, x=stack(biovars, roads.d, rails.d))
# m.top2.ro <- maxent(p=train, x=stack(biovars[[1]], biovars[[19]], roads.d))
m.top3.ro <- maxent(p=train, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d))
# m.top3.rora <- maxent(p=train, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d, rails.d))
# m2.top3.ro <- maxent(p=toh.xy, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d))
# m2.top3.rora <- maxent(p=toh.xy, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d, rails.d))
# m.top4.ro <- maxent(p=train, x=stack(biovars[[1]], biovars[[10]], biovars[[14]], biovars[[19]], roads.d))
# pm.all <- predict(m.all, x=biovars)
# pm.all.ro <- predict(m.all.ro, stack(biovars, roads.d))
# pm.all.rora <- predict(m.all.rora, stack(biovars, roads.d, rails.d))
# pm.top2.ro <- predict(m.top2.ro, x=stack(biovars[[1]], biovars[[19]], roads.d))
pm.top3.ro <- predict(m.top3.ro, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d))
# pm.top3.rora <- predict(m.top3.rora, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d, rails.d))
# pm2.top3.ro <- predict(m2.top3.ro, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d))
# pm2.top3.rora <- predict(m2.top3.rora, x=stack(biovars[[1]], biovars[[14]], biovars[[19]], roads.d, rails.d))
# pm.top4.ro <- predict(m.top4.ro, x=stack(biovars[[1]], biovars[[10]], biovars[[14]], biovars[[19]], roads.d))
slf <- raster('C:\\Users\\bjselige\\Downloads\\slf2019_infested.tif')
slf <- projectRaster(slf, crs=CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 '))
hist(pm.top3.ro[slf$slf2019_infested>0])
# hist(pm.top3.rora[slf$slf2019_infested>0])
# hist(pm2.top3.ro[slf$slf2019_infested>0])
sum(pm.top3.ro[slf$slf2019_infested>0]<.2, na.rm=T)/sum(pm.top3.ro[slf$slf2019_infested>0]>0, na.rm=T)
# sum(pm.top3.rora[slf$slf2019_infested>0]<.2, na.rm=T)/sum(pm.top3.rora[slf$slf2019_infested>0]>0, na.rm=T)
# sum(pm2.top3.ro[slf$slf2019_infested>0]<.2, na.rm=T)/sum(pm.top3.ro[slf$slf2019_infested>0]>0, na.rm=T)
# sum(pm2.top3.rora[slf$slf2019_infested>0]<.2, na.rm=T)/sum(pm.top3.rora[slf$slf2019_infested>0]>0, na.rm=T)
final.out <- pm.top3.ro
final.out[which(values(final.out)<.2)] <- 0
final.out <- rescale(final.out)
final.out <- final.out*100
writeRaster(final.out, 'C:\\Users\\bjselige\\Desktop\\toh.USA.tif', overwrite=T)
toh.usa <- raster('C:\\Users\\bjselige\\Desktop\\toh.USA.tif')
toh.usa <- aggregate(toh.usa, 4)