date()
## [1] "Mon Oct 08 15:12:50 2018"
setwd("~/LabWork/Patuxent/USDA")
Cleaning eBird data using custom “CleanBird” function. Selecting mallard observations for years 2013-2017.
CleanBird("C:/Users/humph173/Documents/LabWork/Patuxent/eBird/ebd_mallar_relFeb-2018/ebd_mallar_relFeb-2018.txt", "MALL", 2000, 2017)
Getting state and county boundaries as well as the snap raster. These county boundaries don’t match the snap raster perfectly…
#Snap raster
SC.r = raster("C:/Users/humph173/Documents/LabWork/Patuxent/Jeff_03092018/INLA_Project/INLA_Project/SouthCarolina_USDA_Layers/MasterSnapRaster1.tif")
proj4string(SC.r)
## [1] "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
#Change to km for modeling
nProj = "+proj=utm +zone=17 +datum=WGS84 +units=km +no_defs +ellps=WGS84 +towgs84=0,0,0"
SC.r = projectRaster(SC.r, crs=CRS(nProj))
#States
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)
SC.st = subset(States, ID == 39)
#Counties
Counties = map("county",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(Counties$names, ":"),
function(x) x[1])
Counties = map2SpatialPolygons(Counties, IDs = IDs,
proj4string = CRS(LL84))
pid = sapply(slot(Counties, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(Counties),
row.names = pid)
Counties = SpatialPolygonsDataFrame(Counties, p.df)
#Select those for SC
SC.st = spTransform(subset(States, ID == 39), nProj)
SC.cnt = crop(spTransform(Counties, nProj), SC.st)
#Dense point grid for later prediction
SC.pnts = rasterToPoints(SC.r, sp=T)
SC.pnts@data = SC.pnts@data %>%
mutate(Long = SC.pnts@coords[,1],
Lat = SC.pnts@coords[,2]) %>%
select(-MasterSnapRaster1)
MALL.pnt = SpatialPointsDataFrame(MALL[,2:3], MALL)
proj4string(MALL.pnt) = LL84
MALL.pnt = spTransform(MALL.pnt, nProj)
MALL.pnt$SC = is.na(over(MALL.pnt, SC.st))[,1]
MALL.pnt = subset(MALL.pnt, SC == FALSE)
F.plt = levelplot(SC.r,
margin = FALSE,
xlab =NULL,
ylab = NULL,
maxpixels = 1e5,
col.regions = "tan",
colorkey = NULL,
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.5)) +
latticeExtra::layer(sp.polygons(SC.cnt, col = "darkgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(SC.st, col = "darkgray", lwd = 1)) +
latticeExtra::layer(sp.polygons(MALL.pnt, col = "red", cex = 1, pch=19))
F.plt
Defining study area.
State.b1 = gBuffer(SC.st, width = 10)
bdry = inla.sp2segment(State.b1)
MaxEdge = 5
mesh2D = inla.mesh.2d(boundary = bdry,
cutoff = 1,
max.edge = MaxEdge,
min.angle = 22)
mesh2D$n #Count of verticies
## [1] 9502
plot(mesh2D, rgl = F, main = " ")
dd = as.data.frame(cbind(mesh2D$loc[,1],
mesh2D$loc[,2]))
names(dd) = c("Long", "Lat")
dd$OBS = 0
dd$Season = "none"
dd$Spp = "node"
Obs = MALL.pnt@data %>%
mutate(OBS =1,
Spp = "Bird",
Long = MALL.pnt@coords[,1],
Lat = MALL.pnt@coords[,2]) %>%
select(Long, Lat, OBS, Season, Spp)
DG = SC.pnts@data %>%
mutate(OBS = 0,
Spp = "Grid",
Season = "none")
Mod.pnts = rbind(Obs, dd, DG)
Mod.pnts = SpatialPointsDataFrame(Mod.pnts[,c("Long","Lat")], Mod.pnts)
proj4string(Mod.pnts) = proj4string(MALL.pnt)
SetRast = raster("C:/Users/humph173/Documents/LandFire/US_140BPS_20180618/grid/us_140bps")
LF.files = list.files(path="C:/Users/humph173/Documents/LandFire",
pattern="*.ovr", full.names=T, recursive=T)
for(i in 1:length(LF.files)){
tmp.lf = raster(LF.files[i])
extent(tmp.lf) = extent(SetRast)
crs(tmp.lf)= CRS(proj4string(SetRast))
Tmp.vals = extract(tmp.lf,
spTransform(Mod.pnts,
CRS(proj4string(tmp.lf))),
method="simple")
if(i == 1){Names.lvl = names(tmp.lf)
hld.tmp = Tmp.vals}
else{Names.lvl = c(Names.lvl, names(tmp.lf))
hld.tmp = cbind(hld.tmp, Tmp.vals)}
}
Mod.cov.df = as.data.frame(hld.tmp)
names(Mod.cov.df) = Names.lvl
http://www.worldclim.org/bioclim
BioFiles = list.files(path="C:/Users/humph173/Documents/Michigan_State/Marten/Marten_SDM/Env/wc2.0_30s_bio",
pattern="*.tif", full.names=T, recursive=FALSE)
Bio.stk = stack(BioFiles)
names(Bio.stk) = paste("bio_", 1:19, sep="")
BioM = as.data.frame(
extract(Bio.stk,
spTransform(Mod.pnts,
CRS(proj4string(Bio.stk))),
method = "simple"))
Mod.cov.df = cbind(Mod.cov.df, BioM)
EnvFiles = list.files(path="C:/Users/humph173/Documents/ENVIREM",
pattern="*.tif", full.names=T, recursive=FALSE)
ENV.stk = stack(EnvFiles)
names(ENV.stk) = paste("env_", 1:18, sep="")
EnvM = as.data.frame(
extract(ENV.stk,
spTransform(Mod.pnts,
CRS(proj4string(ENV.stk))),
method = "simple"))
Mod.cov.df = cbind(Mod.cov.df, EnvM)
Cov.scale = Mod.cov.df
for(i in 1:ncol(Cov.scale)){
V.mean = mean(Cov.scale[,i], na.rm=T)
Cov.scale[,i] = ifelse(is.na(Cov.scale[,i]), V.mean, Cov.scale[,i])
Cov.scale[,i] = as.numeric(scale(Cov.scale[,i], scale = T, center = T))
}
names(Cov.scale) = paste("s_", names(Mod.cov.df), sep="")
Mod.cov.df = cbind(Mod.cov.df, Cov.scale)
Mod.pnts@data = cbind(Mod.pnts@data, Mod.cov.df)
#Optimize
dapc.e = dapc(Cov.scale,
Mod.pnts$OBS,
n.da=5, n.pca=50)
temp = a.score(dapc.e)
temp$pop.score
temp$mean
optim.a.score(dapc.e, plot=FALSE)
###Conduct DAPC
set.seed(1976)
dapc.e = dapc(Cov.scale,
Mod.pnts$OBS,
n.da=5, n.pca=49)
summary(dapc.e)
Mod.pnts$IndF = as.numeric(dapc.e$ind.coord)
dapc.er = rasterize(Mod.pnts,
SC.r,
"IndF",
background = NA)
http://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev10
Pop = raster("C:/Users/humph173/Documents/Michigan_State/Marten/Marten_SDM/Env/PopDensity.tif")/100
Pop[is.na(Pop)] = mean(values(Pop), na.rm=T)
Mod.pnts$Pop = extract(Pop, spTransform(Mod.pnts, CRS(proj4string(Pop))), method = "simple")
Mod.pnts$Pop[is.na(Mod.pnts$Pop)] = mean(Mod.pnts$Pop, na.rm=T)
Mod.fit = subset(Mod.pnts, Spp == "Bird" | Spp == "node")
MALL.wint = subset(Mod.fit, Season == "Winter" | Season == "none")
MALL.spr = subset(Mod.fit, Season == "Spring" | Season == "none")
MALL.brd = subset(Mod.fit, Season == "Breeding" | Season == "none")
MALL.fall = subset(Mod.fit, Season == "Fall" | Season == "none")
Grd.pnts = subset(Mod.pnts, Spp == "Grid")
FE.df = MALL.wint@data
locs = cbind(FE.df$Long, FE.df$Lat)
A.base = inla.spde.make.A(mesh2D,
alpha = 2,
loc=locs)
A.base2 = A.base
spde1 = inla.spde2.pcmatern(mesh2D, alpha = 2,
prior.range=c(1500, 0.01),
prior.sigma=c(1, 0.01),
constr = TRUE)
field1 = inla.spde.make.index("field1", spde1$n.spde)
field1c = inla.spde.make.index("field1c", spde1$n.spde)
FE.lst = list(c(field1,
list(intercept1 = 1)),
list(XX = FE.df[,"Long"]))
pp.stk = inla.stack(data = list(Y = cbind(FE.df$OBS, NA),
e = rep(0, dim(FE.df)[1])),
A = list(A.base2, 1),
effects = FE.lst,
tag = "pp.0")
#2nd Tier
FE.lst2 = list(c(field1c,
list(intercept2 = 1)),
list(IndF = FE.df[,"IndF"],
Pop = FE.df[,"Pop"]))
FE.df$Exp = ifelse(FE.df$OBS==1, 0, 1)
Base.stk = inla.stack(data = list(PA = FE.df$OBS,
e = FE.df$Exp),
A = list(A.base, 1),
effects = FE.lst,
tag = "base.0")
resp.stk = inla.stack(data = list(Y = cbind(NA, FE.df$OBS),
e = FE.df$Exp),
A = list(A.base, 1),
effects = FE.lst2,
tag = "resp.0")
Joint.stk = inla.stack(pp.stk, resp.stk)
Frm1 = Y ~ -1 + intercept1 +
intercept2 +
f(field1,
model=spde1) +
f(field1c,
copy = "field1",
fixed = FALSE) +
IndF + Pop
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(2.9424376, 5.8718553, 0.9807331, 6.8721504)
Model1 = inla(Frm1,
data = inla.stack.data(Joint.stk),
family = c("gaussian", "poisson"),
verbose = TRUE,
E = inla.stack.data(Joint.stk)$e,
control.fixed = list(prec = 0.001,
prec.intercept = 0.0001),
control.predictor = list(
A = inla.stack.A(Joint.stk),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy="gaussian",
int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE))
summary(Model1)
##
## Call:
## c("inla(formula = Frm1, family = c(\"gaussian\", \"poisson\"), data = inla.stack.data(Joint.stk), ", " E = inla.stack.data(Joint.stk)$e, verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Joint.stk), ", " compute = TRUE, link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.results = list(return.marginals.random = TRUE, ", " return.marginals.predictor = TRUE), control.fixed = list(prec = 0.001, ", " prec.intercept = 1e-04), control.mode = list(restart = TRUE, ", " theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.7480 493.4520 3.6988 501.8989
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept1 0.0619 0.0023 0.0574 0.0619 0.0664 0.0619 0
## intercept2 -3.2863 0.0474 -3.3795 -3.2863 -3.1933 -3.2863 0
## IndF 0.2345 0.0261 0.1832 0.2345 0.2858 0.2345 0
## Pop -0.0540 0.0170 -0.0873 -0.0540 -0.0207 -0.0540 0
##
## Random effects:
## Name Model
## field1 SPDE2 model
## field1c Copy
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for the Gaussian observations 18.964 0.2756 18.429
## Range for field1 363.736 81.2004 226.674
## Stdev for field1 2.732 0.6086 1.705
## Beta for field1c 6.872 0.1308 6.625
## 0.5quant 0.975quant mode
## Precision for the Gaussian observations 18.962 19.513 18.956
## Range for field1 356.049 547.523 341.894
## Stdev for field1 2.675 4.109 2.569
## Beta for field1c 6.868 7.140 6.852
##
## Expected number of effective parameters(std dev): 1710.73(0.00)
## Number of equivalent replicates : 12.81
##
## Deviance Information Criterion (DIC) ...............: 1490.69
## Deviance Information Criterion (DIC, saturated) ....: NaN
## Effective number of parameters .....................: 1727.33
##
## Watanabe-Akaike information criterion (WAIC) ...: 1464.31
## Effective number of parameters .................: 1465.92
##
## Marginal log-Likelihood: -1774.86
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computed
tapply(Model1$dic$local.dic, Model1$dic$family, sum)
## 1 2
## 331.1733 1159.5145
tapply(Model1$waic$local.waic, Model1$dic$family, sum)
## 1 2
## -290.7136 1755.0284
Pred.pnts = Grd.pnts
ModResult = Model1
ModMesh = mesh2D
pLoc = cbind(Pred.pnts@coords[,1], Pred.pnts@coords[,2])
Ap = inla.spde.make.A(ModMesh, loc=pLoc)
Pred.pnts$RF = drop(Ap %*% ModResult$summary.random$field1$mean)
Pred.pnts$RFc = drop(Ap %*% ModResult$summary.random$field1c$mean)
mFixed = ModResult$summary.fix
Pred.pnts$Int = mFixed[1,1]
Pred.pnts$Int2 = mFixed[2,1]
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = RF + RFc +Int + Int2 +
mFixed[3,1]*IndF +
mFixed[4,1]*Pop,
EstF = exp(LP))
Pred.pnts@data = Pred.pnts@data %>%
mutate(LP = #RF + RFc +Int + Int2 +
mFixed[3,1]*IndF +
mFixed[4,1]*Pop,
EstP = exp(LP))
tmp.r = rasterize(Pred.pnts,
SC.r,
"EstF",
background = NA)
Est.r = rasterize(Pred.pnts,
SC.r,
"EstP",
background = NA)
rng = seq(0, 6, 0.01)
Est.r2 = Est.r #Cleaning-up for plot
Est.r2[Est.r < 1] = 0
mCols = brewer.pal(9, "YlOrBr")
cr0 = colorRampPalette((mCols))(n = 256)
cr = colorRampPalette(c(cr0),
bias = 1.75, space = "rgb")
levelplot(Est.r2,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
main = "Relative Abundance (Winter)",
col.regions = cr, at = rng,
colorkey = list(space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(SC.cnt, col = "darkgray", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(SC.st, col = "black", lwd = 0.5))
Wint.bird = subset(MALL.wint, OBS==1)
levelplot(Est.r2,
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
main = "Relative Abundance (Winter)",
col.regions = cr, at = rng,
colorkey = list(space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(SC.cnt, col = "darkgray", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(SC.st, col = "black", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(Wint.bird, col = "black", cex = 1, pch=1))
The code above is repeated for each season below; I’ve hidden the code from view since it is redundant.
Looks like the same locations are used throughout the year, but, the number of birds using them changes.
rng = seq(0, 19, 1)
Results.stk2 = Results.stk
Results.stk2[Results.stk2 < 1] = 0
mCols = brewer.pal(9, "YlOrBr")
cr0 = colorRampPalette((mCols))(n = 256)
cr = colorRampPalette(c(cr0),
bias = 1.75, space = "rgb")
levelplot(Results.stk2,
layout = c(2,2),
margin = FALSE,
xlab = NULL,
ylab = NULL,
maxpixels = 1e5,
main = " ",
col.regions = cr, at = rng,
colorkey = list(space = "bottom"),
par.strip.text = list(fontface='bold', cex=1.5),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(SC.cnt, col = "darkgray", lwd = 0.5)) +
latticeExtra::layer(sp.polygons(SC.st, col = "black", lwd = 0.5))
max(values(Results.stk[[1]]), na.rm=T) #Maximum relative abundance for Winter
## [1] 5.971139
max(values(Results.stk[[2]]), na.rm=T) #Maximum relative abundance for Spring
## [1] 9.382087
max(values(Results.stk[[3]]), na.rm=T) #Maximum relative abundance for Breeding
## [1] 8.129089
max(values(Results.stk[[4]]), na.rm=T) #Maximum relative abundance for Fall
## [1] 18.4683
writeRaster(Results.stk[[1]], "./Winter.tif")
writeRaster(Results.stk[[2]], "./Spring.tif")
writeRaster(Results.stk[[3]], "./Breeding.tif")
writeRaster(Results.stk[[4]], "./Fall.tif")