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