Building data for fur seal analysis in: 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_data_assimilation.R

library(sp)
library(rgdal)
library(raster)
library(gstat)
library(rgeos)
library(automap)
library(adehabitatHR)
source("stpp_rsf_helper.R")

# Read in telemtry

tracks <- read.csv("nfs_analysis/2010_gps_filtered.csv")
tracks$gmt <- as.POSIXct(strptime(tracks$gmt, '%m/%d/%Y %H:%M'), tz="GMT")
tracks$id.trip<- paste(tracks$Id, tracks$tripno, sep="-")
coordinates(tracks) <- ~longitude + latitude
proj4string(tracks) <- CRS("+proj=longlat")
azed <- CRS("+proj=aeqd +lat_0=57.125 +lon_0=-170.284167")
tracks <- spTransform(tracks, azed)
tracks <- tracks[!tracks$Id%in%c("SP1008","SP1022"),]
# plot(tracks, pch=20, cex=0.3)


# Read in dive data
# dive <- read.csv("nfs_analysis/2010_dive_all.csv")
# dive$gmt <- as.POSIXct(strptime(dive$date_time, '%m/%d/%Y %H:%M'), tz="GMT")
# dive <- dive[(dive$depth>5 & dive$wigg>=5) | dive$depth>10,]

# Read in habitat boundary data

ak<-readOGR(dsn="/Users/Devin.Johnson/work/base_gis/alaska",layer="alaska_dnr")
ak<-spTransform(ak,CRS(proj4string(tracks)))
plot(ak, add=TRUE, col=gray(0.7))

shelf <- readOGR(dsn="/Users/Devin.Johnson/work/base_gis/ebsshelf7",layer="ebsshelf7")
shelf <-spTransform(shelf,CRS(proj4string(tracks)))
# plot(shelf, add=TRUE)

#clip track data to self animals that started trips in August
tracks.sp <- tracks[tracks$island=="SP",]
ind <- as.vector(is.na(over(tracks.sp, shelf)) & coordinates(tracks.sp)[,2]< -7145.89)
outs <- unique(tracks.sp[ind,]@data$id.trip)
tracks.sp <- tracks.sp[!tracks.sp@data$id.trip%in%outs,]
ins <- unique(tracks.sp@data$id.trip[months(as.POSIXlt(tracks.sp@data$gmt))=="August"])
tracks.sp <- tracks.sp[tracks.sp@data$id.trip%in%ins,]
tracks.sp@data$time.hr <- as.numeric(tracks.sp@data$gmt)/3600
tracks.sp <- tracks.sp[tracks.sp@data$tripno==1,]

plot(shelf)
plot(tracks.sp, pch=20, cex=0.4, add=TRUE, col="red")
plot(ak, add=TRUE, col=gray(0.3))

# Extract "foraging period" locations
id <- unique(tracks.sp@data[,c("Id","id.trip")])
track.i <- tracks.sp[tracks.sp@data$id.trip==id[1,2],]
dive.i <- dive[dive$Id==id[1,1] & dive$gmt>=min(track.i@data$gmt) & dive$gmt<=max(track.i@data$gmt),]
dive.hours <- unique(trunc(dive.i$gmt, "hour"))
idx <- trunc(track.i$gmt,"hour")%in%dive.hours
tracks.sp.dive <- track.i[idx,]
for(i in 2:nrow(id)){
  track.i <- tracks.sp[tracks.sp@data$id.trip==id[i,2],]
  dive.i <- dive[dive$Id==id[i,1] & dive$gmt>=min(track.i@data$gmt) & dive$gmt<=max(track.i@data$gmt),]
  dive.hours <- unique(trunc(dive.i$gmt, "hour"))
  idx <- trunc(track.i$gmt,"hour")%in%dive.hours
  tracks.sp.dive <- rbind(tracks.sp.dive, track.i[idx,])
}

# Read in som environmental data

grdpts <- makegrid(shelf, cellsize=c(9000,9000))
coordinates(grdpts) <- ~x1+x2
proj4string(grdpts) <- CRS(proj4string(tracks))
grdpts <- raster(as(grdpts, "SpatialPixels"))

#Net Primary productivity
aug.npp <- raster("nfs_analysis/aug2010npp.txt", crs="+proj=longlat")
aug.npp.vals <- getValues(aug.npp)
values(aug.npp) <- ifelse(aug.npp.vals< -5000, NA, aug.npp.vals)
aug.npp <- projectRaster(from=aug.npp, to=grdpts)

# SST
aug.sst <- raster("nfs_analysis/aug2010sst.txt", crs="+proj=longlat +lon_wrap=180")
aug.sst.vals <- getValues(aug.sst)
values(aug.sst) <- ifelse(aug.sst.vals< -4e+33, NA, aug.sst.vals)
aug.sst <- projectRaster(from=aug.sst, to=grdpts)

