Animation of a simple distance sampling survey.

distanim.R

# animate a (very boring) distance sampling survey
# watch a gif at: http://converged.yt/distance-animation.gif
# written by David Lawrence Miller, with tweaks from Eric Rexstad
# released under the MIT license

library(animation)

# make sure that we get the same result every time...
set.seed(1234)

# need to wrap everything in a function to pass to saveGIF, below...
plotter<-function(){

  # how many animals are there to detect?
  n <- 250
  # generate animals unformly over (0,1)x(0,1)
  xy.dat <- data.frame(x=runif(n),
                       y=runif(n))

  # sort the points by y (bottom to top)
  # this is the direction in which we go along the transect
  xy.dat <- xy.dat[order(xy.dat$y),]

  # give each point an id
  xy.dat$id <- 1:n

  # list of transects
  n.transects <- 10
  transects <- lapply(as.list(seq(0.05,0.95,len=n.transects)),
                      function(x){
                        list(x=c(x,x),
                             y=c(0,1))
                      })

  # vector to hold indices of detected points
  detected <- rep(FALSE,nrow(xy.dat))

  # vector to hold histogram of distances
  histdat <- c(0)

  # side-by-side survey area and histogram
  par(mfrow=c(1,2))

  # loop over transects
  for(tr.id in 1:n.transects){

    # pick out the current transect
    this.tr <- transects[[tr.id]]

    # calculate all distances to this transect first
    ddist <- abs(this.tr$x[1]-xy.dat$x)

    # loop over the possibly observed animals
    for(j in seq_along(xy.dat[,1])){

      # take this observation
      obs <- xy.dat[j,]
      # was the animal detected?
      #detected[obs$id] <- TRUE # strip transects
      # half-nformal detection function with sigma=0.04
      # do rejection sampling
      sigma <- 0.02
      this.detected <- (exp(-ddist[j]^2/(2*sigma^2)) > runif(1))
      # maybe it was already detected? don't un-detect it!
      detected[obs$id] <- detected[obs$id] | this.detected

      # only do this animation stuff when we detect something
      if(this.detected){
        # plot whole study area and all "animals"
        plot(xy.dat$x, xy.dat$y, axes=FALSE, asp=1,
             xlab="", ylab="",
             main="Survey area",
             xlim=c(-0.01,1.01), ylim=c(-0.01,1.01),
             pch=20,cex=0.7,col=grey(0.7))

        # plot transect centre lines for transects already covered
        lapply(transects[1:tr.id],lines, lty=2, col="grey")

        if(!is.null(detected)){
          # plot the animals that have been detected so far in blue
          points(xy.dat[detected,c("x","y")], pch=20,cex=0.7,col="blue")
        }

        # draw perpendicular line -- from transect to animal
        lines(x = c(this.tr$x[1], obs[1]),
              y = c(obs[2], obs[2]),col="red")
        # add this observation to the histogram
        histdat <- c(histdat,ddist[j])

        # draw "along transect" line showing how far we come so far
        lines(x=this.tr$x,y=c(this.tr$y[1],obs[2]),col="red")

        ## update the histogram plot
        xmax <- 3.5 * sigma
        hist(histdat,ylim=c(0,35),xlim=c(0,xmax),
             main="Histogram of observed distances",xlab="Distance",
             breaks=seq(0,xmax,len=10))
#        Sys.sleep(0.5)
        ani.pause()
      } # end if detected stuff
    } # end of current transect loop
  } # end loop over all transects
} # end plotter function

# test it out here, need to replace ani.pause with Sys.sleep(0.5)
#plotter()


saveGIF(plotter(),"distance-animation.gif", interval = 0.2, ani.width = 800,ani.height = 400, outdir=getwd())

mrds-anim.R

# Mark-recapture distance sampling animation
#
# David L Miller 2015
# Released under MIT license
# this is parameterised in terms of "whales" but you can imagine any animal

# watch a gif at: http://converged.yt/mrds-animation.gif

