Spatial and spatio-temporal objects in google charts

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

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?

An air quality example Geo chart

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.

Irish wind station mean and standard deviations

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)

Time lines

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).

Iris Wind in Time Annotation

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)

Irish wind in a Line Chart

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)

Rural PM10 data from spacetime: Annotation Chart

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.

Motion Chart

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)