SatTagSim: Simulation from Electronic Tags

Documented in make.sim.track.par2

#' Simulate tracks in parallel
#'
#' Function to simulate a single track using monthly parameters, bathymetry and World Ocean Atlas SST in parallel
#' @param par_array array of temporal and spatially explcit simulation parameters @seealso xxxx
#' @param morder month order for simulation (e.g. start in June vs start in January)
#' @param sp list of starting points
#' @param sstmat list of lon, lat and 3D matrix of monthly sea surface temperature. see details
#' @param boxmat raster of spatial strata used when generating movement parameters @seealso xxxx
#' @param seaslen length of each month to simulate. Defaults to 30 (days per month)
#' @param sstol tolerance for sst matching. see details
#' @param mcoptions This is a hidden variable that must be present to run in parallel
#' @details sstmat in this package is a running average of index suitable temperature for each month. Each month has values 0-3. A value of 3, for example, represents areas suitable in the current month, previous month and next month. The sstol flag can make matching criteria more or less restrictive. The default value means that at least two months of suitability, centered on the current month, must be valid for a simulation point.
#' Tracks may be generated by increasing morder (e.g. rep(morder, 2) for a 24 month track)
#' Function uses mcoptions from \code{\link{setup.parallel}} to run in parallel.
#' @return List of simulated tracks. Tracks are data frames with columns: lon lat month
#' @export
#' @author Benjamin Galuardi
#' @examples
#' see vignette
make.sim.track.par2 <- function(par_array = par_array, simorder = simorder, sp = spts, bath = bath, sstmat = sstmat, boxmat = boxmat, seaslen = 30, sstol = 2, mcoptions = setup.parallel(), ...)

  foreach(i = sp, .options.multicore = mcoptions) %dopar% #

  {
    # require(analyzepsat)
    require(SatTagSim)
    ## Subfunctions
    .get.samp <- function (vec, npoints, ci = 0.95)
    {
      vec = as.numeric(vec)
      Sigma <- matrix(vec[1:4], 2, 2) * ci
      mu <- c(vec[5:6])
      if (sum(Sigma) > 0) {
        ndata <- mvrnorm(npoints, mu, Sigma)
        return(ndata)
      }
    }
    .get.bath <- function (lon, lat, BATH)
    {
      X = as.vector(BATH$lon)
      Y = as.vector(BATH$lat)
      xidx = which.min((lon - X)^2)
      yidx = which.min((lat - Y)^2)
      BATH$data[yidx, xidx]
    }
    getsp <- function(sp, var=c(1,0,0,1), bath=bath){
      vec=c(var, sp)
      x = .get.samp(vec, 1)
      while(.get.bath(x[1], x[2], BATH=bath)>0){
        x = .get.samp(vec, 1)
        vec=c(var+c(.5, 0, 0, .5), sp)
      }
      x
    }
    # FUnCTION TO GET A SST MASK VALUE FROM A 3D ARRAY OF WORLD OCEAN ATLAS SST
    get.sst.mask.val <- function (lon, lat, mask, month)
    {
      X = as.vector(mask$lon)
      Y = as.vector(mask$lat)
      xidx = which.min((lon - X)^2)
      yidx = which.min((lat - Y)^2)
      mask$data[xidx, yidx, month] # flagged off for additive months
      # mask$data[xidx, yidx]
    }
    # Convert degrees to radians
    deg2rad <- function(deg) return(deg*pi/180)

    # Calculates the geodesic distance between two points specified by radian latitude/longitude using the
    # Haversine formula (hf)
    gcd.hf <- function(long1, lat1, long2, lat2) {
      R <- 6371 # Earth mean radius [km]
      delta.long <- (long2 - long1)
      delta.lat <- (lat2 - lat1)
      a <- sin(delta.lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta.long/2)^2
      c <- 2 * asin(min(1,sqrt(a)))
      d = R * c
      return(d) # Distance in km
    }

    find.next.sst <- function(lon, lat, sstdf, sstol, expand = 10) {
      # sstdf = data.frame(expand.grid(sstmat$lon, sstmat$lat), sst = as.vector(sstmat$data[,,8]))
      # which(sstdf[,3]>= sstol)
      names(sstdf) = c('lon','lat','sst')

      xlim = c(lon-expand, lon+expand)
      ylim = c(lat-expand, lat+expand)
      xidx  = sstdf$lon>=xlim[1]&sstdf$lon<=xlim[2]
      yidx  = sstdf$lat>=ylim[1]&sstdf$lat<=ylim[2]
      xyidx = which((xidx+yidx)==2)
      pt = c(lon, lat)
      sstdf_sub = sstdf[xyidx, ]
      gidx = which(sstdf_sub[, 3]>= sstol)
      if(length(gidx) > 0){
        # geosphere::distGeo(pt, sstdf_sub[gidx,1:2])
        # dists = geosphere::distGeo(pt, sstdf_sub[gidx,1:2])
        dists = apply(sstdf_sub[gidx,1:2], 1, function(x) gcd.hf(pt[1], pt[1], x[1], x[2]))
        # idxmin = which.min(geosphere::distGeo(pt, sstdf_sub[gidx,1:2]))
        idxmin = sample(which(dists<=quantile(dists, .1)), 1)
        as.numeric(sstdf_sub[gidx[idxmin],1:2])
      }else{
        c(lon, lat)
        }
    }


    # names(temp) = c('lon,','lat','Month')
    if(!is.null(bath)){
      msp = getsp(i, bath = bath)
    }else{
      msp = as.numeric(i)
    }
    uvmult = 30/seaslen

    # object will be the simulated track
    temp = NULL

    # for(n in 1:nrow(sp)){

    # msp = sp

  for(seas in simorder){   # winter, spring, summer, fall

      # browser()
    if(!is.null(sstmat)){
      sstdf = data.frame(expand.grid(sstmat$lon, sstmat$lat), sst = as.vector(sstmat$data[,,seas]))
    }
      # midx = tpar$Month==seas
      # u = c(tpar[midx, 2], tpar[1, 5])*uvmult
      # v = c(tpar[midx, 3], tpar[1, 6])*uvmult
      # D = c(tpar[midx, 4], tpar[1, 7])*uvmult
      # Dorig = D
      # uorig = u
      # vorig = v


    for(j in 1:seaslen) {

      # msps = SpatialPoints(t(as.matrix(msp)))
      # tbox = raster::extract(rasbox, msps)


      xidx = which.min((msp[1] - boxmat$lon)^2)
      yidx = which.min((msp[2] - boxmat$lat)^2)
      tbox = boxmat$box[xidx, yidx]

      if(is.na(tbox) & j == 1){
        if(!is.null(bath)){
          msp = getsp(i, bath = bath)
        }else{
          msp = as.numeric(i)
        }
        xidx = which.min((msp[1] - boxmat$lon)^2)
        yidx = which.min((msp[2] - boxmat$lat)^2)
        tbox = boxmat$box[xidx, yidx]
      }

      if(is.na(tbox) & j > 1){
        msp = temp[j-1, ]
        xidx = which.min((msp[1] - boxmat$lon)^2)
        yidx = which.min((msp[2] - boxmat$lat)^2)
        tbox = boxmat$box[xidx, yidx]
      }

      parbox = as.numeric(attributes(par_array)$dimnames[[1]])
      pbidx = which(parbox==tbox)

      u = c(par_array[pbidx, seas, 1])*uvmult
      v = c(par_array[pbidx, seas, 2])*uvmult
      D = c(par_array[pbidx, seas, 3])#*uvmult
      usd = c(par_array[pbidx, seas, 4])
      vsd = c(par_array[pbidx, seas, 5])
      Dsd = c(par_array[pbidx, seas, 6])

      u = c(u, usd)
      v = c(v, vsd)
      D = c(D, Dsd)

      Dorig = D
      uorig = u
      vorig = v

      ulim = c(-50, 50)*uvmult
      vlim = c(-50, 50)*uvmult
      Dlim = c(0, 5000)#*uvmult

      t1 = SatTagSim::simm.kf(2, u, v, D, msp, ulim, vlim, Dlim)[2,]
#
#         t1s = SpatialPoints(t(as.matrix(t1)))
#         tbox = raster::extract(rasbox, t1s)
#         parbox = as.numeric(attributes(par_array)$dimnames[[1]])
#
#         pbidx = parbox[parbox==tbox]
#
#         u = c(par_array[pbidx, seas, 1])*uvmult
#         v = c(par_array[pbidx, seas, 2])*uvmult
#         D = c(par_array[pbidx, seas, 3])*uvmult
#         usd = c(par_array[pbidx, seas, 4])
#         vsd = c(par_array[pbidx, seas, 5])
#         Dsd = c(par_array[pbidx, seas, 6])
#
#         u = c(u, usd)
#         v = c(v, vsd)
#         D = c(D, Dsd)
#
        # else{
        # ii = 1
        # while(is.na(get.sst.mask.val(t1[1], t1[2], sstmat, seas))){
        # 					cat(seas)
        # t1 = simm.kf(2, u, v, D, msp)[2,]
        # ii = ii+1
        # if(ii>5) {
        # temp = temp[1:(nrow(temp)-j-1),]
        # j = j-1
        # t1 = NA
        # if(is.na(get.sst.mask.val(t1[1], t1[2], sstmat, seas))){
        # pd = pointDistance(t1, xyz[gidx, 1:2], lonlat = T)
        # t1 = as.numeric(xyz[gidx[which.min(pd)],1:2])
        # 		         while(is.na(get.sst.mask.val(t1[1], t1[2], sstmat, seas))){
        #   							ttmp = cbind(rep(t1[1], 100), rep(t1[2], 100))
        # 						    tsamp = t(apply(ttmp, 1, function(x)
        # 						        simm.kf(n = 2, sp = as.numeric(x), u = u, v = v, D = D)[2,]))
        #   							tidx = apply(tsamp, 1, function(y) get.sst.mask.val(y[1], y[2], sstmat, seas))
        # 						    t1 = as.numeric(tsamp[sample(which(!is.na(tidx)), 1),])
        if(!is.null(sstmat)){

          ii=1

        # First, try reversing the advective space.
          if(get.sst.mask.val(t1[1], t1[2], sstmat, seas) < sstol){
            t1 = SatTagSim::simm.kf(2, u = c(-1*u[1], u[2]), v = c(-1*v[1], v[2]), D = c(D[1], 1000), msp, ulim, vlim, Dlim)[2,]
            t1 = as.numeric(t1)
          }
        # Secondly, try using a minimum distance algorithm. This is slow..
          while((get.sst.mask.val(t1[1], t1[2], sstmat, seas)) < sstol){
            tsamp = find.next.sst(t1[1], t1[2], sstdf, sstol, expand = 5*ii)
            t1 = as.numeric(tsamp)
            ii = ii+1
          }
                        # print(paste0(month.name[seas], 'day ', j))
                        # if(is.na(get.sst.mask.val(t1[1], t1[2], sstmat, seas))){
            # t1 = SatTagSim::simm.kf(2, u = c(-1*u[1], u[2]), v = c(-1*v[1], v[2]), D = c(D[1], 1000), msp)[2,]

                        ### Miminum distance

                        # tsamp2 = as.data.frame(rbind(c(1,1,2001,t1), c(1,1+uvmult,2001, as.numeric(tsamp))))
                        # print(tsamp2)
                        # newpar = get.uv(tsamp2)
                        # print(newpar)
                        # t1 = SatTagSim::simm.kf(2, u = c(newpar[1], 0), v = c(newpar[2], 0), D, tsamp2[1,4:5])[2,]
                        # print(t1)

                        #              tsamp = t(sapply((rep(2, 0+ii)),function(x) SatTagSim::simm.kf(x, u, v, D, msp)[2,]))
                        #             tval = apply(tsamp, 1, function(y) get.sst.mask.val(y[1], y[2], sstmat, seas))
                        #             tidx = which(tval>=sstol)
                        #
                        #               if(length(tidx) == 0){
                        #                 tsamp = t(sapply((rep(2, 0+ii)),function(x) SatTagSim::simm.kf(x, u = c(-1*u[1], u[2]), v = c(-1*v[1], v[2]), D = c(D[1], 1000), msp)[2,]))
                        #                           tval = apply(tsamp, 1, function(y) get.sst.mask.val(y[1], y[2], sstmat, seas))
                        #                           tidx = which(tval>=sstol)
                        #               }

                        #   							ttmp = cbind(rep(t1[1], 100), rep(t1[2], 100))
                        # 						    tsamp = t(apply(ttmp, 1, function(x)
                        # 						          simm.kf(n = 2, sp = as.numeric(x), u = u, v = v, D = D)[2,]))
                        #   							tval = apply(tsamp, 1, function(y) get.sst.mask.val(y[1], y[2], sstmat))
                        #   							tidx = which(tval>=3)
                        # 						    t1 = as.numeric(tsamp[sample(tidx, 1, prob = tval[tidx]),])
                        # break
                        # ii = ii+100

                        # t1 = as.numeric(tsamp[sample(tidx, 1, prob = tval[tidx]),])
                        # t1 = as.numeric(tsamp)
        }

        if(!is.null(bath)){
          while(.get.bath(t1[1], t1[2], bath)>0){
            t1 = simm.kf(2, u, v, D, msp, ulim, vlim, Dlim)[2,]
          }
        }
        # if(ii>5) break
        # }
        #             points(t1[1], t1[2], pch = 3, cex=.2, col=2)
        #             ii = ii+100
        # 					 t1 = get.nn.mask(t1, sst[[seas]], tol = 5)
        # }
        temp = rbind(temp,  cbind(t(t1), seas))
        msp = ifelse(is.na(t1), msp, t1)
        # lines(temp[,1:2], col = 2)
      }
    }
    # }
    tsim = as.data.frame(temp)
    names(tsim) = c('lon,','lat','Month')
    # lines(tsim[,1:2])
    tsim
  }