library(animation)

# library to get the whale picture
#library(rphylopic)
library(devtools)
load_all("~/current/rphylopic")

# get a (sperm) whale! (this is "Physeter catodon" by Noah Schlottman)
# http://phylopic.org/image/dc76cbdb-dba5-4d8f-8cf3-809515c30dbd/
whale <- image_data("dc76cbdb-dba5-4d8f-8cf3-809515c30dbd", size=512)

# make sure that we get the same result every time...
set.seed(1234)

t.scale <- 0.01

# whale group locations
whale1 <- matrix(c(0.1, 0.4), 1, 2)

whale2 <- cbind(c(0.75, 0.73, 0.71),
                c(0.78, 0.74, 0.76))

whale2_centroid <- colMeans(whale2)

# move the boats, whales stay stationary
boat_anim <- function(from, to){
  for(i in seq(from, to, t.scale)){

    # survey area
    #plot(0:1, 0:1, asp=1, type="n", axes=FALSE, ylab="", xlab="")
    plot(0:1, 0:1, asp=1, type="n", ylab="", xlab="")

    # transect
    abline(v=0.5, lty=2)

    # a boat
    points(x=0.5, y=i, cex=2, pch=2)

    # some whales
    # group 1
    add_phylopic_base(whale[[1]], x=whale1[,1], y=whale1[,2], ysize=0.075, alpha=1)
    # group 2
    add_phylopic_base(whale[[1]], x=whale2[1,1], y=whale2[1,2], ysize=0.05, alpha=0.8)
    add_phylopic_base(whale[[1]], x=whale2[2,1], y=whale2[2,2], ysize=0.05, alpha=0.8)
    add_phylopic_base(whale[[1]], x=whale2[3,1], y=whale2[3,2], ysize=0.05, alpha=0.8)
  }
}

observe_anim <- function(boat, whale, obs1, obs2){
  # boat   - boat location
  # whale  - whale location
  # obs1/2 - did observer 1/2 see it?

  # if there is more than 1 whale, get centroid
  if(nrow(whale)>1){
    whale <- matrix(colMeans(whale),1,2)
  }

  boatwhale <- rbind(boat, whale)

  # lines to observers
  if(obs1)
    lines(boatwhale+matrix(c(0,0,0.01,0.01),2,2), col="blue")
  if(obs2)
    lines(boatwhale-matrix(c(0,0,0.01,0.01),2,2), col="red")

  # labels
  if(obs1){
    text(matrix(colMeans(boatwhale)+c(0,0.1), 1, 2),
         label="Observer 1: detection", col="blue")
  }else{
    text(matrix(colMeans(boatwhale)+c(0,0.1), 1, 2),
         label="Observer 1: no detection", col="blue")
  }
  if(obs2){
    text(matrix(colMeans(boatwhale)-c(0,0.1), 1, 2),
         label="Observer 2: detection", col="red")
  }else{
    text(matrix(colMeans(boatwhale)-c(0,0.1), 1, 2),
         label="Observer 2: no detection", col="red")
  }
}

### actual animation instructions
plotter <- function(){

  # run the boat to first observation
  boat_anim(0,0.3)

  # observe the whales in group 1
  for(i in 1:30){
    boat_anim(0.3,0.3)
    observe_anim(matrix(c(0.5, 0.3),1,2), whale1, TRUE, FALSE)
    ani.pause()
  }


  # run the boat onward to second observation
  boat_anim(0.3, 0.77)

  # observe the whales in group 1
  for(i in 1:30){
    boat_anim(0.77,0.77)
    observe_anim(matrix(c(0.5, 0.77),1,2), whale2, TRUE, TRUE)
    ani.pause()
  }
  # and on!
  boat_anim(0.77, 1)

} # end plotter function

# make a gif
saveGIF(plotter(), "mrds-animation.gif", ani.width = 500, interval=0.05,
        ani.height = 500, outdir=getwd())