Fur seal analysis from Devin S. Johnson, Mevin B. Hooten, and Carey E. Kuhn (2013) Estimating animal resource selection from telemetry data using point process models. Journal of Animal Ecology 82:1155--1164.

run_nfs_analysis_pop_INLA.R

library(sp)
library(rgdal)
library(RColorBrewer)
library(spdep)
library(INLA)
library(grid)
library(gridExtra)
Greys <- colorRampPalette(brewer.pal(9,"Greys"))
#source("stpp_rsf_helper.R")
load("nfs_analysis/JAE_analysis_data_spatial.RData")
proj4string(aug.hab) <-  "+proj=aeqd +lat_0=57.125 +lon_0=-170.284167 +ellps=WGS84"
aug.hab$dist.stp <- spDistsN1(aug.hab, c(0,0))/1000
aug.hab.poly <- as(aug.hab, "SpatialPolygons")
tracks.cnt <- SpatialPointsDataFrame(as(tracks.sp, "SpatialPoints"), data=data.frame(counts=rep(1,nrow(tracks.sp))))
proj4string(tracks.cnt) <-  "+proj=aeqd +lat_0=57.125 +lon_0=-170.284167 +ellps=WGS84"
aug.hab$count <- over(aug.hab.poly, tracks.cnt, fn=sum)[,1] 
aug.hab$cell <- 1:nrow(aug.hab)
aug.hab$ln.kern.res <- lm(pmax(-30, log(kern$pop.kern)) ~ ln.npp + sst + trans.pollock + dist.stp, data=aug.hab@data)$residuals
image(aug.hab, "ln.kern.res")
nb2INLA(file="aug.hab.graph", poly2nb(aug.hab.poly))



### ICAR model
formula.car <- count ~ ln.npp + sst + trans.pollock + ln.kern.res + dist.stp +
              f(cell,model="besag",graph="aug.hab.graph") 
model.car <- inla(formula.car, family="poisson", data=aug.hab@data, 
              control.inla=list(h=0.05),
              control.predictor=list(compute=TRUE))
aug.hab$avail <- model.car$summary.random$cell[,4] + 
                 as.matrix(aug.hab@data[,c("ln.kern.res","dist.stp")])%*%model.car$summary.fixed[5:6,4]
aug.hab$select <- model.car$summary.fixed[1,4] + as.matrix(aug.hab@data[,c("ln.npp","sst","trans.pollock")])%*%model.car$summary.fixed[2:4,4]
image(aug.hab[as.vector(aug.hab@data$select>0),], "select", col=Greys(1000))
plot(ak,add=TRUE, col=1)
model.car$summary.fixed

### Nonspatial model
model.ns <- inla(count ~ ln.npp + sst + trans.pollock + dist.stp, family="poisson", data=aug.hab@data)
model.ns$summary.fixed

#custom function to create the clipping region
PolygonFromBbox<-function(b) {
  x <- c(b[1,1],b[1,2],b[1,2],b[1,1],b[1,1])
  y <- c(b[2,2],b[2,2],b[2,1],b[2,1],b[2,2])
  xy<-cbind(x,y)
  Sr1<-Polygon(coords=xy)
  Srs1<-Polygons(list(Sr1),"s1")
  SpP <- SpatialPolygons(list(Srs1),1:1)
  return(SpP)
}

#custom function to simplify a polygon and eliminate very small polys
fortifi <- function(poly, tol, minarea=NA) {
  require(rgeos)
  require(sp)
  require(ggplot2)
  if(!is.na(minarea)) {
    poly<-poly[gArea(poly,byid=TRUE)>minarea]
  }
  poly<-gSimplify(poly,tol=tol,topologyPreserve=TRUE)
  l<-length(poly)
  poly<-SpatialPolygonsDataFrame(poly,data=data.frame(seq(1,l,1)),match.ID=FALSE)
  names(poly)[1]<-'region'
  slot(poly, "polygons") <- lapply(slot(poly, "polygons"), checkPolygonsHoles)
  poly<-fortify(poly)
  return(poly)
}


### SOME PLOTS
saPoly <- gUnaryUnion(as(aug.hab, "SpatialPolygons"))
box <- bbox(aug.hab)
box[,1] <- box[,1] - 90000
box[,2] <- box[,2] + 90000
clip<-PolygonFromBbox(box)
xlim <- box[1,]
ylim <- box[2,]
proj4string(clip)<-CRS(proj4string(aug.hab))
proj4string(ak) <-  "+proj=aeqd +lat_0=57.125 +lon_0=-170.284167 +ellps=WGS84"
ak_clip<-gIntersection(clip,ak)
ak_fort<-fortifi(ak_clip,tol=1000,minarea=1000000)
selectPoly <- gUnaryUnion(as(aug.hab[as.vector(aug.hab@data$select>0),], "SpatialPolygons"))

### Some colors for mapping
Reds <- colorRampPalette(brewer.pal(9,"Reds"))
Greys <- colorRampPalette(brewer.pal(9,"Greys"))
Blues <- colorRampPalette(brewer.pal(9,"Blues"))
Oranges <- colorRampPalette(brewer.pal(9,"Oranges"))
blueGreen <- colorRampPalette(brewer.pal(3,"BuGn"))
heat <- colorRampPalette(heat_hcl(10))
cmcol <- diverge_hcl(9, h = c(180, 330), c = 59, l = c(75, 95))

