Last active
May 2, 2022 18:51
-
-
Save SwampThingPaul/cc582ac9f277ba6a34b24532176b1084 to your computer and use it in GitHub Desktop.
GAM animation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment