Mallard Relative Abundance (South Carolina)

Log-Gaussian Cox Process under Preferential Sampling

date() 
## [1] "Mon Oct 08 15:12:50 2018"

Set Directory

setwd("~/LabWork/Patuxent/USDA")

Clean eBird Data

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)

Get Country Boundaries

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)

Convert Observations to Points and Re-project

MALL.pnt = SpatialPointsDataFrame(MALL[,2:3], MALL)
proj4string(MALL.pnt) = LL84 

MALL.pnt = spTransform(MALL.pnt, nProj)

Select Locations in South Carolina

MALL.pnt$SC = is.na(over(MALL.pnt, SC.st))[,1]
MALL.pnt = subset(MALL.pnt, SC == FALSE)

View Data

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

Spatial Model

Construct Mesh

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) 

Inspect Mesh

mesh2D$n #Count of verticies
## [1] 9502
plot(mesh2D, rgl = F, main = " ")

Get Mesh Nodes

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"

Combine Observations, Nodes, and Dense Grid

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)

Model Covariates

LandFire

https://www.landfire.gov/

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

BioClimate

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)

ENVIREM

http://envirem.github.io/

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)

Scale and Center

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)

Decomposing Variables using Discriminant Analysis

#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)

Human Population Density

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)

Seperate Point Sets by Season

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")

Set-up Winter Model

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)  

Construct Data Stacks for LGCP

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)

LGCP Under Preferential Sampling

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

Relative Abundance Estimates

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)

View Results

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)) 

Same Plot as above with mallard locations layered on top.

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))

Run Model for Spring, Breeding, and Fall Seasons

The code above is repeated for each season below; I’ve hidden the code from view since it is redundant.

View Results

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))  

Seasonal Relative Abundance

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

Save Copies of raster results

writeRaster(Results.stk[[1]], "./Winter.tif")
writeRaster(Results.stk[[2]], "./Spring.tif")
writeRaster(Results.stk[[3]], "./Breeding.tif")
writeRaster(Results.stk[[4]], "./Fall.tif")