Instantly share code, notes, and snippets.

# dill/distanim.R

Last active November 25, 2023 06:20
Show Gist options
• Save dill/2b0c120d5484d338d8ef to your computer and use it in GitHub Desktop.
Animation of a simple distance sampling survey.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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())
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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())