Uncertainty map animation

thingo.R

# little animation to illustrate uncertainty in a density map

library(dsm)
library(mvtnorm)
library(ggplot2)
library(animation)
library(viridis)
library(cowplot)

set.seed(1997)


# load the models and prediction data
# using data from the Duke spatial distance sampling course
# http://distancesampling.org/workshops/duke-spatial-2015/practicals.html

# models of sperm whales off the east coast of the US
load("sperm-data/count-models.RData")
load("sperm-data/predgrid.RData")


# ignoring detection function uncertainty for now...

# get the Lp matrix
Lp <- predict(dsm.nb.xy, newdata=predgrid, type="lpmatrix")

# how many realisations do we want?
frames <- 100

# generate the betas from the GAM "posterior"
betas <- rmvnorm(frames, coef(dsm.nb.xy), vcov(dsm.nb.xy))

library(maptools)
library(rgeos)
library(maps)
library(mapdata)
coastline <- map("world", c("USA"), fill=TRUE, plot=FALSE)
coastline <- map2SpatialPolygons(coastline, ID=coastline$name,# ID=IDs,
                                 proj4string=CRS("+proj=longlat +datum=WGS84"))
# AEA projection, as in EEZ models of Roberts et al
Projection <- CRS("+proj=aea +lat_1=38 +lat_2=30 +lat_0=34 +lon_0=-73 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
# project!
coastline <- spTransform(coastline, Projection)
bnd <- as.data.frame(coastline@polygons[[1]]@Polygons[[1]]@coords)
colnames(bnd) <- c("x", "y")

us <- geom_polygon(aes(x=x, y=y), fill = "#1A9850",
               data=bnd)

mma <- 0
mmi <- 0


# use a function to get animation to play nice...
anim_map <- function(){
  # loop to make plots
  for(frame in 1:frames){

    # make the prediction
    predgrid$preds <- predgrid$off.set * exp(Lp%*%betas[frame,])

mma <<- max(mma, predgrid$preds)
mmi <<- min(mmi, predgrid$preds)

    # plot it (using viridis)
    p <- ggplot(predgrid) +
          geom_tile(aes(x=x,y=y,fill=preds), width=sqrt(1e8), height=sqrt(1e8)) +
          us + coord_equal(xlim=range(predgrid$x), ylim=range(predgrid$y)) +
          scale_fill_viridis(limits=c(0,8.6), option="cividis") +
          theme_cowplot() +
          theme(panel.grid.major = element_line(colour="grey92", size=0.25),
                axis.line=element_blank())
    print(p)
  }
}

# make the animation!
ani.options(interval = 0.15)
saveGIF(anim_map(), outdir = "new")