dat4bugs <- function(data, timestep = 1, tod = F) { # # Code to convert typical Argos data into proper format for WinBUGS # # Created by Ian Jonsen # # created on: 01/04/2022 # last modified on: 02/25/2006 # # # Note data file must include the following named columns in the first 5 positions: date, time, lc, lat, lon # # date must be a numeric object of class "dates times", this can be accomplished as follows: # data$date <- chron(as.character(data$date), format='d/m/y', out.format='m/d/y') # where format (eg, ='d/m/y') must be set to whatever format the date is in, and out.format is set to 'm/d/y' # # time must be converted to "minutes after midnight" format, this can be accomplished as follows: # data$time <- as.numeric(times(as.character(data$time))) * 1440 # # timestep units are in days # # # The tod argument determines whether time of day is relevant, if not, timesteps are converted to minutes from t = # 0. This allows timesteps < 1 day to be used. # # try to convert date and time to correct format if not already so # require(chron) if(is.factor(data$date) || is.character(data$date)) data$date <- chron(as.character(data$date), format = 'd/m/y', out.format='m/d/y') else if(class(data$date)[1] != "dates") {stop('Error - date is not in correct format, see comments in dat4bugs')} if(!is.numeric(data$time)) data$time <- as.numeric(times(as.character(data$time))) * 1440 else if(max(data$time) > 1440 || min(data$time) < 0){stop('Error - time is not in correct format, see comments in dat4bugs')} newdata <- data # set first observation to time = 0, if time of day (tod) is not relevant if(tod == F){ newdata$time <- newdata$time - newdata$time[1] newdata$date[newdata$time < 0] <- newdata$date[newdata$time < 0] - 1 newdata$time[newdata$time < 0] <- newdata$time[newdata$time < 0] + 1440 } time <- newdata$time lon <- newdata$lon lat <- newdata$lat date <- newdata$date lc <- as.character(newdata$lc) lc[lc == 3] <- 99 lc[lc == 2] <- 2 lc[lc == 1] <- 3 lc[lc == 0] <- 4 lc[lc == "A"] <- 5 lc[lc == "B"] <- 6 lc[lc == 99] <- 1 newdata$lc <- as.factor(lc) if(tod == T){ # group observations into 'windows' of duration = timestep and create index (idx) for WinBUGS start.date <- date[1] end.date <- date[nrow(newdata)] numday <- end.date - start.date numsteps <- ceiling(numday / timestep) day.steps <- seq(start.date, dates(as.numeric(start.date) + numsteps * timestep), by = timestep) row.out <- sapply(2:numsteps, function(i){row.names(newdata)[date < day.steps[i] & date >= day.steps[i-1]]}) tmp <- sapply(1:length(row.out), function(i){if(length(row.out[[i]]) == 0) row.out[[i]] <- NA else row.out[[i]]}) idx <- cumsum(c(1, sapply(tmp, length))) newdata <- newdata[unlist(tmp), ] newdata <- data.frame(newdata, row.names = NULL) # determine fraction of timestep at which each observation is made, j j <- c() for(i in 1:(length(idx) - 1)){ t <- newdata$time[idx[i]:(idx[i+1]-1)] d <- as.numeric(newdata$date[idx[i]:(idx[i+1]-1)] - newdata$date[idx[i]]) j1 <- (t + d * 1440) / (timestep * 1440) j <- c(j, j1) } } else if(tod == F){ # group observations into 'windows' of duration = timestep and create index (idx) for WinBUGS start.date <- date[1] # end.date <- date[nrow(newdata)] # numday <- as.numeric(end.date - start.date) dtime <- time + 1440 * (as.numeric(date - start.date)) newdata <- data.frame(newdata, dtime) numsteps <- ceiling(max(dtime) / (1440 * timestep)) steps <- seq(timestep * 1440, timestep * 1440 * ceiling(max(dtime) / (timestep * 1440)), by = (timestep * 1440)) row.out <- sapply(2:(numsteps+1), function(i){if(i==2){row.names(newdata)[dtime 2){row.names(newdata)[dtime < steps[i-1] & dtime >= steps[i-2]]} }) tmp <- sapply(1:length(row.out), function(i){if(length(row.out[[i]]) == 0) row.out[[i]] <- NA else row.out[[i]]}) idx <- cumsum(c(1, sapply(tmp, length))) newdata <- newdata[unlist(tmp), ] newdata <- data.frame(newdata, row.names = NULL) # determine fraction of timestep at which each observation is made, j j <- c() for(i in 1:(length(idx) - 1)){ if(i == 1) j1 <- newdata$dtime[idx[i]:(idx[i+1]-1)] / (timestep * 1400) else if(i > 1){ j1 <- (newdata$dtime[idx[i]:(idx[i+1]-1)] - steps[i-1]) / (timestep * 1440) } j <- c(j, j1) } } # set up t-distribution parameters and create look-up values for the different location classess sigma.lon <- c() sigma.lat <- c() nu.lon <- c() nu.lat <- c() # parameters for t-distributions (obtained by MLE); estimated from Vincent et al. 2002 Argos data # these values are in Appendix A, Table A1 sigma.lon[1] <- 0.2898660 sigma.lon[2] <- 0.3119293 sigma.lon[3] <- 0.9020423 sigma.lon[4] <- 2.1625936 sigma.lon[5] <- 0.5072920 sigma.lon[6] <- 4.2050261 sigma.lat[1] <- 0.1220553 sigma.lat[2] <- 0.2605126 sigma.lat[3] <- 0.4603374 sigma.lat[4] <- 1.607056 sigma.lat[5] <- 0.5105468 sigma.lat[6] <- 3.041276 nu.lon[1] <- 3.070609 nu.lon[2] <- 1.220822 nu.lon[3] <- 2.298819 nu.lon[4] <- 0.9136517 nu.lon[5] <- 0.786954 nu.lon[6] <- 1.079216 nu.lat[1] <- 2.075642 nu.lat[2] <- 6.314726 nu.lat[3] <- 3.896554 nu.lat[4] <- 1.010729 nu.lat[5] <- 1.057779 nu.lat[6] <- 1.331283 # add following code to restrict nu to be >= 2; required for WinBUGS nu.lon[nu.lon < 2] <- 2 nu.lat[nu.lat < 2] <- 2 # convert sigma's from km back to degrees sigma.lon <- (sigma.lon/6366.71 * 180)/pi sigma.lat <- (sigma.lat/6366.71 * 180)/pi # convert standard errors to precisions (se^-2); required for WinBUGS itau2.lon <- sigma.lon[newdata$lc] ^ -2 itau2.lat <- sigma.lat[newdata$lc] ^ -2 # Note following code deals with timesteps that have no observations by interpolating a single observation at mid-day # this could also be addressed by putting priors on the missing data for itau2, nu, and j. # # convert NA's in itau to Argos class B standard errors (= small precisions) itau2.lon[is.na(itau2.lat)] <- min(itau2.lon, na.rm=T) itau2.lat[is.na(itau2.lat)] <- min(itau2.lat, na.rm=T) # converts NA's in nu to smallest nu allowed by Winbugs nu.lon <- nu.lon[newdata$lc] nu.lat <- nu.lat[newdata$lc] nu.lon[is.na(nu.lon)] <- 2 nu.lat[is.na(nu.lat)] <- 2 # set missing observations to be at mid-day j[is.na(j)] <- 0.5 return(list(y = cbind(newdata$lon, newdata$lat), itau2 = cbind(itau2.lon, itau2.lat), nu = cbind(nu.lon, nu.lat), idx = idx, j = j, RegN = length(idx) - 1)) }