This vignette shows how google charts can be used to display some spatio-temporal data sets used in the JSS paper on spacetime.
Google charts are interactive graphs in web pages that use google proprietary code. R pachage googleVis converts R data.frame objects into google charts. This vignette uses knitr and markdown to create a web page from the R-markdown source file. It was inspired by, and copies small sections of, the corresponding googleVis vignette
Set the googleVis options first to change the behaviour of plot.gvis,
so that only the chart component of the HTML file is written into
the output file.
library(googleVis)
op <- options(gvis.plot.tag='chart')
Following plot statements for gvis objects will automatically return
the HTML required for the 'knitted' output.
Geo charts work with country or region (state, administrative regions, e.g. NUTS1, formally ISO 3166-2) data. We will try to make the DE_NUTS1 data in spacetime ready for this. We read the table with names and codes, using help from stackoverflow:
library(XML)
url = "http://en.wikipedia.org/wiki/ISO_3166-2:DE"
ISO_3166_2_DE <- readHTMLTable(url, stringsAsFactors = FALSE)[[1]]
library(spacetime)
data(air) # loads rural and DE_NUTS1
Tbl <- gvisTable(ISO_3166_2_DE, options=list(width=400))
plot(Tbl)
Luckily, the regions in DE_NUTS are in the same (alphabetical) order
as in the table downloaded from wikipedia, so we can simply copy them
without need for matching. We will correct the two bilangual entries:
DE_NUTS1$name = ISO_3166_2_DE[,2]
DE_NUTS1$name[2] = "Bayern" # Not: "Bayern (Bavaria)"
DE_NUTS1$name[9] = "Niedersachsen" # Not: "Niedersachsen (Lower Saxony)"
Plotting Shape_Area, a variable present for all regions in DE_NUTS1:
library(googleVis)
M = gvisGeoMap(DE_NUTS1@data, locationvar = "name", numvar = "Shape_Area",
options=list(region="DE"))
plot(M)
reveals that Brandenburg is not recognized! We will now try with the ISO 3166-2 codes:
DE_NUTS1$code = ISO_3166_2_DE[,1]
M = gvisGeoMap(DE_NUTS1@data, locationvar = "code", numvar = "Shape_Area",
options=list(region="DE"))
plot(M)
which reveals that we now have all regions displayed. If however, we would omit Berlin, as in
noBerlin = DE_NUTS1[ISO_3166_2_DE[,1] != "DE-BE",]
M = gvisGeoMap(noBerlin@data, locationvar = "code", numvar = "Shape_Area",
options=list(region="DE"))
plot(M)
we see that Brandenburg now includes Berlin. Google, did you miss a hole here?
We can compute yearly average PM10 concentration over each of the states:
DE_NUTS1.years = STF(DE_NUTS1, as.Date(c("2008-01-01", "2009-01-01")))
agg = aggregate(rural[,"2008::2009"], DE_NUTS1.years, mean, na.rm=TRUE)
d = agg[,1]@data # select time step one, take attr table of resulting SpatialPolygonsDataFrame object
d$code = ISO_3166_2_DE[,1] # add region codes
M = gvisGeoMap(na.omit(d), locationvar = "code", numvar = "PM10",
options=list(region="DE",height=350)) # drop NA values for Geo chart
Tbl <- gvisTable(d, options=list(height=380, width=200))
plot(gvisMerge(M, Tbl, horizontal=TRUE))
|
|
|
The white states received no value and correspond to NA values
in the R object d; they were omitted by the na.omit call.
The Irish wind data has 12 stations, shown here on a map, colour denoting mean wind speed, symbols size its standard deviation:
library(sp)
data("wind", package = "gstat")
wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]])))
wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]])))
wind.loc$yx = paste(wind.loc$y, wind.loc$x, sep=":")
m = round(apply(wind,2,mean), 3)
sd = round(apply(wind,2,sd), 3)
wind.loc$mean = m[match(wind.loc$Code, names(m))]
wind.loc$sd = sd[match(wind.loc$Code, names(sd))]
M <- gvisGeoChart(wind.loc, "yx", "mean", "sd", "Station",
options = list(region = "IE"))
plot(M)
Two time lines are shown that have interesting interaction possibilities. For each of them, a subset of the full data sets (Irish wind and German air quality) are given, as the full data sets takes rather long to process (Rmd to html) and makes the browser slow loading the graph and when interacting with the data (2014).
We select the last three years of the wind data set, to not overload the browser:
# select:
wind = wind[wind$year > 76,] # select three years, 1976-1978
time = ISOdate(wind$year+1900, wind$month, wind$day)
wind.stack = stack(wind[,4:15])
names(wind.stack) = c("wind_speed", "station")
wind.stack$time = rep(time, 12)
M = gvisAnnotatedTimeLine(wind.stack, "time", "wind_speed", "station",
options = list(width = "1000px", height = "500px",
zoomStartTime=time[length(time)-365], zoomEndTime=max(time)))
plot(M)
The LineChart also supports time, and allows for identifying separate lines, when clicking the line or the station. It does not allow zooming, and does not work well for many observations or many lines.
wind$time = time
M = gvisLineChart(wind[1:200,], "time", names(wind)[4:9])
plot(M)
Line charts deal well with gaps in time series:
wind$time = time
wind[4:5, 7:9] = NA
M = gvisLineChart(wind[1:10,], "time", names(wind)[4:9])
plot(M)
The PM10 data in the spacetime vignette contain missing values. These are connected by straight lines in the Annotated Time Line, rendering this visualisation useless.
In contrast, Annotation Charts deal with gaps (missing values) in time series, and show them as discontinuities in the lines:
library(spacetime)
data(air)
d = as(rural[1:10,"2001"], "data.frame") # daily, 2001
M = gvisAnnotationChart(d, "time", "PM10", "sp.ID",
options = list(width = "1000px", height = "500px"))
# zoomStartTime=d$time[1], zoomEndTime=d$time[1]+365))
plot(M)
The result is still a bit flaky; just try it for the full time series.
The following example shows a motion chart of the Produc data used
in the spacetime manual, but does not convert the data. Time is simply
year (integer).
data("Produc", package = "plm")
M <- gvisMotionChart(Produc, "state", "year")
plot(M)
(Please note that the Motion Chart is only displayed when hosted on a web server, or if placed in a directory which has been added to the trusted sources in the security settings of Macromedia. See the googleVis package vignette for more details. )
## Set options back to original options
options(op)