# Animations in R: Time series of erythemal irradiance in the St. Lawrence

## The “need” for animations

As part of my PhD thesis, I am currently characterizing the intensity of multiple stressors in the estuary and gulf of St. Lawrence (see ResearchGate project for more details). I have recently needed (read: thought it would be cool) to create an animation of the temporal variations in ultra-violet intensity in the St. Lawrence. Here is how I did it.

## Setting up R

sessionInfo()[[1]]$version.string #R> [1] "R version 3.5.1 (2018-12-12)" ### Installing imagemagick The package animation, which I used to create the animations in R, works with the application imagemagick, which can be installed using Homebrew on MacOS. brew install imagemagick ### Installing the package animation install.packages('animation', repos = 'http://yihui.name/xran') ## Data There is a incredible amount of data available on the NASA website. For this post, we downloaded all available data from NASA’s Ozone Monitoring Instrument (OMI) aboard the Earth Observing System’s (EOS) Aura satellite. More specifically, we used OMI/Aura Surface UVB Irradiance and Erythemal Dose Daily at a scale of 0.25 x 0.25 degree and selected the daily erythemal irradiance (mW/m2), i.e. the potentially harmful range of UV radiations, from October 1st 2004 to April 5th 2017. We then selected only the years with 12 months worth of data to create the following figures, i.e. from January 1st 2005 to December 31st 2016. The maps generated in this post require the user to load ocean shapefile from Natural Earth. # UV data UV <- readRDS('MonthOMI.rds') # Projection spatProj <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs' # Ocean shapefile ocean <- readOGR(dsn = "shpfile", layer = "ne_10m_ocean") %>% spTransform(CRSobj = CRS(spatProj)) #R> OGR data source with driver: ESRI Shapefile #R> Source: "/home/travis/build/inSileco/inSileco.github.io/content/post/animationinR/shpfile", layer: "ne_10m_ocean" #R> with 1 features #R> It has 2 fields class(ocean) #R> [1] "SpatialPolygonsDataFrame" #R> attr(,"package") #R> [1] "sp" ## Set parameters and functions We start be setting the different parameters and functions required to build the animation. # Graphical parameters min_ <- min(UV[[1]], na.rm = T) %>% plyr::round_any(., 100, floor) # Minimum values in UV dataset max_ <- max(UV[[1]], na.rm = T) %>% plyr::round_any(., 100, ceiling) # Minimum values in UV dataset ext <- c(-72, -54, 44, 53) # c(xmin, xmax, ymin, ymax) rbPal <- colorRampPalette(c('#2f6eb9','#2aadba','#c4c726','#b45f5f')) # color palette labels <- matrix(seq(4,(151-4)), ncol = 12, nrow = 12, dimnames = list(2005:2016, month.abb), byrow = T) # Select only full years (start: 2004/10; end: 2017/04) # Monthly trends monthMean <- apply(X = UV[[1]], MARGIN = 3, FUN = mean, na.rm = T) # monthly mean monthSd <- apply(X = UV[[1]], MARGIN = 3, FUN = sd, na.rm = T) # monthly sd years <- 2005:2016 lowCI <- monthMean - (1.96*monthSd) # Lower confidence interval upCI <- monthMean + (1.96*monthSd) # Upper confidence interval # Empty plot function eplot <- function(x = 1, y = 1, xmin = 0, xmax = 1, ymin = 0, ymax = 1) { plot(x = x,y = y,bty = "n",ann = FALSE,xaxt = "n",yaxt = "n",type = "n",bg = 'grey',ylim = c(ymin,ymax),xlim = c(xmin,xmax)) } # Legend function colorBar <- function(colRamp, min, max, nticks=11, ticks=seq(min, max, len=nticks)) { scale = (length(colRamp)-1)/(max-min) par(mar = c(5,0.5,5,4)) plot(c(0,10), c(min_,max_), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='') axis(4, ticks, las=1, cex = 1.5) for (i in 1:(length(colRamp)-1)) { y = (i-1)/scale + min rect(0,y,10,y+1/scale, col=colRamp[i], border=NA) } } ## Create first animation Let’s start with a simple animation of erythemal irradiance for 2005. saveGIF({ for(i in 1:ncol(labels)){ par(mar = c(0,0,0,0), pch = 21, xaxs = "i", yaxs = "i", family = "serif") r <- raster::raster(UV[[1]][,,labels[1,i]], xmn = ext[1], xmx = ext[2], ymn = ext[3], ymx = ext[4], crs = spatProj) %>% mask(ocean) eplot(xmin = ext[1], xmax = ext[2], ymin = ext[3], ymax = ext[4]) plot(r, col= rbPal(100), breaks = seq(min_,max_,length.out=100), legend = F, add = T) plot(ocean, add = T, lwd = 2) text(x = -54.5, y = 45, labels = paste(colnames(labels)[1], rownames(labels)[i]), adj = 1, cex = 1.5) } }, movie.name = "UV2005animate.gif") ## Create complex animation Now let’s make a more informative and complex animation that allows to visualize monthly erythemal irradiance through the years as well as annual trends. ani.options(ani.height = 720, ani.width = 720) saveGIF({ for(i in 1:nrow(labels)) { #layout mat <- matrix(c(13:15,1:3,16:18,4:6,19:21,7:9,22:24,10:12), ncol = 3, nrow = 8, byrow = T) %>% cbind(c(0,0,25,25,25,25,0,0)) %>% rbind(26, .) layout(mat, widths = c(3,3,3,1.25), heights = c(1,rep(c(2,3),4))) # Rasters for(j in 1:ncol(labels)) { par(pch = 21, xaxs = "i", yaxs = "i", family = "serif", mar = c(0,0,0,0)) r <- raster::raster(UV[[1]][,,labels[i,j]], xmn = ext[1], xmx = ext[2], ymn = ext[3], ymx = ext[4], crs = spatProj) %>% mask(ocean) eplot(xmin = ext[1], xmax = ext[2], ymin = ext[3], ymax = ext[4]) image(r, col= rbPal(100), breaks = seq(0,max_,length.out=101), legend = F, add = T) plot(ocean, add = T) text(x = -54.5, y = 45, labels = paste(colnames(labels)[j], rownames(labels)[i]), adj = 1, cex = 1.5) } # Plot trends for(j in 1:ncol(labels)) { temporal <- data.frame(years = years[1:i], monthMean = monthMean[labels[1:i,j]], lowCI = lowCI[labels[1:i,j]], upCI = upCI[labels[1:i,j]]) %>% na.omit() polyY <- c(temporal$upCI, rev(temporal$lowCI)) # Y coordinates for confidence interval polygon polyX <- c(temporal$years, rev(temporal$years)) # X coordinates for confidence interval polygon par(mar = c(0,3,1,1), family = "serif") eplot(xmin = 2005, xmax = 2016, ymin = 0, ymax = max(upCI)) polygon(polyX, polyY, border = 'transparent', col = "#d9e5f4") lines(temporal$monthMean ~ temporal\$years, lwd = 0.5, col = '#0057bb')
axis(side = 2, at = seq(0, plyr::round_any(max(upCI), 1000, ceiling), by = 2500), las = 1, pos = 2004.75, cex =  1.5)
}

# Plot legend
colorBar(rbPal(101), min = 0, max = max_)

# Graph title
par(mar = c(0,0,0,0))
eplot(xmin = 0, xmax = 10, ymin = 0, ymax = 10)
text(x = 5, y = 5, labels = bquote("Erythemal UV Irradiance (mW/m" ^2 *")"), cex = 2)

}
}, movie.name = "UVanimate.gif")

## Concluding remarks

As you can see, it is quite straightforward to create powerful and informative animations in R.