Skip to content

Instantly share code, notes, and snippets.

@dill
Last active November 25, 2023 06:20
Show Gist options
  • Save dill/2b0c120d5484d338d8ef to your computer and use it in GitHub Desktop.
Save dill/2b0c120d5484d338d8ef to your computer and use it in GitHub Desktop.
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[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())
# 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())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment