Skip to content

Instantly share code, notes, and snippets.

# dill/distanim.R Last active Feb 2, 2017

Animation of a simple distance sampling survey.
 # 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-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, obs), y = c(obs, obs),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,obs),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())
 # 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[], x=whale1[,1], y=whale1[,2], ysize=0.075, alpha=1) # group 2 add_phylopic_base(whale[], x=whale2[1,1], y=whale2[1,2], ysize=0.05, alpha=0.8) add_phylopic_base(whale[], x=whale2[2,1], y=whale2[2,2], ysize=0.05, alpha=0.8) add_phylopic_base(whale[], 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())
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.