GAM animation

SpatialAnimation.R

# This is a very general script to perform gif animations from 
# space-time Generalized Additive Model predictions

# GIS libraries 
library(rgdal)
library(rgeos)
library(raster)
library(tmap)

# GAM libraries
library(mgcv)
library(gratia)
library(DHARMa)

# Definate a geographic projection
nad83.pro=CRS("+init=epsg:4269")

# -------------------------------------------------------------------------
## spatial polygon data...added for example purposes
# shore2.sim=readOGR(...)
# CRE.nnc=readOGR(...)

# -------------------------------------------------------------------------
## Some big GAM model with time and space (lat/long) interactions
# m1=bam(...)

# get extent of area of interest
# in this example lat/longs are in decimal degrees
reg.ext=extent(AOI); 

# use expand.grid() to make a data.frame for predicitions
pdat.sp=data.frame(expand.grid(
  DD_long=seq(reg.ext[1],reg.ext[2],by=0.0025),
  DD_lat=seq(reg.ext[3],reg.ext[4],by=0.0025),
  DOY=seq(1,365,20),
  CY=2010
))

fit.mod <- predict(m1, pdat.sp)
pred.mod <- cbind(pdat.sp, Fitted = fit.mod)
# -------------------------------------------------------------------------
# make a series of raster files from GAM predictions

time.vals=expand.grid(
  DOY=seq(1,365,20),
  CY=2010)
  
for(i in 1:nrow(time.vals)){
  # subset pred.mod for each time step
  tmp=subset(pred.mod,CY==time.vals$CY[i]&DOY==time.vals$DOY[i])[,c("Fitted","DD_long","DD_lat")]
  coordinates(tmp)<-~DD_long + DD_lat # define coordinates
  gridded(tmp)<-TRUE # need this to make it a raster
  tmp=as(tmp,"RasterLayer")
  proj4string(tmp)<-nad83.pro
  tmp.m=raster::mask(tmp,CRE.nnc) # clipped raster to specific area in this case CRE.nnc
  
  # assign name of tmp.m ... I've heard pros/cons on this
  # use at your own risk
  assign(paste0("GAM.Kd.",time.vals$CY[i],"_",time.vals$DOY[i]),tmp.m)
  # print interation to track progress
  print(i)
}

#a series of spatial GAM predictions 
# make it a raster stack for tmap to do the animation
# GAM.stack=stack(GAM.Kd.2010_1,
                  GAM.Kd.2010_21,
                  ...)
# stack all your GAM rasters here
                  
# -------------------------------------------------------------------------
# Set tmap to plot
tmap_mode("plot")

# Define your bounding box
bbox.lims=bbox(CRE.nnc)

# tmap (like ggplot) map
GAM.kd.map=tm_shape(shore2.sim,bbox=bbox.lims)+tm_fill(col="cornsilk")+
  tm_shape(GAM.kd.stack2)+
  tm_raster(title="",palette="-viridis",
            breaks=c(0,0.25,0.5,0.75,1,1.5,2,3),
            labels=c("< 0.25","0.25 - 0.50","0.50 - 0.75","0.75 - 1.00","1.00 - 1.50","1.50 - 2.00",">2.00"))+
  tm_shape(shore2.sim)+tm_fill(col="cornsilk")+ # add shore twice for looks only.
  tm_borders(col="grey30",lwd=0.1)+
  tm_facets(free.scales=FALSE,nrow=1,ncol=1)+
  tm_layout(panel.labels=format(as.Date(seq(1,365,20),origin="2010-01-01"),"%b %d %Y"),
            fontfamily = "serif",bg.color="lightblue", 
            panel.label.size=0.8,legend.title.size=0.8,legend.text.size=0.6)+
  tm_legend(title="Light Attenuation Kd (m\u207B\u00B9)",legend.outside=T)

# This is what does the animation part. 
tmap_animation(GAM.kd.map,filename="./Plots/Kd_GAM.gif",delay=100,width=550,height=250,loop=T,dpi=150)



# -------------------------------------------------------------------------
# You can also use the gifski package if you have a series of plots 
# exported to a fold (temporary or "permanent")
# if you " lift the hood" of tmap_animation(...) it uses gifski similar
# to what is below

# a list of file paths to already generated maps/plots
map.files=paste0("./Plots/tmp/GAM.Kd.",time.vals$CY,"_",time.vals$DOY,".png")

gifski::gifski(map.files, 
               delay = 100 / 100, 
               gif_file  = paste0(plot.path,"Plots/Kd_GAM2.gif"),
               loop = T)

## END
## Varsågod