| # 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