# RACE survey stuff
ebs <- read.csv("nfs_analysis/ebs2009_2011.csv")
ebs <- ebs[ebs$YEAR==2010,]
bsslope <- read.csv("nfs_analysis/bsslope2002_2010.csv")
bsslope <- bsslope[bsslope$YEAR==2010,]
nbs <- read.csv("nfs_analysis/nbs1982_2010.csv")
nbs <- nbs[nbs$YEAR==2010,]
combo <- rbind(ebs, nbs, bsslope)
sites <- combo[!duplicated(combo$STATION),c(1:6, 12:17)]
sites.sp <- SpatialPointsDataFrame(sites[,c("LONGITUDE","LATITUDE")], data=sites[,-c(1:2)], 
                                   proj4string=CRS("+proj=longlat"))
sites.sp <- spTransform(sites.sp, CRS=CRS(proj4string(tracks.sp)))
raceHull <- gSimplify(gBuffer(sites.sp, width=40000), tol=12800)


# Walleye pollock
pollock <- combo[combo$COMMON=="walleye pollock",c(1:6,8)]
pollock <- merge(pollock, sites, all=TRUE)
pollock$NUMCPUE[is.na(pollock$NUMCPUE)] <- 0
coordinates(pollock) <- ~LONGITUDE+LATITUDE
proj4string(pollock) <- CRS("+proj=longlat")
pollock <- spTransform(pollock, CRS(proj4string(tracks.sp)))
pollock@data$trans.pollock <- pollock@data$NUMCPUE^0.25
mod <- vgm(psill=var(pollock@data$trans.pollock),model="Mat",range=200000,nugget=0, kappa=2)
fit.geo <- fit.variogram(variogram(trans.pollock~1,data=pollock), model=mod)
trans.pollock <- autoKrige(trans.pollock~1, input_data=pollock, new_data=as(grdpts, "SpatialGrid"))[[1]] 


# Bottom temperature
sites.sp@data$BOT_TEMP[sites.sp@data$BOT_TEMP==-9999] <- NA
mod <- vgm(psill=var(sites.sp@data$BOT_TEMP, na.rm=TRUE), model="Sph",range=300000,nugget=0)
bot.temp <- autoKrige(BOT_TEMP~1, input_data=sites.sp[!is.na(sites.sp@data$BOT_TEMP),],
                  new_data=as(grdpts, "SpatialGrid"))[[1]]


# Create SpatialPixel objects 
aug.hab <- as(aug.sst, "SpatialGridDataFrame")
aug.hab$npp <-getValues(aug.npp) 
aug.hab$ln.npp <- log(getValues(aug.npp))
aug.hab$trans.pollock <- trans.pollock$var1.pred
aug.hab$bot.temp <- bot.temp$var1.pred
names(aug.hab@data)[1] <- "sst"
aug.hab <- as(aug.hab, "SpatialPointsDataFrame")
aug.hab <- aug.hab[!is.na(over(aug.hab, raceHull)),]
aug.hab <- as(aug.hab, "SpatialPixelsDataFrame")
aug.hab <- aug.hab[!is.na(aug.hab@data$ln.npp),]
aug.hab <- aug.hab[!is.na(aug.hab@data$sst),]



# Obtain space-time quadrature points and tile areas
idx <- (tracks.sp@data$tripno==1) & (!tracks.sp@data$Id%in%c("SP1008","SP1022"))
tracks.sp <- tracks.sp[idx,]
anim.trip <- unique(tracks.sp@data[,c("Id","id.trip","tripno")])
stQuadList <- list()
for(i in 1:nrow(anim.trip)){
cat("\n Animal-trip: ", i, "/", nrow(anim.trip), "\n")
track.i <- tracks.sp[tracks.sp@data$Id==anim.trip[i,1],]
track.i <- track.i[c(TRUE, apply(diff(coordinates(track.i))==0, 1, sum)!=2),]
stQuadList[[i]] <- SpatTempQuadrature( track.data=track.i, 
                                       environ.data=aug.hab, time.col="time.hr", time.int=1, tile.dim=rep(4000,2))
}

save(tracks.sp, aug.hab, stQuadList, ak, shelf, raceHull, file="nfs_analysis/JAE_analysis_data_st_20120720.RData")


# Obtain spatial quadrature points
kern <- makeKernHR(track.data=tracks.sp, environ.data=aug.hab, h=10800/1.96, time.col="Time", id.col="Id")
spatQuadData <- SpatQuadrature(track.data=tracks.sp, environ.data=aug.hab)
save(tracks.sp, aug.hab, kern, spatQuadData, ak, shelf, raceHull, file="nfs_analysis/JAE_analysis_data_spatial.RData")