In galuardi/SatTagSim: Simulation from Electronic Tags
Simulation of Nova Scotia tagged Atlantic bluefin tuna \emph{(Thunnus thynnus)} and get seasonal transition matrices using an 11-box spatial stratification.
knitr::opts_chunk$set(echo = FALSE, warning = FALSE,
message = FALSE, cache = FALSE,
progress = TRUE, verbose = FALSE
, comment = FALSE
, fig.width = 8
, fig.height = 6
, dev = 'pdf')
load background data and set up parallel
library(SatTagSim)
# library(analyzepsat)
library(rworldmap)
library(rworldxtra)
library(knitr)
data(box7)
data(box11)
data(rmask)
data(nsfish)
data(myramps, package = 'analyzepsat')
map = rworldmap::getMap(resolution = 'high')
# sst = rmask2array(rmask)
# mcoptions = setup.parallel()
Tag data:
par(mar = c(3,3,3,6))
plot(nsfish, col = 1, pch = 19, cex = .6, ylim = c(20, 50), xlim = c(-100, 0))
sapply(1:12, function(x) points(nsfish[nsfish$Month==x,1:2], col = month.colors[x==month.colors[,1],2], pch = 19))
plot(map, add=T, col = 'cornsilk4', fill=T)
degAxis(1)
degAxis(2)
box()
# analyzepsat::.add.month.scale()
\newpage
# length of deployments
ddply(as.data.frame(nsfish), c('TagID'), function(x) data.frame(nmonths = length(unique(x$Month))))
\newpage
make.tag.raster <- function(simdat, xmn = -100, xmx = 30, ymn = 0, ymx = 60, boxsize= 60){
# require(raster)
if(class(simdat)=="SpatialPointsDataFrame"){
simdat = as.data.frame(simdat)
}
r = raster(nrow=(ymx-ymn)*60/boxsize, ncol= (xmx-xmn)*60/boxsize, xmn = xmn, xmx = xmx, ymn = ymn, ymx = ymx)
simdat$CID = 1
coordinates(simdat) = ~Longitude+Latitude
sr = rasterize(simdat, r, field = 'CID', fun='count') #[[1]]
return(sr)
}
r = dlply(as.data.frame(nsfish), 'Month', function(x) make.tag.raster(x))
r = make.tag.raster(nsfish)
library(viridis)
plot(r, col = rev(viridis(100)), ylim = c(20, 50), xlim = c(-100, 0), axes =F)
world(add=T, col = 'cornsilk4', fill=T)
degAxis(1)
degAxis(2)
box()
# pdat = ddply(as.data.frame(nsfish), .(TagID, Length), .fun=function(x) data.frame(DAL = nrow(x)))
#
# ggplot(pdat, aes(TagID, DAL))+geom_bar(stat = 'identity')
#
# ggplot(as.data.frame(nsfish), aes(Longitude, Latitude, group = TagID, colour = MaxTemp))+geom_point(stat = 'identity')+coord_map()+scale_color_continuous(low = 'lightblue', high = 'red4', na.value = 'white')
\newpage
Spatial strata in this package
# library(fields)
# library(raster)
# library(SatTagSim)
# data(box7)
plot(box7, xlim = c(-100, 40), ylim = c(0,65))
world(add=T)
plot(box7, add = T, border = 'salmon', lwd=2)
text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data[,3], cex = 2, col = 'darkred')
degAxis(1)
degAxis(2)
box()
title('7-box')
\newpage
plot(box11, xlim = c(-100, 40), ylim = c(-50,80))
world(add=T)
plot(box11, add = T, border = 'salmon', lwd=2)
lpts = coordinates(box11)
lpts[4,2] = 35
lpts[9,] = c(-25, 25)
text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'darkred')
degAxis(1)
degAxis(2)
box()
title('11-box')
\newpage
Environmental data in this package
library(fields)
library(raster)
library(matlab)
data(rmask)
data(woasst)
# par(mfrow=c(1,2))
par(mar = c(2,2,2,4))
# image.plot(woasst$lon, woasst$lat, apply(woasst$sst, 1:2, mean, na.rm = T), axes = F)
# map(add = T, col = 'white', fill = T)
# # world(add=T)
# degAxis(1)
# degAxis(2)
# box()
# title('World Ocean Atlas mean temperature')
# image.plot(sstmat$lon, sstmat$lat, apply(sstmat$data, 1:2, sum, na.rm=T), xlim = c(-100, 40), ylim = c(0,65))
# world(add=T, col = 'white')
# title('Number of months with suitable surface temperature')
# Do it as a raster!!
sstr = flip(raster(xmn = min(woasst$lon), xmx = max(woasst$lon), ymn = min(woasst$lat), ymx = max(woasst$lat), resolution = c(diff(woasst$lon)[1], diff(woasst$lat)[1]), vals = t(apply(woasst$sst, 1:2, mean, na.rm = T))), 'y')
plot(sstr, col = tim.colors(256), axes = F)
map(add = T, col = 'white', fill = T)
# world(add=T)
degAxis(1)
degAxis(2)
box()
title('World Ocean Atlas mean temperature')
\newpage
Set up the temperature constraint field.
This can be done several ways. The preferred method is to use a running average of suitability, by month, for 12 months. In this manner, there is not a hard limit as to whether a fish may be in the area or not, but does provide a selection criteria that is weighted toward the suitability in that month.
The temperature field provided data(rmask) is a surface temperature suitability for bluefin tuna, based on tag measured temperatures.
dat = nsfish
par(mfrow=c(3,4))
par(mar = c(4,4,4,4))
for(i in 1:12){
plot(ecdf(dat$MaxTemp[dat$Month==i]), main = month.name[i], xlab = '' ,ylab = '')
# hist(dat$MaxTemp[seasidx[[i]]], breaks=30, main=seasons[i], xlim = c(0, 40))
mtemp = mean(dat$MaxTemp[dat$Month==i], na.rm=T)
sdtemp = sd(dat$MaxTemp[dat$Month==i], na.rm=T)
abline(v=mtemp, lwd=2, col=2)
abline(v=quantile(dat$MaxTemp[dat$Month==i], c(0.025, .975),na.rm=T), lwd=2, col=4)
# print(quantile(dat$MaxTemp[dat$Month==i], c(0.025, .975),na.rm=T))
}
Based on the plots, a range of $10-28\,^{\circ}\mathrm{C}$ were chosen
\newpage
Temprature preference raster mask
sstmat = list(lon = sort(unique(coordinates(rmask[[1]])[,1]))
,lat = sort(unique(coordinates(rmask[[1]])[,2]))
# ,data = as.array(rmask)
)
sstdata = array(NA, dim = dim(rmask)[c(2,1,3)])
#------------------------------------#
# mask out gomex in summer, GSL in winter
for(i in 7:9){
rmask[[i]] = mask(rmask[[i]], box11[1,], inverse = T)
rmask[[i]] = mask(rmask[[i]], box11[2,], inverse = T)
}
for(i in 1:3){
rmask[[i]] = mask(rmask[[i]], box11[3,], inverse = T)
}
#------------------------------------#
for(i in 1:(dim(rmask)[3])) sstdata[,,i] = rot90(as.array(rmask)[,,i],3)
sstdata[is.na(sstdata)] = 0
# sstmat$data = apply(sstdata, 1:2, sum, na.rm = T)
# sstdata[is.na(sstdata)] = 0
sstmat$data = sstdata
# image.plot(sstmat$lon, sstmat$lat, apply(sstmat$data, 1:2, sum, na.rm=T), xlim = c(-100, 40), ylim = c(0,65), col = viridis(12), axes = F)
sstmatr = t(flip(raster(apply(sstmat$data, 1:2, sum, na.rm = T)), 'x'))
extent(sstmatr) = extent(sstr)
#raster(xmn = min(sstmat$lon), xmx = max(sstmat$lon), ymn = min(sstmat$lat), ymx = max(sstmat$lat), resolution = c(560, 300), vals = apply(sstmat$data, 1:2, mean, na.rm = T))
par(mfrow=c(1,1))
plot(sstmatr, xlim = c(-100, 40), ylim = c(0,65), col = viridis(12), axes = F)
plot(map, add=T, col = 'white')
# world(add=T, col = 'white')
degAxis(1)
degAxis(2)
box()
title('Number of months with suitable surface temperature')
# rm(sstdata, sstdata2)
\newpage
Monthly temperature mask
#
# sstdata2 = sstdata
# for(i in 1:12){
# idx = c(i-1, i, i+1)
# if(i == 1){
# idx = c(12, 1, 2)
# }
# if(i == 12){
# idx = c(11, 12, 1)
# }
# # else idx = c(i-1, i, i+1)
# sstdata2[,,i] = apply(sstdata[,,idx], 1:2, sum, na.rm = T)
# }
sstdata[sstdata==0] = NA
sstdata2 = sstdata
for(i in 1:12){
idx = c(i, i+1)
if(i == 1){
idx = c(12, 1, 2)
}
if(i == 12){
idx = c(11, 12, 1)
}
else idx = c(i-1, i, i+1)
sstdata2[,,i] = apply(sstdata[,,idx], 1:2, sum, na.rm = T)+sstdata[,,i]
}
sstdata2[is.na(sstdata2)] = 0
sstmat$data = sstdata2
par(mfrow=c(3,4))
for(i in 1:12){
par(mar=c(1,1,2,5))
plot(map, add=F, col = 0, border = 0, xlim = c(-95, 30), ylim = c(20,60), axes = F)
image(sstmat$lon, sstmat$lat, sstmat$data[,,i], xlim = c(-100, 40), ylim = c(0,65), col = viridis(max(sstmat$data[,,i])+1), axes = F, add=T)
plot(map, add=T, col = 'white', axes = F)
image.plot(sstmat$data[,,i], legend.only = T, col = viridis(max(sstmat$data[,,i])+1), nlevel = 3)
box()
# world(add=T, col = 'white')
degAxis(1)
degAxis(2)
# degAxis(1, tck = 0.02, las = 2, hadj = -1)
# degAxis(2, tck = 0.02, las = 1, hadj = -1)
title(month.name[i])
}
# title('Number of months with suitable surface temperature')
\newpage
# SET UP THE WOA SST AS A 3D ARRAY
# Function does not work!!!
# sstmat = rmask2array(rmask)
sstmat = list(lon = sort(unique(coordinates(rmask[[1]])[,1]))
,lat = sort(unique(coordinates(rmask[[1]])[,2]))
# ,data = as.array(rmask)
)
sstdata = array(NA, dim = dim(rmask)[c(2,1,3)])
for(i in 1:(dim(rmask)[3])) sstdata[,,i] = rot90(as.array(rmask)[,,i],3)
# sstmat$data = apply(sstdata, 1:2, sum, na.rm = T)
sstmat$data = sstdata
## If using the running average suitability
{
sstm = sstdata
for(i in c(2:6, 10:11)){
sstm[,,i] = apply(sstdata[,,c(i-1,i,i+1)], 1:2, sum, na.rm = T)
}
sstm[,,1] = apply(sstdata[,,c(12, 1, 2)], 1:2, sum, na.rm = T)
sstm[,,12] = apply(sstdata[,,c(11, 12, 1)], 1:2, sum, na.rm = T)
sstm[,,7] = sstm[,,7]*2 # purely for the summer in the gulf of mexico..
sstm[,,8] = sstm[,,8]*2
sstm[,,9] = sstm[,,9]*2
for(i in 7:9){
sstm[,,i][is.na(sstm[,,i])] = 0
}
sstmat$data = sstm
rm(sstdata, sstm)
}
par(mfrow=c(3,4))
for(i in 1:12){
par(mar=c(2,2,2,2))
image(sstmat$lon, sstmat$lat, sstmat$data[,,i], col = tim.colors(4))
world(add=T)
}
Get parameters from tag data
#-----------------------------------------------------------------#
# add spatial overlay to tag data
nsfish@proj4string = box11@proj4string
nsfish$box = over(nsfish, box11)$ID
par2 = daply(as.data.frame(nsfish), c('box', 'TagID', 'Month'), function(x) get.uv(x[,c('Day','Month','Year','Longitude','Latitude')]))
tagwts = daply(as.data.frame(nsfish), c('box', 'Month'), function(x) nrow(x))
ubox = apply(par2[,,,1], 3, rowMeans, na.rm=T)*-1
vbox = apply(par2[,,,2], 3, rowMeans, na.rm=T)
#-----------------------------------------------------------------#
par.ns = get.allpar(as.data.frame(nsfish))
d.ns = get.kfD(as.data.frame(nsfish))
simpar = merge.par(par.ns, d.ns, return.mean = T)
simpar = make.par.array(tracks = nsfish, inbox = box11, rasbox = NULL, rrows = 26*5, rcols = 29*5, use_wts = NULL, missvec = c(3,7,10,11), fillvec = c(4, 8, 9, 9))
# par.ns = get.allpar(as.data.frame(nsfish))
# d.ns = get.kfD(as.data.frame(nsfish))
# simpar = merge.par(par.ns, d.ns, return.mean = T)
knitr::kable(simpar, caption = 'Advection and Diffusion parameters for simulation. For these fish, a fixed diffusion was used during track estimation so there is no variance in D')
\newpage
Examine the tag based Diffusion parameters
pdat = ddply(as.data.frame(nsfish), .(TagID), .fun=function(x) c(x$Length[1], x$D[1], x$Dsd[1]))
names(pdat)[2:4] = c('CFL','D','Dsd')
ggplot(pdat, aes(x=CFL, y=D)) +
geom_errorbar(aes(ymin=D-Dsd, ymax=D+Dsd), colour="black", width=.1) +
geom_abline(aes(intercept=mean(D, na.rm=T), slope=0, width = 2), colour='red')+
# geom_abline(aes(intercept=mean(D+Dsd, na.rm=T), slope=0, colour='blue', width = 2))+
# geom_abline(aes(intercept=mean(D-Dsd, na.rm=T), slope=0, colour='blue', width = 2))+
geom_point( size=3, shape=21, fill="white") + # 21 is filled circle
# theme(axis.text.x = element_blank())+
ggtitle(paste0("D distributions for ",nrow(pdat) ," tagged fish") )
\newpage
Examine the tag based Advection parameters
par.ns = get.allpar(as.data.frame(nsfish))
plot.uv.density(par.ns)
\newpage
Setup simulation parameters
#------------------------------------------------------------------------#
# Set number of sims, simulation steps per month and number of years
msims = 50 # simulations to be started in each month
npmon = 4 # number of simulation steps per month.
nyears = 2 # number of years to simulate for each track
nreps = 167 # number of replicates for variance calculation
sstol = 2 # number of suitable months for each daily simulation step
#------------------------------------------------------------------------#
simvals = data.frame(sim_param = c('msims','npmon','nyears','sstol'), value = c(msims, npmon, nyears, sstol), description = c('simulations to be started in each month','number of simulation steps per month','number of years to simulate for each track','number of suitable months for each daily simulation step'))
kable(simvals, caption = "Control parameters for simulation")
#------------------------------------------------------------------------#
# Set month order of simulaitons to correspond to release month
morder = array(rep(1:12, 12), dim=c(12, 12))
for(i in 1:12){
if(i==1) morder[i,] = 1:12
else morder[i,] = c(i:12,1:((1:12)[i]-1))
}
#------------------------------------------------------------------------#
morder = data.frame(morder, row.names = month.abb)
names(morder) = paste0('m',1:12)
kable(morder, caption = 'Month order for multi-month start simulation', row.names = T)
\newpage
Set up starting points for multi-start month simulation
Set any indices for size or other criteria here. This may also be done prior to the parameter calculation step.
ns = as.data.frame(nsfish)
ns = ns[,c('TagID','Day','Month','Year','Longitude','Latitude')]
spts = get.start.pts(ns, msims, months = 1:12, posnames = c('Longitude','Latitude'))
# str(spts)
# plot(map, xlim = c(-100, 0), ylim = c(20, 50))
# for(i in 1:12) points(spts[[i]], col = month.colors[which(as.numeric(month.colors[,1])==i),2], pch = 19)
# degAxis(1)
# degAxis(2)
# box()
\newpage
rasbox = make.rasbox(box11, raster = T)
# boxmat = make.rasbox(box11, raster = F)
# fix the bay of biscay
rasbox[cellFromXY(rasbox, c(-2.5, 42.5))] = 8
rasbox[cellFromXY(rasbox, c(-2.5, 45.5))] = 8
# make a matrix for speed
boxmat = list()
boxmat$lon = unique(coordinates(rasbox)[,1])
boxmat$lat = sort(unique(coordinates(rasbox)[,2]))
boxmat$box = t(as.matrix(flip(rasbox, 2)))
\newpage
Run a single simulation in parallel
mcoptions = setup.parallel()
#---------------------------------------------------------------#
# SIMULATIONS
print(paste0('simulating ', length(spts)*msims, ' tracks for ', nyears, ' years'))
stime = Sys.time()
simdat = list()
ncores = detectCores()
cl = makeCluster(ncores)
registerDoParallel(cl, cores = ncores)
for(i in 1:12){
subsp = spts[[i]][sample(1:nrow(spts[[i]]), msims, replace = T),]
subsp$row = 1:nrow(subsp)
sp = dlply(subsp, 'row', function(x) x[,1:2])
simorder = as.numeric(morder[i,])
print(paste0('simulating ', length(sp), ' tracks starting in ', month.name[i]))
# test = make.sim.track.par(tpar = simpar, morder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions)
test = make.sim.track.par2(par_array = simpar, boxmat = boxmat, simorder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions)
simdat[[i]] = test
runtime = Sys.time()-stime
print(paste0('elapsed time: ', runtime))
}
stopCluster(cl)
rm(test, sp, subsp)
simdat = unlist(simdat, recursive = F)
stime = Sys.time()
mcoptions = setup.parallel()
for(j in 58:nreps){
simdat = list()
#---------------------------------------------------------------#
# SIMULATIONS
# print(paste0('simulating ', length(spts)*msims, ' tracks for ', nyears, ' years'))
print(paste0('simulating ', msims*12, ' tracks: ',j ,' of ', nreps, ' times'))
# stime = Sys.time()
simdat = list()
ncores = detectCores()
cl = makeCluster(ncores)
registerDoParallel(cl, cores = ncores)
for(i in 1:12){
subsp = spts[[i]][sample(1:nrow(spts[[i]]), msims, replace = T),]
subsp$row = 1:nrow(subsp)
sp = dlply(subsp, 'row', function(x) x[,1:2])
simorder = as.numeric(morder[i,])
print(paste0('simulating ', length(sp), ' tracks starting in ', month.name[i]))
# test = make.sim.track.par(tpar = simpar, morder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions)
test = make.sim.track.par2(par_array = simpar, boxmat = boxmat, simorder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions)
simdat[[i]] = test
runtime = Sys.time()-stime
print(paste0('elapsed time: ', runtime))
simdat[[i]] = test
}
stopCluster(cl)
simdat = unlist(simdat, recursive = F)
# save('simdat' , file = paste0('VAR_SIMS_WABFT_LARGE/WABFT_1008_', j, '.Rdata'))
save('simdat' , file = paste0('SIM_NS_PC', j, '.Rdata'))
runtime = Sys.time()-stime
print(paste0('elapsed time: ', runtime))
}
runtime = Sys.time()-stime
write.table(runtime, file = 'pc_runtime.txt')
\newpage
## Make a dataframe and raster of results
```r
# MAKE A DATA FRAME OF THE RESULTS
simdatdf = ldply(simdat)
seasons = c("Winter", "Spring", "Summer", "Fall")
# seasidx = make.seas.idx(simdatdf)
simdatdf$seas[simdatdf$Month%in%c(1,2,3)] = 1
simdatdf$seas[simdatdf$Month%in%c(4,5,6)] = 2
simdatdf$seas[simdatdf$Month%in%c(7,8,9)] = 3
simdatdf$seas[simdatdf$Month%in%c(10,11,12)] = 4
# rm(seasidx)
names(simdatdf)[which(names(simdatdf)%in%c('lon,','lat'))] = c('lon','lat')
#---------------------------------------------------------------#
# Make Rasters of Results
mycol = colorRampPalette(c("lightcyan", "royalblue", "blue", "lemonchiffon", "orange", "red"), space = "Lab")
sr = dlply(simdatdf, 'seas', function(x) make.sim.raster(x, boxsize = 150))
par(mfrow=c(2,2))
par(mar=c(2,2,2,2))
for(i in 1:4){
plot(sr[[i]]/max(values(sr[[i]]), na.rm=T), zlim = c(0.01,1), interp=F, col=mycol(256), axes=F)
rpts = xyFromCell(sr[[i]], which(!is.na(values(sr[[i]]))))
world(add=T, col='grey80', fill=F)
# plot(box7, add=T, border='salmon', lwd=2)
# text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data$ID, font=2, cex=1.2)
plot(box11, add=T, border='salmon', lwd=2)
# text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2)
text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black')
degAxis(1)
degAxis(2)
title(seasons[i])
box()
}
\newpage
Get seasonal transition matrices using the 11 box model
# try new version of function... it works!!
datbox = get.first.box(simdat, 2000, box11, seas.len = npmon*3) # season length is number per month*3
boxtrans = get.trans.prob(datbox, nyears=100, adims = c(11,11,4), perc = T)
names(boxtrans) = seasons
plot.boxtrans(boxtrans)
# try new version of function... it works!!
datbox = get.first.box(simdat, 2000, box11, seas.len = npmon*3) # season length is number per month*3
boxtrans = get.trans.prob(datbox, nyears=100, adims = c(11,11,4), perc = F)
names(boxtrans) = seasons
plot.boxtrans(boxtrans)
Plot simulation results by starting month
par(mfrow=c(3,4))
k = 1:msims
for(i in 1:12){
par(mar = c(2,2,2,2))
plot(box11, add=F, border='salmon', lwd=2, axes = F, xlim = c(-100, 40), ylim = c(0, 70))
world(add = T, col = 'grey50')
# text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2)
# text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black')
degAxis(1)
degAxis(2)
title(month.abb[i])
# cidx = which(month.colors$Month==i)
# lapply(k, function(j) lines(simdat[[j]][,1:2], col = as.character(month.colors[cidx,2])))
lapply(k, function(j) lines(simdat[[j]][,1:2], col = hsv(.80, .9, .9, alpha = .15)))
lapply(k, function(j) points(simdat[[j]][1,1:2], col = hsv(.5, .5, .9, alpha = .5), pch = 19, cex = 1))
k = k+50
}
plot raster of results by month
sr = dlply(simdatdf, 'Month', function(x) make.sim.raster(x, boxsize = 150))
par(mfrow=c(3, 4))
par(mar=c(2,2,2,2))
for(i in 1:12){
plot(sr[[i]]/max(values(sr[[i]]), na.rm=T), zlim = c(0.01,1), interp=F, col=mycol(256), axes=F)
rpts = xyFromCell(sr[[i]], which(!is.na(values(sr[[i]]))))
world(add=T, col='grey80', fill=F)
# plot(box7, add=T, border='salmon', lwd=2)
# text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data$ID, font=2, cex=1.2)
plot(box11, add=T, border='salmon', lwd=2)
# text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2)
text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black')
degAxis(1)
degAxis(2)
title(month.abb[i])
box()
}
Compare seasonal changes, spatially
sr = dlply(simdatdf, 'seas', function(x) make.sim.raster(x, boxsize = 150))
plot.seas.change(sr, col = (fields::tim.colors(256)), mar = c(2,2,2,4), box = box7)
Compare monthly changes, spatially
sr = dlply(simdatdf, 'Month', function(x) make.sim.raster(x, boxsize = 150))
plot.seas.change(sr, col = (fields::tim.colors(256)), par = c(3,4), mar = c(2,2,2,4), box = box7)
galuardi/SatTagSim documentation built on Nov. 15, 2020, 6:28 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Simulation of Nova Scotia tagged Atlantic bluefin tuna \emph{(Thunnus thynnus)} and get seasonal transition matrices using an 11-box spatial stratification.
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, cache = FALSE, progress = TRUE, verbose = FALSE , comment = FALSE , fig.width = 8 , fig.height = 6 , dev = 'pdf')
load background data and set up parallel
library(SatTagSim) # library(analyzepsat) library(rworldmap) library(rworldxtra) library(knitr) data(box7) data(box11) data(rmask) data(nsfish) data(myramps, package = 'analyzepsat') map = rworldmap::getMap(resolution = 'high') # sst = rmask2array(rmask) # mcoptions = setup.parallel()
Tag data:
par(mar = c(3,3,3,6)) plot(nsfish, col = 1, pch = 19, cex = .6, ylim = c(20, 50), xlim = c(-100, 0)) sapply(1:12, function(x) points(nsfish[nsfish$Month==x,1:2], col = month.colors[x==month.colors[,1],2], pch = 19)) plot(map, add=T, col = 'cornsilk4', fill=T) degAxis(1) degAxis(2) box() # analyzepsat::.add.month.scale()
\newpage
# length of deployments ddply(as.data.frame(nsfish), c('TagID'), function(x) data.frame(nmonths = length(unique(x$Month))))
\newpage
make.tag.raster <- function(simdat, xmn = -100, xmx = 30, ymn = 0, ymx = 60, boxsize= 60){ # require(raster) if(class(simdat)=="SpatialPointsDataFrame"){ simdat = as.data.frame(simdat) } r = raster(nrow=(ymx-ymn)*60/boxsize, ncol= (xmx-xmn)*60/boxsize, xmn = xmn, xmx = xmx, ymn = ymn, ymx = ymx) simdat$CID = 1 coordinates(simdat) = ~Longitude+Latitude sr = rasterize(simdat, r, field = 'CID', fun='count') #[[1]] return(sr) } r = dlply(as.data.frame(nsfish), 'Month', function(x) make.tag.raster(x)) r = make.tag.raster(nsfish) library(viridis) plot(r, col = rev(viridis(100)), ylim = c(20, 50), xlim = c(-100, 0), axes =F) world(add=T, col = 'cornsilk4', fill=T) degAxis(1) degAxis(2) box()
# pdat = ddply(as.data.frame(nsfish), .(TagID, Length), .fun=function(x) data.frame(DAL = nrow(x))) # # ggplot(pdat, aes(TagID, DAL))+geom_bar(stat = 'identity') # # ggplot(as.data.frame(nsfish), aes(Longitude, Latitude, group = TagID, colour = MaxTemp))+geom_point(stat = 'identity')+coord_map()+scale_color_continuous(low = 'lightblue', high = 'red4', na.value = 'white')
\newpage
Spatial strata in this package
# library(fields) # library(raster) # library(SatTagSim) # data(box7) plot(box7, xlim = c(-100, 40), ylim = c(0,65)) world(add=T) plot(box7, add = T, border = 'salmon', lwd=2) text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data[,3], cex = 2, col = 'darkred') degAxis(1) degAxis(2) box() title('7-box')
\newpage
plot(box11, xlim = c(-100, 40), ylim = c(-50,80)) world(add=T) plot(box11, add = T, border = 'salmon', lwd=2) lpts = coordinates(box11) lpts[4,2] = 35 lpts[9,] = c(-25, 25) text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'darkred') degAxis(1) degAxis(2) box() title('11-box')
\newpage
Environmental data in this package
library(fields) library(raster) library(matlab) data(rmask) data(woasst) # par(mfrow=c(1,2)) par(mar = c(2,2,2,4)) # image.plot(woasst$lon, woasst$lat, apply(woasst$sst, 1:2, mean, na.rm = T), axes = F) # map(add = T, col = 'white', fill = T) # # world(add=T) # degAxis(1) # degAxis(2) # box() # title('World Ocean Atlas mean temperature') # image.plot(sstmat$lon, sstmat$lat, apply(sstmat$data, 1:2, sum, na.rm=T), xlim = c(-100, 40), ylim = c(0,65)) # world(add=T, col = 'white') # title('Number of months with suitable surface temperature') # Do it as a raster!! sstr = flip(raster(xmn = min(woasst$lon), xmx = max(woasst$lon), ymn = min(woasst$lat), ymx = max(woasst$lat), resolution = c(diff(woasst$lon)[1], diff(woasst$lat)[1]), vals = t(apply(woasst$sst, 1:2, mean, na.rm = T))), 'y') plot(sstr, col = tim.colors(256), axes = F) map(add = T, col = 'white', fill = T) # world(add=T) degAxis(1) degAxis(2) box() title('World Ocean Atlas mean temperature')
\newpage
Set up the temperature constraint field.
This can be done several ways. The preferred method is to use a running average of suitability, by month, for 12 months. In this manner, there is not a hard limit as to whether a fish may be in the area or not, but does provide a selection criteria that is weighted toward the suitability in that month.
The temperature field provided data(rmask) is a surface temperature suitability for bluefin tuna, based on tag measured temperatures.
dat = nsfish par(mfrow=c(3,4)) par(mar = c(4,4,4,4)) for(i in 1:12){ plot(ecdf(dat$MaxTemp[dat$Month==i]), main = month.name[i], xlab = '' ,ylab = '') # hist(dat$MaxTemp[seasidx[[i]]], breaks=30, main=seasons[i], xlim = c(0, 40)) mtemp = mean(dat$MaxTemp[dat$Month==i], na.rm=T) sdtemp = sd(dat$MaxTemp[dat$Month==i], na.rm=T) abline(v=mtemp, lwd=2, col=2) abline(v=quantile(dat$MaxTemp[dat$Month==i], c(0.025, .975),na.rm=T), lwd=2, col=4) # print(quantile(dat$MaxTemp[dat$Month==i], c(0.025, .975),na.rm=T)) }
Based on the plots, a range of $10-28\,^{\circ}\mathrm{C}$ were chosen
\newpage
Temprature preference raster mask
sstmat = list(lon = sort(unique(coordinates(rmask[[1]])[,1])) ,lat = sort(unique(coordinates(rmask[[1]])[,2])) # ,data = as.array(rmask) ) sstdata = array(NA, dim = dim(rmask)[c(2,1,3)]) #------------------------------------# # mask out gomex in summer, GSL in winter for(i in 7:9){ rmask[[i]] = mask(rmask[[i]], box11[1,], inverse = T) rmask[[i]] = mask(rmask[[i]], box11[2,], inverse = T) } for(i in 1:3){ rmask[[i]] = mask(rmask[[i]], box11[3,], inverse = T) } #------------------------------------# for(i in 1:(dim(rmask)[3])) sstdata[,,i] = rot90(as.array(rmask)[,,i],3) sstdata[is.na(sstdata)] = 0 # sstmat$data = apply(sstdata, 1:2, sum, na.rm = T) # sstdata[is.na(sstdata)] = 0 sstmat$data = sstdata # image.plot(sstmat$lon, sstmat$lat, apply(sstmat$data, 1:2, sum, na.rm=T), xlim = c(-100, 40), ylim = c(0,65), col = viridis(12), axes = F) sstmatr = t(flip(raster(apply(sstmat$data, 1:2, sum, na.rm = T)), 'x')) extent(sstmatr) = extent(sstr) #raster(xmn = min(sstmat$lon), xmx = max(sstmat$lon), ymn = min(sstmat$lat), ymx = max(sstmat$lat), resolution = c(560, 300), vals = apply(sstmat$data, 1:2, mean, na.rm = T)) par(mfrow=c(1,1)) plot(sstmatr, xlim = c(-100, 40), ylim = c(0,65), col = viridis(12), axes = F) plot(map, add=T, col = 'white') # world(add=T, col = 'white') degAxis(1) degAxis(2) box() title('Number of months with suitable surface temperature') # rm(sstdata, sstdata2)
\newpage
Monthly temperature mask
# # sstdata2 = sstdata # for(i in 1:12){ # idx = c(i-1, i, i+1) # if(i == 1){ # idx = c(12, 1, 2) # } # if(i == 12){ # idx = c(11, 12, 1) # } # # else idx = c(i-1, i, i+1) # sstdata2[,,i] = apply(sstdata[,,idx], 1:2, sum, na.rm = T) # } sstdata[sstdata==0] = NA sstdata2 = sstdata for(i in 1:12){ idx = c(i, i+1) if(i == 1){ idx = c(12, 1, 2) } if(i == 12){ idx = c(11, 12, 1) } else idx = c(i-1, i, i+1) sstdata2[,,i] = apply(sstdata[,,idx], 1:2, sum, na.rm = T)+sstdata[,,i] } sstdata2[is.na(sstdata2)] = 0 sstmat$data = sstdata2 par(mfrow=c(3,4)) for(i in 1:12){ par(mar=c(1,1,2,5)) plot(map, add=F, col = 0, border = 0, xlim = c(-95, 30), ylim = c(20,60), axes = F) image(sstmat$lon, sstmat$lat, sstmat$data[,,i], xlim = c(-100, 40), ylim = c(0,65), col = viridis(max(sstmat$data[,,i])+1), axes = F, add=T) plot(map, add=T, col = 'white', axes = F) image.plot(sstmat$data[,,i], legend.only = T, col = viridis(max(sstmat$data[,,i])+1), nlevel = 3) box() # world(add=T, col = 'white') degAxis(1) degAxis(2) # degAxis(1, tck = 0.02, las = 2, hadj = -1) # degAxis(2, tck = 0.02, las = 1, hadj = -1) title(month.name[i]) } # title('Number of months with suitable surface temperature')
\newpage
# SET UP THE WOA SST AS A 3D ARRAY # Function does not work!!! # sstmat = rmask2array(rmask) sstmat = list(lon = sort(unique(coordinates(rmask[[1]])[,1])) ,lat = sort(unique(coordinates(rmask[[1]])[,2])) # ,data = as.array(rmask) ) sstdata = array(NA, dim = dim(rmask)[c(2,1,3)]) for(i in 1:(dim(rmask)[3])) sstdata[,,i] = rot90(as.array(rmask)[,,i],3) # sstmat$data = apply(sstdata, 1:2, sum, na.rm = T) sstmat$data = sstdata ## If using the running average suitability { sstm = sstdata for(i in c(2:6, 10:11)){ sstm[,,i] = apply(sstdata[,,c(i-1,i,i+1)], 1:2, sum, na.rm = T) } sstm[,,1] = apply(sstdata[,,c(12, 1, 2)], 1:2, sum, na.rm = T) sstm[,,12] = apply(sstdata[,,c(11, 12, 1)], 1:2, sum, na.rm = T) sstm[,,7] = sstm[,,7]*2 # purely for the summer in the gulf of mexico.. sstm[,,8] = sstm[,,8]*2 sstm[,,9] = sstm[,,9]*2 for(i in 7:9){ sstm[,,i][is.na(sstm[,,i])] = 0 } sstmat$data = sstm rm(sstdata, sstm) } par(mfrow=c(3,4)) for(i in 1:12){ par(mar=c(2,2,2,2)) image(sstmat$lon, sstmat$lat, sstmat$data[,,i], col = tim.colors(4)) world(add=T) }
Get parameters from tag data
#-----------------------------------------------------------------# # add spatial overlay to tag data nsfish@proj4string = box11@proj4string nsfish$box = over(nsfish, box11)$ID par2 = daply(as.data.frame(nsfish), c('box', 'TagID', 'Month'), function(x) get.uv(x[,c('Day','Month','Year','Longitude','Latitude')])) tagwts = daply(as.data.frame(nsfish), c('box', 'Month'), function(x) nrow(x)) ubox = apply(par2[,,,1], 3, rowMeans, na.rm=T)*-1 vbox = apply(par2[,,,2], 3, rowMeans, na.rm=T) #-----------------------------------------------------------------# par.ns = get.allpar(as.data.frame(nsfish)) d.ns = get.kfD(as.data.frame(nsfish)) simpar = merge.par(par.ns, d.ns, return.mean = T)
simpar = make.par.array(tracks = nsfish, inbox = box11, rasbox = NULL, rrows = 26*5, rcols = 29*5, use_wts = NULL, missvec = c(3,7,10,11), fillvec = c(4, 8, 9, 9)) # par.ns = get.allpar(as.data.frame(nsfish)) # d.ns = get.kfD(as.data.frame(nsfish)) # simpar = merge.par(par.ns, d.ns, return.mean = T)
knitr::kable(simpar, caption = 'Advection and Diffusion parameters for simulation. For these fish, a fixed diffusion was used during track estimation so there is no variance in D')
\newpage
Examine the tag based Diffusion parameters
pdat = ddply(as.data.frame(nsfish), .(TagID), .fun=function(x) c(x$Length[1], x$D[1], x$Dsd[1])) names(pdat)[2:4] = c('CFL','D','Dsd') ggplot(pdat, aes(x=CFL, y=D)) + geom_errorbar(aes(ymin=D-Dsd, ymax=D+Dsd), colour="black", width=.1) + geom_abline(aes(intercept=mean(D, na.rm=T), slope=0, width = 2), colour='red')+ # geom_abline(aes(intercept=mean(D+Dsd, na.rm=T), slope=0, colour='blue', width = 2))+ # geom_abline(aes(intercept=mean(D-Dsd, na.rm=T), slope=0, colour='blue', width = 2))+ geom_point( size=3, shape=21, fill="white") + # 21 is filled circle # theme(axis.text.x = element_blank())+ ggtitle(paste0("D distributions for ",nrow(pdat) ," tagged fish") )
\newpage
Examine the tag based Advection parameters
par.ns = get.allpar(as.data.frame(nsfish)) plot.uv.density(par.ns)
\newpage
Setup simulation parameters
#------------------------------------------------------------------------# # Set number of sims, simulation steps per month and number of years msims = 50 # simulations to be started in each month npmon = 4 # number of simulation steps per month. nyears = 2 # number of years to simulate for each track nreps = 167 # number of replicates for variance calculation sstol = 2 # number of suitable months for each daily simulation step #------------------------------------------------------------------------# simvals = data.frame(sim_param = c('msims','npmon','nyears','sstol'), value = c(msims, npmon, nyears, sstol), description = c('simulations to be started in each month','number of simulation steps per month','number of years to simulate for each track','number of suitable months for each daily simulation step')) kable(simvals, caption = "Control parameters for simulation") #------------------------------------------------------------------------# # Set month order of simulaitons to correspond to release month morder = array(rep(1:12, 12), dim=c(12, 12)) for(i in 1:12){ if(i==1) morder[i,] = 1:12 else morder[i,] = c(i:12,1:((1:12)[i]-1)) } #------------------------------------------------------------------------# morder = data.frame(morder, row.names = month.abb) names(morder) = paste0('m',1:12) kable(morder, caption = 'Month order for multi-month start simulation', row.names = T)
\newpage
Set up starting points for multi-start month simulation
Set any indices for size or other criteria here. This may also be done prior to the parameter calculation step.
ns = as.data.frame(nsfish) ns = ns[,c('TagID','Day','Month','Year','Longitude','Latitude')] spts = get.start.pts(ns, msims, months = 1:12, posnames = c('Longitude','Latitude')) # str(spts) # plot(map, xlim = c(-100, 0), ylim = c(20, 50)) # for(i in 1:12) points(spts[[i]], col = month.colors[which(as.numeric(month.colors[,1])==i),2], pch = 19) # degAxis(1) # degAxis(2) # box()
\newpage
rasbox = make.rasbox(box11, raster = T) # boxmat = make.rasbox(box11, raster = F) # fix the bay of biscay rasbox[cellFromXY(rasbox, c(-2.5, 42.5))] = 8 rasbox[cellFromXY(rasbox, c(-2.5, 45.5))] = 8 # make a matrix for speed boxmat = list() boxmat$lon = unique(coordinates(rasbox)[,1]) boxmat$lat = sort(unique(coordinates(rasbox)[,2])) boxmat$box = t(as.matrix(flip(rasbox, 2)))
\newpage
Run a single simulation in parallel
mcoptions = setup.parallel() #---------------------------------------------------------------# # SIMULATIONS print(paste0('simulating ', length(spts)*msims, ' tracks for ', nyears, ' years')) stime = Sys.time() simdat = list() ncores = detectCores() cl = makeCluster(ncores) registerDoParallel(cl, cores = ncores) for(i in 1:12){ subsp = spts[[i]][sample(1:nrow(spts[[i]]), msims, replace = T),] subsp$row = 1:nrow(subsp) sp = dlply(subsp, 'row', function(x) x[,1:2]) simorder = as.numeric(morder[i,]) print(paste0('simulating ', length(sp), ' tracks starting in ', month.name[i])) # test = make.sim.track.par(tpar = simpar, morder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions) test = make.sim.track.par2(par_array = simpar, boxmat = boxmat, simorder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions) simdat[[i]] = test runtime = Sys.time()-stime print(paste0('elapsed time: ', runtime)) } stopCluster(cl) rm(test, sp, subsp) simdat = unlist(simdat, recursive = F)
stime = Sys.time() mcoptions = setup.parallel() for(j in 58:nreps){ simdat = list() #---------------------------------------------------------------# # SIMULATIONS # print(paste0('simulating ', length(spts)*msims, ' tracks for ', nyears, ' years')) print(paste0('simulating ', msims*12, ' tracks: ',j ,' of ', nreps, ' times')) # stime = Sys.time() simdat = list() ncores = detectCores() cl = makeCluster(ncores) registerDoParallel(cl, cores = ncores) for(i in 1:12){ subsp = spts[[i]][sample(1:nrow(spts[[i]]), msims, replace = T),] subsp$row = 1:nrow(subsp) sp = dlply(subsp, 'row', function(x) x[,1:2]) simorder = as.numeric(morder[i,]) print(paste0('simulating ', length(sp), ' tracks starting in ', month.name[i])) # test = make.sim.track.par(tpar = simpar, morder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions) test = make.sim.track.par2(par_array = simpar, boxmat = boxmat, simorder = rep(simorder, nyears), sp = sp, bath = bath, sstmat = sstmat, seaslen = npmon, sstol = sstol, mcoptions = mcoptions) simdat[[i]] = test runtime = Sys.time()-stime print(paste0('elapsed time: ', runtime)) simdat[[i]] = test } stopCluster(cl) simdat = unlist(simdat, recursive = F) # save('simdat' , file = paste0('VAR_SIMS_WABFT_LARGE/WABFT_1008_', j, '.Rdata')) save('simdat' , file = paste0('SIM_NS_PC', j, '.Rdata')) runtime = Sys.time()-stime print(paste0('elapsed time: ', runtime)) } runtime = Sys.time()-stime write.table(runtime, file = 'pc_runtime.txt')
\newpage ## Make a dataframe and raster of results ```r # MAKE A DATA FRAME OF THE RESULTS simdatdf = ldply(simdat) seasons = c("Winter", "Spring", "Summer", "Fall") # seasidx = make.seas.idx(simdatdf) simdatdf$seas[simdatdf$Month%in%c(1,2,3)] = 1 simdatdf$seas[simdatdf$Month%in%c(4,5,6)] = 2 simdatdf$seas[simdatdf$Month%in%c(7,8,9)] = 3 simdatdf$seas[simdatdf$Month%in%c(10,11,12)] = 4 # rm(seasidx) names(simdatdf)[which(names(simdatdf)%in%c('lon,','lat'))] = c('lon','lat') #---------------------------------------------------------------# # Make Rasters of Results mycol = colorRampPalette(c("lightcyan", "royalblue", "blue", "lemonchiffon", "orange", "red"), space = "Lab") sr = dlply(simdatdf, 'seas', function(x) make.sim.raster(x, boxsize = 150)) par(mfrow=c(2,2)) par(mar=c(2,2,2,2)) for(i in 1:4){ plot(sr[[i]]/max(values(sr[[i]]), na.rm=T), zlim = c(0.01,1), interp=F, col=mycol(256), axes=F) rpts = xyFromCell(sr[[i]], which(!is.na(values(sr[[i]])))) world(add=T, col='grey80', fill=F) # plot(box7, add=T, border='salmon', lwd=2) # text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data$ID, font=2, cex=1.2) plot(box11, add=T, border='salmon', lwd=2) # text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2) text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black') degAxis(1) degAxis(2) title(seasons[i]) box() }
\newpage
Get seasonal transition matrices using the 11 box model
# try new version of function... it works!! datbox = get.first.box(simdat, 2000, box11, seas.len = npmon*3) # season length is number per month*3 boxtrans = get.trans.prob(datbox, nyears=100, adims = c(11,11,4), perc = T) names(boxtrans) = seasons plot.boxtrans(boxtrans)
# try new version of function... it works!! datbox = get.first.box(simdat, 2000, box11, seas.len = npmon*3) # season length is number per month*3 boxtrans = get.trans.prob(datbox, nyears=100, adims = c(11,11,4), perc = F) names(boxtrans) = seasons plot.boxtrans(boxtrans)
Plot simulation results by starting month
par(mfrow=c(3,4)) k = 1:msims for(i in 1:12){ par(mar = c(2,2,2,2)) plot(box11, add=F, border='salmon', lwd=2, axes = F, xlim = c(-100, 40), ylim = c(0, 70)) world(add = T, col = 'grey50') # text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2) # text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black') degAxis(1) degAxis(2) title(month.abb[i]) # cidx = which(month.colors$Month==i) # lapply(k, function(j) lines(simdat[[j]][,1:2], col = as.character(month.colors[cidx,2]))) lapply(k, function(j) lines(simdat[[j]][,1:2], col = hsv(.80, .9, .9, alpha = .15))) lapply(k, function(j) points(simdat[[j]][1,1:2], col = hsv(.5, .5, .9, alpha = .5), pch = 19, cex = 1)) k = k+50 }
plot raster of results by month
sr = dlply(simdatdf, 'Month', function(x) make.sim.raster(x, boxsize = 150)) par(mfrow=c(3, 4)) par(mar=c(2,2,2,2)) for(i in 1:12){ plot(sr[[i]]/max(values(sr[[i]]), na.rm=T), zlim = c(0.01,1), interp=F, col=mycol(256), axes=F) rpts = xyFromCell(sr[[i]], which(!is.na(values(sr[[i]])))) world(add=T, col='grey80', fill=F) # plot(box7, add=T, border='salmon', lwd=2) # text(coordinates(box7)[,1], coordinates(box7)[,2], box7@data$ID, font=2, cex=1.2) plot(box11, add=T, border='salmon', lwd=2) # text(coordinates(box11)[,1], coordinates(box11)[,2], box11@data$ID, font=2, cex=1.2) text(lpts[,1], lpts[,2], box11@data[,3], cex = 2, col = 'black') degAxis(1) degAxis(2) title(month.abb[i]) box() }
Compare seasonal changes, spatially
sr = dlply(simdatdf, 'seas', function(x) make.sim.raster(x, boxsize = 150)) plot.seas.change(sr, col = (fields::tim.colors(256)), mar = c(2,2,2,4), box = box7)
Compare monthly changes, spatially
sr = dlply(simdatdf, 'Month', function(x) make.sim.raster(x, boxsize = 150)) plot.seas.change(sr, col = (fields::tim.colors(256)), par = c(3,4), mar = c(2,2,2,4), box = box7)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.