Time Series - Displaying time series, spatial and space-time data with R
Table of Contents
1. Time on the horizontal axis
1.1. Time Graph of Variables with Different Scales
1.2. Time Series of Variables with the Same Scale
1.3. The Horizon Graph
1.4. Stacked Graphs
1.5. Interactive
1.5.1. dygraphs
1.5.2. highcharter
1.5.3. plotly
1.5.4. streamgraph
2. Time as a conditioning or grouping variable
2.1. Scatterplot Matrix: Time as a Grouping Variable
2.2. Scatterplot with Time as a Conditioning Variable
3. Time as a complementary variable
3.1. Polylines
3.2. Labels to Show Time Information
3.3. Using Variable Size to Encode an Additional Variable
3.4. Interactive graphics: animation
3.4.1. plotly
3.4.2. gganimate
3.4.3. gridSVG
4. Code
4.1. Time on the horizontal axis
################################################################## ## Initial configuration ################################################################## ## Clone or download the repository and set the working directory ## with setwd to the folder where the repository is located. library("lattice") library("ggplot2") ## latticeExtra must be loaded after ggplot2 to prevent masking of its ## `layer` function. library("latticeExtra") library("RColorBrewer") source('configLattice.R') ################################################################## ################################################################## ## Time graph of variables with different scales ################################################################## library("zoo") load('data/TimeSeries/aranjuez.RData') ## The layout argument arranges panels in rows xyplot(aranjuez, layout = c(1, ncol(aranjuez))) autoplot(aranjuez) + facet_free() ################################################################## ## Annotations to enhance the time graph ################################################################## ## lattice version ## Auxiliary function to extract the year value of a POSIXct time ## index Year <- function(x)format(x, "%Y") xyplot(aranjuez, layout = c(1, ncol(aranjuez)), strip = FALSE, scales = list(y = list(cex = 0.6, rot = 0)), panel = function(x, y, ...){ ## Alternation of years panel.xblocks(x, Year, col = c("lightgray", "white"), border = "darkgray") ## Values under the average highlighted with red regions panel.xblocks(x, y < mean(y, na.rm = TRUE), col = "indianred1", height = unit(0.1, 'npc')) ## Time series panel.lines(x, y, col = 'royalblue4', lwd = 0.5, ...) ## Label of each time series panel.text(x[1], min(y, na.rm = TRUE), names(aranjuez)[panel.number()], cex = 0.6, adj = c(0, 0), srt = 90, ...) ## Triangles to point the maxima and minima idxMax <- which.max(y) panel.points(x[idxMax], y[idxMax], col = 'black', fill = 'lightblue', pch = 24) idxMin <- which.min(y) panel.points(x[idxMin], y[idxMin], col = 'black', fill = 'lightblue', pch = 25) }) ## ggplot2 version timeIdx <- index(aranjuez) aranjuezLong <- fortify(aranjuez, melt = TRUE) summary(aranjuezLong) ## Values below mean are negative after being centered scaled <- fortify(scale(aranjuez, scale = FALSE), melt = TRUE) ## The 'scaled' column is the result of the centering. ## The new 'Value' column store the original values. scaled <- transform(scaled, scaled = Value, Value = aranjuezLong$Value) underIdx <- which(scaled$scaled <= 0) ## 'under' is the subset of values below the average under <- scaled[underIdx,] library("xts") ep <- endpoints(timeIdx, on = 'years') ep <- ep[-1] N <- length(ep) ## 'tsp' is start and 'tep' is the end of each band. One of each two ## years are selected. tep <- timeIdx[ep[seq(1, N, 2)] + 1] tsp <- timeIdx[ep[seq(2, N, 2)]] minIdx <- timeIdx[apply(aranjuez, 2, which.min)] minVals <- apply(aranjuez, 2, min, na.rm = TRUE) mins <- data.frame(Index = minIdx, Value = minVals, Series = names(aranjuez)) maxIdx <- timeIdx[apply(aranjuez, 2, which.max)] maxVals <- apply(aranjuez, 2, max, na.rm = TRUE) maxs <- data.frame(Index = maxIdx, Value = maxVals, Series = names(aranjuez)) ggplot(data = aranjuezLong, aes(Index, Value)) + ## Time series of each variable geom_line(colour = "royalblue4", lwd = 0.5) + ## Year bands annotate('rect', xmin = tsp, xmax = tep, ymin = -Inf, ymax = Inf, alpha = 0.4) + ## Values below average geom_rug(data = under, sides = 'b', col = 'indianred1') + ## Minima geom_point(data = mins, pch = 25) + ## Maxima geom_point(data = maxs, pch = 24) + ## Axis labels and theme definition labs(x = 'Time', y = NULL) + theme_bw() + ## Each series is displayed in a different panel with an ## independent y scale facet_free() ################################################################## ## Time series of variables with the same scale ################################################################## load('data/TimeSeries/navarra.RData') avRad <- zoo(rowMeans(navarra, na.rm = 1), index(navarra)) pNavarra <- xyplot(navarra - avRad, superpose = TRUE, auto.key = FALSE, lwd = 0.5, alpha = 0.3, col = 'midnightblue') pNavarra ################################################################## ## Aspect ratio and rate of change ################################################################## xyplot(navarra - avRad, aspect = 'xy', cut = list(n = 3, overlap = 0.1), strip = FALSE, superpose = TRUE, auto.key = FALSE, lwd = 0.5, alpha = 0.3, col = 'midnightblue') ################################################################## ## The horizon graph ################################################################## horizonplot(navarra - avRad, layout = c(1, ncol(navarra)), origin = 0, ## Deviations in each panel are calculated ## from this value colorkey = TRUE, col.regions = brewer.pal(6, "RdBu")) ################################################################## ## Time graph of the differences between a time series and a reference ################################################################## Ta <- aranjuez$TempAvg timeIndex <- index(aranjuez) longTa <- ave(Ta, format(timeIndex, '%j')) diffTa <- (Ta - longTa) xyplot(cbind(Ta, longTa, diffTa), col = c('darkgray', 'red', 'midnightblue'), superpose = TRUE, auto.key = list(space = 'right'), screens = c(rep('Average Temperature', 2), 'Differences')) years <- unique(format(timeIndex, '%Y')) horizonplot(diffTa, cut = list(n = 8, overlap = 0), colorkey = TRUE, col.regions = brewer.pal(6, "RdBu"), layout = c(1, 8), scales = list(draw = FALSE, y = list(relation = 'same')), origin = 0, strip.left = FALSE) + layer(grid.text(years[panel.number()], x = 0, y = 0.1, gp = gpar(cex = 0.8), just = "left")) year <- function(x)as.numeric(format(x, '%Y')) day <- function(x)as.numeric(format(x, '%d')) month <- function(x)as.numeric(format(x, '%m')) myTheme <- modifyList(custom.theme(region = brewer.pal(9, 'RdBu')), list( strip.background = list(col = 'gray'), panel.background = list(col = 'gray'))) maxZ <- max(abs(diffTa)) levelplot(diffTa ~ day(timeIndex) * year(timeIndex) | factor(month(timeIndex)), at = pretty(c(-maxZ, maxZ), n = 8), colorkey = list(height = 0.3), layout = c(1, 12), strip = FALSE, strip.left = TRUE, xlab = 'Day', ylab = 'Month', par.settings = myTheme) df <- data.frame(Vals = diffTa, Day = day(timeIndex), Year = year(timeIndex), Month = month(timeIndex)) library("scales") ## The packages scales is needed for the breaks_pretty function. ggplot(data = df, aes(fill = Vals, x = Day, y = Year)) + facet_wrap(~ Month, ncol = 1, strip.position = 'left') + scale_y_continuous(breaks = breaks_pretty()) + scale_fill_distiller(palette = 'RdBu', direction = 1) + geom_raster() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) ################################################################## ## Stacked graphs ################################################################## load('data/TimeSeries/unemployUSA.RData') xyplot(unemployUSA, superpose = TRUE, par.settings = custom.theme, auto.key = list(space = 'right')) autoplot(unemployUSA, facets = NULL) + geom_area(aes(fill = Series)) + scale_x_yearmon() library("ggstream") sep2008 <- as.numeric(as.yearmon('2008-09')) autoplot(unemployUSA, facets = NULL) + geom_stream(aes(fill = Series), color = "black", lwd = 0.25, bw = 0.5) + scale_x_yearmon() + geom_vline(xintercept = sep2008, lwd = 0.25, color = "gray50") + theme_bw() library("grid") panel.flow <- function(x, y, groups, origin, ...) { dat <- data.frame(x = x, y = y, groups = groups) (dataframe) nVars <- nlevels(groups) groupLevels <- levels(groups) ## From long to wide yWide <- unstack(dat, y~groups) (yunstack) ## Where are the maxima of each variable located? We will use ## them to position labels. idxMaxes <- apply(yWide, 2, which.max) ##Origin calculated following Havr.eHetzler.ea2002 if (origin=='themeRiver') origin = -1/2*rowSums(yWide) else origin = 0 yWide <- cbind(origin = origin, yWide) ## Cumulative sums to define the polygon yCumSum <- t(apply(yWide, 1, cumsum)) (ycumsum) Y <- as.data.frame(sapply(seq_len(nVars), (yPolygon) function(iCol)c(yCumSum[,iCol+1], rev(yCumSum[,iCol])))) names(Y) <- levels(groups) ## Back to long format, since xyplot works that way y <- stack(Y)$values (yStack) ## Similar but easier for x xWide <- unstack(dat, x~groups) (xunstack) x <- rep(c(xWide[,1], rev(xWide[,1])), nVars) (xPolygon) ## Groups repeated twice (upper and lower limits of the polygon) groups <- rep(groups, each = 2) (groups) ## Graphical parameters superpose.polygon <- trellis.par.get("superpose.polygon") (gpar) col = superpose.polygon$col border = superpose.polygon$border lwd = superpose.polygon$lwd ## Draw polygons for (i in seq_len(nVars)){ (forVars) xi <- x[groups==groupLevels[i]] yi <- y[groups==groupLevels[i]] panel.polygon(xi, yi, border = border, lwd = lwd, col = col[i]) } ## Print labels for (i in seq_len(nVars)){ xi <- x[groups==groupLevels[i]] yi <- y[groups==groupLevels[i]] N <- length(xi)/2 ## Height available for the label h <- unit(yi[idxMaxes[i]], 'native') - unit(yi[idxMaxes[i] + 2*(N-idxMaxes[i]) +1], 'native') ##...converted to "char" units hChar <- convertHeight(h, 'char', TRUE) ## If there is enough space and we are not at the first or ## last variable, then the label is printed inside the polygon. if((hChar >= 1) && !(i %in% c(1, nVars))){ (hChar) grid.text(groupLevels[i], xi[idxMaxes[i]], (yi[idxMaxes[i]] + yi[idxMaxes[i] + 2*(N-idxMaxes[i]) +1])/2, gp = gpar(col = 'white', alpha = 0.7, cex = 0.7), default.units = 'native') } else { ## Elsewhere, the label is printed outside grid.text(groupLevels[i], xi[N], (yi[N] + yi[N+1])/2, gp = gpar(col = col[i], cex = 0.7), just = 'left', default.units = 'native') } } } prepanel.flow <- function(x, y, groups, origin,...) { dat <- data.frame(x = x, y = y, groups = groups) nVars <- nlevels(groups) groupLevels <- levels(groups) yWide <- unstack(dat, y~groups) if (origin=='themeRiver') origin = -1/2*rowSums(yWide) else origin = 0 yWide <- cbind(origin = origin, yWide) yCumSum <- t(apply(yWide, 1, cumsum)) list(xlim = range(x), (xlim) ylim = c(min(yCumSum[,1]), max(yCumSum[,nVars+1])), (ylim) dx = diff(x), (dx) dy = diff(c(yCumSum[,-1]))) } library("colorspace") ## We will use a qualitative palette from colorspace nCols <- ncol(unemployUSA) pal <- rainbow_hcl(nCols, c = 70, l = 75, start = 30, end = 300) myTheme <- custom.theme(fill = pal, lwd = 0.2) sep2008 <- as.numeric(as.yearmon('2008-09')) xyplot(unemployUSA, superpose = TRUE, auto.key = FALSE, panel = panel.flow, prepanel = prepanel.flow, origin = 'themeRiver', scales = list(y = list(draw = FALSE)), par.settings = myTheme) + layer(panel.abline(v = sep2008, col = 'gray', lwd = 0.7)) ################################################################## ## Panel and prepanel functions to implement the ThemeRiver with =xyplot= ################################################################## library("dygraphs") dyTemp <- dygraph(aranjuez[, c("TempMin", "TempAvg", "TempMax")], main = "Temperature in Aranjuez", ylab = "ºC") dyTemp dyTemp %>% dyHighlight(highlightSeriesBackgroundAlpha = 0.2, highlightSeriesOpts = list(strokeWidth = 2)) dygraph(aranjuez[, c("TempMin", "TempAvg", "TempMax")], main = "Temperature in Aranjuez", ylab = "ºC") %>% dySeries(c("TempMin", "TempAvg", "TempMax"), label = "Temperature") library("highcharter") library("xts") aranjuezXTS <- as.xts(aranjuez) highchart(type = "stock") %>% hc_add_series(name = 'TempMax', aranjuezXTS[, "TempMax"]) %>% hc_add_series(name = 'TempMin', aranjuezXTS[, "TempMin"]) %>% hc_add_series(name = 'TempAvg', aranjuezXTS[, "TempAvg"]) aranjuezDF <- fortify(aranjuez[, c("TempMax", "TempAvg", "TempMin")], melt = TRUE) summary(aranjuezDF) library("plotly") plot_ly(aranjuezDF) %>% add_lines(x = ~ Index, y = ~ Value, color = ~ Series) ## remotes::install_github("hrbrmstr/streamgraph") library("streamgraph") unemployDF <- fortify(unemployUSA, melt = TRUE) ## streamgraph does not work with the yearmon class unemployDF$Index <- as.Date(unemployDF$Index) head(unemployDF) streamgraph(unemployDF, key = "Series", value = "Value", date = "Index") %>% sg_axis_x(1, "year", "%Y") %>% sg_fill_brewer("Set1")
4.2. Time as a conditioning or grouping variable
################################################################## ## Initial configuration ################################################################## ## Clone or download the repository and set the working directory ## with setwd to the folder where the repository is located. library("lattice") library("ggplot2") ## latticeExtra must be loaded after ggplot2 to prevent masking of its ## `layer` function. library("latticeExtra") library("RColorBrewer") source('configLattice.R') ################################################################## ################################################################## ## Scatterplot matrix: time as a grouping variable ################################################################## library("zoo") load('data/TimeSeries/aranjuez.RData') aranjuezDF <- as.data.frame(aranjuez) aranjuezDF$Month <- format(index(aranjuez), '%m') ## Red-Blue palette with black added (12 colors) colors <- c(brewer.pal(n = 11, 'RdBu'), '#000000') ## Rearrange according to months (darkest for summer) colors <- colors[c(6:1, 12:7)] splom(~ aranjuezDF[1:10], groups = aranjuezDF$Month, auto.key = list(space = 'right', title = 'Month', cex.title = 1), pscale = 0, varname.cex = 0.7, xlab = '', par.settings = custom.theme(symbol = colors, pch = 19), cex = 0.3, alpha = 0.1) trellis.focus('panel', 1, 1) idx <- panel.link.splom(pch = 13, cex = 0.6, col = 'green') aranjuez[idx,] library("GGally") ggpairs(aranjuezDF, columns = 1:10, ## Do not include "Month" upper = list(continuous = wrap("points", alpha = 0.7)), lower = list(continuous = wrap("points", alpha = 0.7)), diag = list(continuous = wrap("densityDiag", alpha = 0.7)), mapping = aes(colour = Month), legend = 1) + scale_fill_manual(values = colors) + scale_colour_manual(values = colors) + theme(axis.text = element_text(size = 6), axis.text.x = element_text(angle = 45), strip.text = element_text(size = 7)) ################################################################## ## Hexagonal binning ################################################################## library("hexbin") splom(~as.data.frame(aranjuez), panel = panel.hexbinplot, diag.panel = function(x, ...){ yrng <- current.panel.limits()$ylim d <- density(x, na.rm = TRUE) d$y <- with(d, yrng[1] + 0.95 * diff(yrng) * y / max(y)) panel.lines(d) diag.panel.splom(x, ...) }, lower.panel = function(x, y, ...){ panel.hexbinplot(x, y, ...) panel.loess(x, y, ..., col = 'red') }, xlab = '', pscale = 0, varname.cex = 0.7) library("reshape2") aranjuezRshp <- melt(aranjuezDF, measure.vars = c('TempMax', 'TempAvg', 'TempMin'), variable.name = 'Statistic', value.name = 'Temperature') summary(aranjuezRshp) hexbinplot(Radiation ~ Temperature | Statistic, data = aranjuezRshp, layout = c(1, 3)) + layer(panel.loess(..., col = 'red')) ggplot(data = aranjuezRshp, aes(Temperature, Radiation)) + geom_hex() + stat_smooth(se = FALSE, method = 'loess', col = 'red') + facet_wrap(~ Statistic, ncol = 1) + theme_bw() ################################################################## ## Scatterplot with time as a conditioning variable ################################################################## ggplot(data = aranjuezRshp, aes(Radiation, Temperature)) + facet_grid(Statistic ~ Month) + geom_point(col = 'skyblue4', pch = 19, cex = 0.5, alpha = 0.3) + geom_rug() + stat_smooth(se = FALSE, method = 'loess', col = 'indianred1', lwd = 1.2) + theme_bw() + theme(axis.text.x = element_text(size = 7, angle = 45)) useOuterStrips( xyplot(Temperature ~ Radiation | Month * Statistic, data = aranjuezRshp, between = list(x = 0), col = 'skyblue4', pch = 19, cex = 0.5, alpha = 0.3)) + layer({ panel.rug(..., col.line = 'indianred1', end = 0.05, alpha = 0.6) panel.loess(..., col = 'indianred1', lwd = 1.5, alpha = 1) })
4.3. Time as a complementary variable
################################################################## ## Initial configuration ################################################################## ## Clone or download the repository and set the working directory ## with setwd to the folder where the repository is located. library("lattice") library("ggplot2") ## latticeExtra must be loaded after ggplot2 to prevent masking of its ## `layer` function. library("latticeExtra") library("RColorBrewer") source('configLattice.R') ################################################################## ################################################################## ## Polylines ################################################################## library("zoo") load('data/TimeSeries/CO2.RData') ## lattice version xyplot(GNI.capita ~ CO2.capita, data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, type = 'b') ## ggplot2 version ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, color = Country.Name)) + xlab("Carbon dioxide emissions (metric tons per capita)") + ylab("GNI per capita, PPP (current international $)") + geom_point() + geom_path() + theme_bw() ################################################################## ## Choosing colors ################################################################## nCountries <- nlevels(CO2data$Country.Name) pal <- brewer.pal(n = nCountries, 'Set1') ## Rank of average values of CO2 per capita CO2mean <- aggregate(CO2.capita ~ Country.Name, data = CO2data, FUN = mean) palOrdered <- pal[rank(CO2mean$CO2.capita)] library("reshape2") CO2capita <- CO2data[, c('Country.Name', 'Year', 'CO2.capita')] CO2capita <- dcast(CO2capita, Country.Name ~ Year) summary(CO2capita) hCO2 <- hclust(dist(CO2capita[, -1])) oldpar <- par(mar = c(0, 2, 0, 0) + .1) plot(hCO2, labels = CO2capita$Country.Name, xlab = '', ylab = '', sub = '', main = '') par(oldpar) idx <- match(levels(CO2data$Country.Name), CO2capita$Country.Name[hCO2$order]) palOrdered <- pal[idx] ## simpleTheme encapsulates the palette in a new theme for xyplot myTheme <- simpleTheme(pch = 19, cex = 0.6, col = palOrdered) ## lattice version pCO2.capita <- xyplot(GNI.capita ~ CO2.capita, data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, par.settings = myTheme, type = 'b') ## ggplot2 version gCO2.capita <- ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, color = Country.Name)) + geom_point() + geom_path() + scale_color_manual(values = palOrdered, guide = "none") + xlab('CO2 emissions (metric tons per capita)') + ylab('GNI per capita, PPP (current international $)') + theme_bw() ################################################################## ## Labels to show time information ################################################################## xyplot(GNI.capita ~ CO2.capita, data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, par.settings = myTheme, type = 'b', panel = function(x, y, ..., subscripts, groups){ panel.text(x, y, ..., labels = CO2data$Year[subscripts], pos = 2, cex = 0.5, col = 'gray') panel.superpose(x, y, subscripts, groups,...) }) ## lattice version pCO2.capita <- pCO2.capita + glayer_(panel.text(..., labels = CO2data$Year[subscripts], pos = 2, cex = 0.5, col = 'gray')) ## ggplot2 version gCO2.capita <- gCO2.capita + geom_text(aes(label = Year), colour = 'gray', size = 2.5, hjust = 0, vjust = 0) ################################################################## ## Country names: positioning labels ################################################################## library("directlabels") ## lattice version direct.label(pCO2.capita, method = 'extreme.grid') ## ggplot2 version direct.label(gCO2.capita, method = 'extreme.grid') ################################################################## ## A panel for each year ################################################################## ## lattice version xyplot(GNI.capita ~ CO2.capita | factor(Year), data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, type = 'p', auto.key = list(space = 'right')) ## ggplot2 version ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, colour = Country.Name)) + facet_wrap(~ Year) + geom_point(pch = 19) + xlab('CO2 emissions (metric tons per capita)') + ylab('GNI per capita, PPP (current international $)') + theme_bw() ## lattice version pl <- xyplot(GNI.capita ~ CO2.capita | factor(Year), data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, type = 'b', par.settings = myTheme) ## Extracted from the qp.labels help page dlmethod <- list("first.points", cex = 0.5, "calc.boxes", "enlarge.box", qp.labels("y","bottom","top")) direct.label(pl, dlmethod) ## ggplot2 version pg <- ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, colour = Country.Name)) + facet_wrap(~ Year) + geom_point(pch = 19) + scale_color_manual(values = palOrdered, guide = "none") + xlab('CO2 emissions (metric tons per capita)') + ylab('GNI per capita, PPP (current international $)') + theme_bw() direct.label(pg, dlmethod) ################################################################## ## Using variable size to encode an additional variable ################################################################## library("classInt") CO2data$CO2.PPP <- with(CO2data, CO2/GNI.PPP * 1e9) intervals <- classIntervals(CO2data$CO2.PPP, n = 4, style = 'fisher') nInt <- length(intervals$brks) - 1 cex.key <- seq(0.5, 1.8, length = nInt) idx <- findCols(intervals) CO2data$cexPoints <- cex.key[idx] ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, colour = Country.Name)) + facet_wrap(~ Year) + geom_point(aes(size = cexPoints), pch = 19) + xlab('Carbon dioxide emissions (metric tons per capita)') + ylab('GNI per capita, PPP (current international $)') + theme_bw() op <- options(digits = 2) tab <- print(intervals) options(op) key <- list(space = 'right', title = expression(CO[2]/GNI.PPP), cex.title = 1, ## Labels of the key are the intervals strings text = list(labels = names(tab), cex = 0.85), ## Points sizes are defined with cex.key points = list(col = 'black', pch = 19, cex = cex.key, alpha = 0.7) ) pl2 <- xyplot(GNI.capita ~ CO2.capita|factor(Year), data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", groups = Country.Name, alpha = 0.7, panel = panel.superpose, panel.groups = function(x, y, subscripts, group.number, group.value, ...){ panel.xyplot(x, y, col = palOrdered[group.number], cex = CO2data$cexPoints[subscripts]) }) ## Add labels... plLabel <- direct.label(pl2, "dlmethod") ## ... and the key (after the call to direct.labels because this ## function removes the legend) update(plLabel, key = key) library("plotly") p <- plot_ly(CO2data, x = ~CO2.capita, y = ~GNI.capita, sizes = c(10, 100), marker = list(opacity = 0.7, sizemode = 'diameter')) p <- add_markers(p, size = ~CO2.PPP, color = ~Country.Name, text = ~Country.Name, hoverinfo = "text", ids = ~Country.Name, frame = ~Year, showlegend = FALSE) p <- animation_opts(p, frame = 1000, transition = 800, redraw = FALSE) p <- animation_slider(p, currentvalue = list(prefix = "Year ")) p ################################################################## ## gganimate ################################################################## library("gganimate") pgvis <- ggplot(data = CO2data, aes(x = CO2.capita, y = GNI.capita, size = CO2.PPP, color = Country.Name)) + geom_point(pch = 19) + scale_size(range = c(1, 12)) + scale_color_viridis_d() + theme_bw() + labs(title = 'Year: {trunc(frame_time)}', x = 'CO2 emissions (metric tons per capita)', y = 'GNI per capita, PPP (current international $)') pgAnim <- pgvis + transition_time(Year) + ease_aes("linear") + shadow_wake(wake_length = 0.3) animate(pgAnim, height = 1080, width = 1080, res = 150, units = "px") anim_save("figs/TimeSeries/CO2_gganimate.gif") library("gridSVG") library("grid") xyplot(GNI.capita ~ CO2.capita, data = CO2data, xlab = "Carbon dioxide emissions (metric tons per capita)", ylab = "GNI per capita, PPP (current international $)", subset = Year==2000, groups = Country.Name, ## The limits of the graphic are defined ## with the entire dataset xlim = extendrange(CO2data$CO2.capita), ylim = extendrange(CO2data$GNI.capita), panel = function(x, y, ..., subscripts, groups) { color <- palOrdered[groups[subscripts]] radius <- CO2data$CO2.PPP[subscripts] ## Size of labels cex <- 1.1*sqrt(radius) ## Bubbles grid.circle(x, y, default.units = "native", r = radius*unit(.25, "inch"), name = trellis.grobname("points", type = "panel"), gp = gpar(col = color, ## Fill color ligther than border fill = adjustcolor(color, alpha = .5), lwd = 2)) ## Country labels grid.text(label = groups[subscripts], x = unit(x, 'native'), ## Labels above each bubble y = unit(y, 'native') + 1.5 * radius *unit(.25, 'inch'), name = trellis.grobname('labels', type = 'panel'), gp = gpar(col = color, cex = cex)) }) ## Duration in seconds of the animation duration <- 20 nCountries <- nlevels(CO2data$Country.Name) years <- unique(CO2data$Year) nYears <- length(years) ## Intermediate positions of the bubbles x_points <- animUnit(unit(CO2data$CO2.capita, 'native'), id = rep(seq_len(nCountries), each = nYears)) y_points <- animUnit(unit(CO2data$GNI.capita, 'native'), id = rep(seq_len(nCountries), each = nYears)) ## Intermediate positions of the labels y_labels <- animUnit(unit(CO2data$GNI.capita, 'native') + 1.5 * CO2data$CO2.PPP * unit(.25, 'inch'), id = rep(seq_len(nCountries), each = nYears)) ## Intermediate sizes of the bubbles size <- animUnit(CO2data$CO2.PPP * unit(.25, 'inch'), id = rep(seq_len(nCountries), each = nYears)) grid.animate(trellis.grobname("points", type = "panel", row = 1, col = 1), duration = duration, x = x_points, y = y_points, r = size, rep = TRUE) grid.animate(trellis.grobname("labels", type = "panel", row = 1, col = 1), duration = duration, x = x_points, y = y_labels, rep = TRUE) countries <- unique(CO2data$Country.Name) URL <- paste('http://en.wikipedia.org/wiki/', countries, sep = '') grid.hyperlink(trellis.grobname('points', type = 'panel', row = 1, col = 1), URL, group = FALSE) visibility <- matrix("hidden", nrow = nYears, ncol = nYears) diag(visibility) <- "visible" yearText <- animateGrob(garnishGrob(textGrob(years, .9, .15, name = "year", gp = gpar(cex = 2, col = "grey")), visibility = "hidden"), duration = 20, visibility = visibility, rep = TRUE) grid.draw(yearText) grid.export("figs/TimeSeries/bubbles.svg")
4.4. Data
################################################################## ## SIAR ################################################################## ################################################################## ## Daily data of different meteorological variables ################################################################## library("zoo") aranjuez <- read.zoo("data/TimeSeries/aranjuez.gz", index.column = 3, format = "%d/%m/%Y", fileEncoding = 'UTF-16LE', header = TRUE, fill = TRUE, sep = ';', dec = ",", as.is = TRUE) aranjuez <- aranjuez[, -c(1:4)] names(aranjuez) <- c('TempAvg', 'TempMax', 'TempMin', 'HumidAvg', 'HumidMax', 'WindAvg', 'WindMax', 'Radiation', 'Rain', 'ET') aranjuezClean <- within(as.data.frame(aranjuez),{ TempMin[TempMin > 40] <- NA HumidMax[HumidMax > 100] <- NA }) aranjuez <- zoo(aranjuezClean, index(aranjuez)) save(aranjuez, file = 'data/aranjuez.RData') ################################################################## ## Solar radiation measurements from different locations ################################################################## library("zoo") load('data/TimeSeries/navarra.RData') ################################################################## ## Unemployment in the United States ################################################################## unemployUSA <- read.csv('data/TimeSeries/unemployUSA.csv') nms <- unemployUSA$Series.ID ##columns of annual summaries annualCols <- 14 + 13*(0:12) ## Transpose. Remove annual summaries unemployUSA <- as.data.frame(t(unemployUSA[,-c(1, annualCols)])) ## First 7 characters can be suppressed names(unemployUSA) <- substring(nms, 7) library("zoo") Sys.setlocale("LC_TIME", 'C') idx <- as.yearmon(row.names(unemployUSA), format = '%b.%Y') unemployUSA <- zoo(unemployUSA, idx) unemployUSA <- unemployUSA[complete.cases(unemployUSA), ] save(unemployUSA, file = 'data/unemployUSA.RData') ################################################################## ## Gross National Income and $CO_2$ emissions ################################################################## library("WDI") countries <- c('CN', 'DE', 'ES', 'FR', 'GR', 'IN', 'NO', 'US') ## GNI, PPP (current international $) ## GNI per capita, PPP (current international $) ## Population WDIdata <- WDI(indicator = c("NY.GNP.MKTP.PP.CD", "NY.GNP.PCAP.PP.CD", "SP.POP.TOTL"), start = 2000, end = 2021, country = countries) names(WDIdata) <- c("Country.Name", "iso2c", "Country.Code", "Year", "GNI.PPP", "GNI.capita", "Population") ## Bulk data downloaded from CW ## MtCO2e of CO2, Total including LUCF CWdata <- read.csv2("data/TimeSeries/CO2_ClimateWatch.csv") ## Filter by countries countries3 <- unique(WDIdata$Country.Code) CWdata <- subset(CWdata, Country %in% countries3) library("reshape2") CWdata <- melt(CWdata, id.vars = "Country", variable.name = "Year", value.name = "CO2") ## Remove the X in year values CWdata$Year <- as.numeric(substring(CWdata$Year, 2, 5)) CO2data <- merge(WDIdata, CWdata, by.x = c("Country.Code", "Year"), by.y = c("Country", "Year")) CO2data$CO2.capita <- with(CO2data, CO2 / Population * 1e6) CO2data$Country.Code <- factor(CO2data$Country.Code) CO2data$Country.Name <- factor(CO2data$Country.Name) save(CO2data, file = 'data/TimeSeries/CO2.RData')