### Themes
## bw theme
theme_ocean_map_bw <- theme_grey() +
  theme(panel.background = element_rect(fill=NA)) +
  theme(panel.border = element_rect(colour = "black",fill=NA)) +
  theme(axis.text.y = element_blank(), axis.text.x = element_blank()) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x=element_text(size = rel(1)), axis.title.y=element_text(size = rel(1))) +
  theme(axis.ticks = element_blank()) +
  theme(plot.title=element_text(size = rel(1.5))) +
  theme(legend.position="bottom", legend.title=element_text(size=rel(1), face="italic")) +
  theme(legend.text=element_text(size=rel(1)))
library(vcd)

p <- ggplot() + geom_raster(data=as(aug.hab[as.vector(aug.hab@data$select>0),],"data.frame"), aes(s1, s2, fill = select)) +
  geom_path( data=fortifi(saPoly, tol=1000,minarea=1000000), aes(x=long, y=lat, group = group, id=id)) + 
  #geom_point(data=as(tracks.sp, "data.frame"), aes(x=longitude, y=latitude)) +
  scale_x_continuous(name='east') + scale_y_continuous('north') + 
  scale_fill_gradientn(colours=rev(heat_hcl(100)), #brewer.pal(9,"Greys"),
                       guide=guide_colorbar(label=TRUE,barwidth=10, title="log w(s)", title.position="top")) + 
  coord_fixed(xlim=xlim, ylim=ylim) +
  geom_polygon( data=ak_fort, aes(x=long, y=lat, group = group, id=id),fill=grey(0) ) +
  theme_ocean_map_bw + theme(plot.margin = unit(c(0,1,0,0), "lines")) + 
  ggtitle("(a) Posterior median selection surface\n")

pbw <- ggplot() + geom_raster(data=as(aug.hab[as.vector(aug.hab@data$select>0),],"data.frame"), aes(s1, s2, fill = select)) +
  geom_path( data=fortifi(saPoly, tol=1000,minarea=1000000), aes(x=long, y=lat, group = group, id=id)) + 
  #geom_point(data=as(tracks.sp, "data.frame"), aes(x=longitude, y=latitude)) +
  scale_x_continuous(name='east') + scale_y_continuous('north') + 
  scale_fill_gradientn(colours=Greys(100), #brewer.pal(9,"Greys"),
                       guide=guide_colorbar(label=TRUE,barwidth=10, title="log w(s)", title.position="top")) + 
  coord_fixed(xlim=xlim, ylim=ylim) +
  geom_polygon( data=ak_fort, aes(x=long, y=lat, group = group, id=id),fill=grey(0) ) +
  theme_ocean_map_bw + theme(plot.margin = unit(c(0,1,0,0), "lines")) + 
  ggtitle("(a) Posterior median selection surface\n")


p2 <- ggplot() + geom_raster(data=as(aug.hab[as.vector(aug.hab@data$avail>-30),],"data.frame"), aes(s1, s2, fill = avail)) +
  geom_path( data=fortifi(saPoly, tol=1000,minarea=1000000), aes(x=long, y=lat, group = group, id=id)) + 
  #geom_point(data=as(tracks.sp, "data.frame"), aes(x=longitude, y=latitude)) +
  scale_x_continuous(name='east') + scale_y_continuous('north') + 
  scale_fill_gradientn(colours=Blues(100),
                       guide=guide_colorbar(label=TRUE,barwidth=10, title="log g(s)", title.position="top")) + 
  coord_fixed(xlim=xlim, ylim=ylim) +
  geom_polygon( data=ak_fort, aes(x=long, y=lat, group = group, id=id),fill=grey(0) ) +
  theme_ocean_map_bw + theme(plot.margin = unit(c(0,1,0,0), "lines")) +
  ggtitle("(b) Posterior median availability surface\n")

p2bw <- ggplot() + geom_raster(data=as(aug.hab[as.vector(aug.hab@data$avail>-30),],"data.frame"), aes(s1, s2, fill = avail)) +
  geom_path( data=fortifi(saPoly, tol=1000,minarea=1000000), aes(x=long, y=lat, group = group, id=id)) + 
  #geom_point(data=as(tracks.sp, "data.frame"), aes(x=longitude, y=latitude)) +
  scale_x_continuous(name='east') + scale_y_continuous('north') + 
  scale_fill_gradientn(colours=Greys(100),
                       guide=guide_colorbar(label=TRUE,barwidth=10, title="log g(s)", title.position="top")) + 
  coord_fixed(xlim=xlim, ylim=ylim) +
  geom_polygon( data=ak_fort, aes(x=long, y=lat, group = group, id=id),fill=grey(0) ) +
  theme_ocean_map_bw + theme(plot.margin = unit(c(0,1,0,0), "lines")) +
  ggtitle("(b) Posterior median availability surface\n")


#png("spatOutput.png", height=400, width=800)
xxx <- arrangeGrob(p,p2,ncol=2)
xxxbw <- arrangeGrob(pbw,p2bw,ncol=2)
ggsave("spatOutput_20130306_color.png", plot=xxx, scale=2.75, width=6, height=3, dpi=1200)
ggsave("spatOutput_20130306_color_loRes.png", plot=xxx, scale=2.75, width=6, height=3, dpi=144)
ggsave("spatOutput_20130306_bw.png", plot=xxxbw, scale=2.75, width=6, height=3, dpi=1200)
ggsave("spatOutput_20130306_bw_loRes.png", plot=xxxbw, scale=2.75, width=6, height=3, dpi=144)