% Discussion, add: % no type (field, object? point pattern or geostat?) % more ST types available (Galton 2008) % Graphs: move up, add all types. % panel data: do plm analysis on STFDF object? % % time from point -> interval; can this be done generic? % figures: % - what is implicitly the time interval? % - how do the classes look like on a s x t plot? % - reduction / simplification if t is regular - how does that work? % as.ts? in zoo? % TODO: % over, aggregate, interpolate? % http://www.data.gov/raw/1424# % Fri Jan 20 15:34:43 CET 2012; review round 1 from JSS. % Sun May 13 22:33:14 CEST 2012; review round 2 from JSS. %\documentclass[nogin,a4paper]{article} \documentclass[article]{jss} \Month{November} \Year{2012} \Volume{51} \Issue{7} \Submitdate{2011-10-05} \Acceptdate{2012-09-26} %\usepackage[OT1]{fontenc} %\usepackage[colorlinks=true,urlcolor=blue]{hyperref} %% need no \usepackage{Sweave.sty} \usepackage[utf8]{inputenc} \usepackage{alltt} \title{\bf \pkg{spacetime}: Spatio-Temporal Data in {\tt R}} \Shorttitle{\pkg{spacetime}: Spatio-Temporal Data in \proglang{R}} \author{ \includegraphics[width=5cm]{ifgi-logo_int} \hspace{.5cm} \includegraphics[width=4cm]{logo52n} \\ {\bf Edzer Pebesma} %Institute for Geoinformatics\\ %University of M\"{u}nster, Germany } %\date{\small \today } \Address{ Edzer Pebesma\\ Institute for Geoinformatics\\ University of M\"{u}nster\\ Weseler Strasse 253\\ 48151 M\"{u}nster, Germany\\ E-mail: \email{edzer.pebesma@uni-muenster.de}\\ URL: \url{http://ifgi.uni-muenster.de/} \\ \\ 52{\textdegree}North Initiative for Geospatial Open Source Software GmbH\\ Martin-Luther-King-Weg 24\\ 48155 Muenster, Germany\\ URL: \url{http://www.52north.org/} } \Abstract{ This document describes classes and methods designed to deal with different types of spatio-temporal data in \proglang{R} implemented in the \proglang{R} package \pkg{spacetime}, and provides examples for analyzing them. It builds upon the classes and methods for spatial data from package \pkg{sp}, and for time series data from package \pkg{xts}. The goal is to cover a number of useful representations for spatio-temporal sensor data, and results from predicting (spatial and/or temporal interpolation or smoothing), aggregating, or subsetting them, and to represent trajectories. The goals of this paper are to explore how spatio-temporal data can be sensibly represented in classes, and to find out which analysis and visualisation methods are useful and feasible. We discuss the time series convention of representing time intervals by their starting time only. This vignette is the main reference for the \proglang{R} package \code{spacetime}; it has been published as \cite{jss816}, but is kept up-to-date with the software.} \Keywords{Time series analysis, spatial data, spatio-temporal statistics, Geographic Information Systems} \begin{document} % \VignetteIndexEntry{ spacetime: Spatio-Temporal Data in R } % jss options, see https://www.jstatsoft.org/style <>= options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) set.seed(1331) @ \maketitle %\tableofcontents \section{Introduction} Spatio-temporal data are abundant, and easily obtained. Examples are satellite images of parts of the earth, temperature readings for a number of nearby stations, election results for voting districts and a number of consecutive elections, trajectories for people or animals possibly with additional sensor readings, disease outbreaks or volcano eruptions. \cite{schabenberger} argue that analysis of spatio-temporal data often happens {\em conditionally}, meaning that either first the spatial aspect is analysed, after which the temporal aspects are analysed, or reversed, but not in a joint, integral modelling approach, where space and time are not separated. As a possible reason they mention the lack of good software, data classes and methods to handle, import, export, display and analyse such data. This \proglang{R} \citep{RR} package is a start to fill this gap. Spatio-temporal data are often relatively abundant in either space, or time, but not in both. Satellite imagery is typically very abundant in space, giving lots of detail in high spatial resolution for large areas, but relatively sparse in time. Analysis of repeated images over time may further be hindered by difference in light conditions, errors in georeferencing resulting in spatial mismatch, and changes in obscured areas due to changed cloud coverage. On the other side, data from fixed sensors give often very detailed signals over time, allowing for elaborate modelling, but relatively little detail in space because a very limited number of sensors is available. The cost of an in situ sensor network typically depends primarily on its spatial density; the choice of the temporal resolution with which the sensors register signals may have little effect on total cost. Although for example \cite{botts} describe a number of open standards that allow the interaction with sensor data (describing sensor characteristics, requesting observed values, planning sensors, and processing raw sensed data to predefined events), the available statistical or geographic information system (GIS) software for this is in an early stage, and scattered. This paper describes an attempt to combine available infrastructure in the \proglang{R} statistical environment with ideas from statistical literature \citep{cressie} and data base literature \citep{rhg} to a set of useful classes and methods for manipulating, plotting and analysing spatio-temporal data. A number of case studies from different application areas will illustrate its use. The paper is structured as follows. Section~\ref{sec:longwide} describes how spatio-temporal data are usually recorded in tables. Section~\ref{layouts} describes a number of useful spatio-temporal layouts. Section~\ref{classes} introduces classes and methods for data, based on these layouts. Section~\ref{plots} presents a number of useful graphs for spatio-temporal data, and implementations for these. Section~\ref{support} discusses the spatial and temporal footprint, or support, of data, and how time intervals are dealt with in practice. Section~\ref{cases} presents a number of worked examples, some of which include statistical analysis on the spatio-temporal data. Section~\ref{further} points to further material, including further vignettes in package \pkg{spacetime} on spatio-temporal overlay and aggregation, and on using proxy data sets to PostgreSQL tables when data are too large for R. Section~\ref{discussion} finishes with a discussion. This paper is also found as a \href{https://cran.r-project.org/package=spacetime}{vignette} in package \pkg{spacetime}, which implements the classes and methods for spatio-temporal data described here. The vignette is kept up-to-date with the software. \section{How spatio-temporal data are recorded in tables} \label{sec:longwide} For reasons of simplicity, spatio-temporal data often come in the form of single tables. If this is the case, they come in one of three forms: \begin{description} \item[time-wide] where different columns reflect different moments in time, \item[space-wide] where different columns reflect different measurement locations or areas, or \item[long formats] where each record reflects a single time and space combination. \end{description} Alternatively, they may be stored in different, related tables, which is more typical for relational data bases, or in tree structures which is typical for xml files. We will now illustrate the different single-table formats with simple examples. \subsection{Time-wide format} Spatio-temporal data for which each location has data for each time can be provided in two so-called {\bf wide formats}. An example where a single column refers to a single moment or period in time is found in the North Carolina Sudden Infant Death Syndrome (sids) data set \citep{nc} available from package \pkg{sf}, which is in the {\bf time-wide format}: <<>>= if (require(foreign, quietly = TRUE) && require(sf, quietly = TRUE)) read.dbf(system.file("shape/nc.dbf", package="sf"))[1:5,c(5,9:14)] @ where {\bf columns} refer to a particular {\bf time}: \code{SID74} contains to the infant death syndrome cases for each county at a particular time period (1974-1984). \subsection{Space-wide format} The Irish wind data \citep{haslett} available from package \pkg{gstat} \citep{pebesma04}, for which the first six records and 9 of the stations (abbreviated by \code{RPT}, \code{VAL}, ...) are shown by <<>>= if (require(gstat, quietly = TRUE)) { data("wind", package = "gstat") wind[1:6,1:12] } @ are in {\bf space-wide format}: each {\em column} refers to another wind measurement {\bf location}, and the rows reflect a single time period; wind was reported as daily average wind speed in knots (1 knot = 0.5418 m/s). \subsection{Long format} Finally, panel data are shown in {\bf long form}, where the full spatio-temporal information is held in a single column, and other columns denote location and time. In the \code{Produc} data set \citep{baltagi}, a panel of 48 observations from 1970 to 1986 available from package \pkg{plm} \citep{croissant}, the first five records and nine columns are <<>>= if (require(plm, quietly = TRUE)) { data("Produc", package = "plm") Produc[1:5,1:9] } @ where the first two columns denote space and time (the default assumption for package \pkg{plm}), and e.g., \code{pcap} reflects private capital stock. None of these examples has strongly {\em referenced} spatial or temporal information: it is from the data alone not clear that the number \code{1970} refers to a year, or that \code{ALABAMA} refers to a state, and where this state is. Section \ref{cases} shows for each of these three cases how the data can be converted into classes with strongly referenced space and time information. \section{Space-time layouts} \label{layouts} In the following we will use the word spatial {\em feature} \citep{sfs} to denote a spatial entity. This can be a particular spatial point (location), a line or set of lines, a polygon or set of polygons, or a pixel (grid or raster cell). For a particular feature, one or more measurements are registered at particular moments in time. Four layouts of space-time data will be discussed next. Two of them reflect lattice layouts, one that is efficient when a particular spatial feature has data values for more than one time point, and one that is most efficient when all spatial feature have data values at each time point. Two others reflect irregular layouts, one of which specializes to trajectories (moving objects). \begin{figure} %[htb] \begin{center} <>= opar = par() par(mfrow=c(2,2)) # 1: s = 1:3 t = c(1, 1.75, 3, 4.5) g = data.frame(rep(t, each=3), rep(s,4)) col = 'blue' pch = 16 plot(g, xaxt = 'n', yaxt = 'n', xlab = "Time points", ylab = "Spatial features", xlim = c(.5,5.5), ylim = c(.5,3.5), pch = pch, col = col) abline(h=s, col = grey(.8)) abline(v=t, col = grey(.8)) points(g) axis(1, at = t, labels = c("1st", "2nd", "3rd", "4th")) axis(2, at = s, labels = c("1st", "2nd", "3rd")) text(g, labels = 1:12, pos=4) title("STF: full grid layout") # 2: s = 1:3 t = c(1, 2.2, 3, 4.5) g = data.frame(rep(t, each=3), rep(s,4)) sel = c(1,2,3,5,6,7,11) plot(g[sel,], xaxt = 'n', yaxt = 'n', xlab = "Time points", ylab = "Spatial features", xlim = c(.5,5.5), ylim = c(.5,3.5), pch = pch, col = col) abline(h=s, col = grey(.8)) abline(v=t, col = grey(.8)) points(g[sel,]) axis(1, at = t, labels = c("1st", "2nd", "3rd", "4th")) axis(2, at = s, labels = c("1st", "2nd", "3rd")) text(g[sel,], labels = paste(1:length(sel), "[",c(1,2,3,2,3,1,2),",",c(1,1,1,2, 2,3,4),"]", sep=""), pos=4) title("STS: sparse grid layout") # 3: s = c(1,2,3,1,4) t = c(1, 2.2, 2.5, 4, 4.5) g = data.frame(t,s) plot(g, xaxt = 'n', yaxt = 'n', xlab = "Time points", ylab = "Spatial features", xlim = c(.5,5.5), ylim = c(.5,4.5), pch = pch, col = col) #abline(h=s, col = grey(.8)) #abline(v=t, col = grey(.8)) arrows(t,s,0.5,s,.1,col='red') arrows(t,s,t,0.5,.1,col='red') points(g) axis(1, at = sort(unique(t)), labels = c("1st", "2nd", "3rd", "4th", "5th")) axis(2, at = sort(unique(s)), labels = c("1st,4th", "2nd", "3rd", "5th")) text(g, labels = 1:5, pos=4) title("STI: irregular layout") # 4: traj ns = 400 nt = 100 s = sort(runif(ns)) t = sort(runif(nt)) g = data.frame(t[1:30],s[1:30]) plot(g, xaxt = 'n', yaxt = 'n', xlab = "Time points", ylab = "Spatial features", type='l', col = 'blue', xlim = c(0,1), ylim = c(0,s[136])) lines(data.frame(t[41:60],s[31:50]), col = 'blue') lines(data.frame(t[91:100],s[51:60]), col = 'blue') lines(data.frame(t[21:40],s[61:80]), col = 'red') lines(data.frame(t[51:90],s[81:120]), col = 'red') lines(data.frame(t[11:25],s[121:135]), col = 'green') #abline(h=s, col = grey(.8)) #abline(v=t, col = grey(.8)) #arrows(t,s,0.5,s,.1,col='red') #arrows(t,s,t,0.5,.1,col='red') axis(1, at = sort(unique(t)), labels = rep("", length(t))) axis(2, at = sort(unique(s)), labels = rep("", length(s))) #text(g, labels = 1:5, pos=4) title("STT: trajectory") opar$cin = opar$cra = opar$csi = opar$cxy = opar$din = opar$page = NULL par(opar) @ \end{center} \caption{Four space-time layouts: (i) the top-left: full grid (\code{STF}) layout stores all space-time combinations; (ii) top-right: the sparse grid (\code{STS}) layout stores only the non-missing space-time combinations on a lattice; (iii) bottom-left: the irregular (\code{STI}) layout: each observation has its spatial feature and time stamp stored, in this example, spatial feature 1 is stored twice -- the fact that observations 1 and 4 have the same feature is not registered; (iv) bottom right: simple trajectories (\code{STT}), plotted against a common time axis. It should be noted that in both {\em gridded} layouts the grid is in space-time, meaning that spatial features {\em can} be gridded, but can also be any other non-gridded type (lines, points, polygons). } \label{fig:st} \end{figure} \subsection{Spatio-temporal full grids} A full space-time grid of observations for spatial features (points, lines, polygons, grid cells)\footnote{note that neither spatial features nor time points need to follow a regular layout} $s_i, i = 1,...,n$ and observation time $t_j, j = 1,...,m$ is obtained when the full set of $n \times m$ set of observations $z_k$ is stored, with $k=1,...,nm$. We choose to cycle spatial features first, so observation $k$ corresponds to feature $s_i$, $i=((k-1)\ \%\ n) + 1$ and with time moment $t_j$, $j=((k-1) / n)+ 1$, with $/$ integer division and \% integer division remainder (modulo). The $t_j$ are assumed to be in time order. In this data class (top left in Figure~\ref{fig:st}), for each spatial feature, the same temporal sequence of data is sampled. Alternatively one could say that for each moment in time, the same set of spatial entities is sampled. Unsampled combinations of (space, time) are stored in this class, but are assigned a missing value \code{NA}. It should be stressed that for this class (and the next) the word {\em grid} in {\em spatio-temporal grid} refers to the layout in space-time, not in space. Examples of phenomena that could well be represented by this class are regular (e.g., hourly) measurements of air quality at a spatially irregular set of points (measurement stations), or yearly disease counts for a set of administrative regions. An example where space is {\em gridded} as well could be a sequence of rainfall images (e.g., monthly sums), interpolated to a spatially regular grid. \subsection{Spatio-temporal sparse grids} A sparse grid has the same general layout, with measurements laid out on a space time grid (top right in Figure~\ref{fig:st}), but instead of storing the full grid, only non-missing valued observations $z_{k}$ are stored. For each $k$, an index $[i,j]$ is stored that refers which spatial feature $i$ and time point $j$ the value belongs to. Storing data this way may be efficient \begin{itemize} \item If full space-time lattices have many missing or trivial values (e.g., when one want to store features or grid cells where fires were registered, discarding those that did not), \item If a limited set of spatial features each have different time instances (e.g., to record the times of crime cases for a set of administrative regions), or, \item If for a limited set of times the set of spatial features varies (e.g., locations of crimes registered per year, or spatially misaligned remote sensing images). \end{itemize} \subsection{Spatio-temporal irregular data} Space-time irregular data cover the case where time and space points of measured values have no apparent organisation: for each measured value the spatial feature and time point is stored, as in the long format. This is equivalent to the (maximally) sparse grid where the index for observation $k$ is $[k,k]$, and hence can be dropped. For these objects, $n=m$ equals the number of records. Spatial features and time points need not be unique, but are replicated in case they are not. Any of the gridded types can be represented by this layout, in principle, but this would have the disadvantages that \begin{itemize} \item Spatial features and time points need to be stored for each data value, and would be redundant, \item The regular layout is lost, and needs be retrieved indirectly, \item Spatial and temporal selection would be inefficient, as the grid structure forms an index. \end{itemize} Examples of phenomena that are best served by this layout could be spatio-temporal point processes, such as crime or disease cases or forest fires. Other phenomena could be measurements from mobile sensors (in case the trajectory sequence is of no importance). \subsection{Interval time, moving objects, trajectories} In their book ``moving objects databases'', \cite{rhg} distinguish 10 different data types in space-time. In particular, they define for point features\footnote{ They obtain 10 types by adding the singular/atomic form of {\bf a} and {\bf b}, and doubling this set of 5 by distinguishing area from point features.}. \begin{itemize} \item[{\bf a}] Sets of events {\em without} temporal duration (time is an instant), e.g., accidents, lightning, birth, death; \item[{\bf b}] Sets of events {\em with} a temporal duration but no movement, e.g., a tree, a (point in the) capital of a country, people's home address; \item[{\bf c}] (Sets of) moving points, e.g., the trajectories of one or more persons, or birds. \end{itemize} To accomodate this typology we distinguish three cases, shown in figure \ref{fig:move}: \begin{itemize} \item[(i)] Time is instant and the feature is not moving (it may only exist at a time instant), \item[(ii)] Time is interval, objects do not move during this interval, \item[(iii)] Time is instant and features move (objects exist between time instants and may move) along a trajectory. \end{itemize} When time reflects intervals, it means that the spatial feature (spatial location or extent of the observation) or its associated data values does not change during this interval, but reflects the value or state {\em during} this interval. An examples is the yearly mean temperature of a country or of a set of locations, or the existence (duration) of a nation with a particular layout of its boundaries. \begin{figure}[htb] \begin{center} \includegraphics[width=.6\columnwidth]{move} \end{center} \caption{Time instant (top left), object movement (top right), time interval with consecutive (bottom left) or arbitrary (bottom right) intervals. $s_1$ refers to the first feature/location, $t_1$ to the first time instance or interval, thick lines indicate time intervals, arrows indicate movement. Filled circles denote start time, empty circles end times, intervals are right-closed. } \label{fig:move} \end{figure} Time {\em instants} can reflect the moments of change (e.g., the {\em start} of the meteorological summer), or events with a zero or negligible duration (e.g., an earthquake, a lightning). Movement reflects the fact that moving objects exist and may change location during a time interval. For moving object data, time instants reflect the location at a particular moment, and movement occurs between registered (time, feature) pairs, and must be continuous. Trajectories cover the case where sets of (irregular) space-time points form sequences, and depict a trajectory. Their grouping may be simple (e.g., the trajectories of two persons on different days), nested (for several objects, a set of trajectories representing different trips) or more complex (e.g., with objects that split, merge, or disappear). Examples of trajectories can be human trajectories, mobile sensor measurements (where the sequence is kept, e.g., to derive the speed and direction of the sensor), or trajectories of tornados where the tornado extent of each time stamp can be reflected by a different polygon. \section{Classes and methods for spatio-temporal data} \label{classes} The different layouts, or types, of spatio-temporal data discussed in Section~\ref{layouts} have been implemented in the \pkg{spacetime} \proglang{R} package, along with methods for import, export, coercion, selection, and visualisation. \subsection{Classes} The classes for the different layouts are shown in Figure~\ref{cls}. Similar to the classes of package \pkg{sp} \citep{pebesma,bivand}, the classes all derive from a base class \code{ST} which is not meant to represent actual data. The first order derived classes specify particular spatio-temporal geometries (i.e., only the spatial and temporal information), the second order derived classes augment each of these with actual data, in the form of a \code{data.frame}. \begin{figure}[htb] \begin{center} \includegraphics[width=.5\columnwidth]{cls} \end{center} \caption{Classes for spatio-temporal data in package \pkg{spacetime}. Arrows denote inheritance, lower side of boxes list slot name and type, green lines indicate supported coercions (both ways). } \label{cls} \end{figure} % is.regular. or not. % is.gridded. is sth different: 2D % what is wide and what is long format. To store temporal information, we chose to use objects of class \code{xts} in package \pkg{xts} \citep{ryan} for time, because \begin{itemize} \item It extends the functionality of package \pkg{zoo} \citep{zeileis}, \item It supports several basic types to represent time or date: \code{Date}, \code{POSIXt}, \code{timeDate}, \code{yearmon}, and \code{yearqtr}, \item It has good tools for {\em aggregation} over time using arbitrary aggregation functions, essentially deriving this from package \code{zoo} \citep{zeileis}, \item It has a flexible syntax to select time periods that adheres to ISO 8601\footnote{\url{http://en.wikipedia.org/wiki/ISO_8601}}. \end{itemize} An overview of the different time classes in \proglang{R} is found in \cite{ripleyhornik}. Further advice on which classes to use is found in \cite{grothendieck}, or in the \href{https://cran.r-project.org/web/views/TimeSeries.html}{CRAN task view on time series analysis}. For spatial interpolation, we used the classes deriving from \code{Spatial} in package \pkg{sp} \citep{pebesma,bivand} because \begin{itemize} \item They are the dominant set of classes in \proglang{R} for dealing with spatial data, \item They are interfaced to key external libraries, and, \item They provide a single interface to dealing with points, lines, polygons and grids. \end{itemize} We do not use \code{xts} or \code{Spatial} objects to {\em store} spatio-temporal data values, but we use \code{data.frame} to store data values. For purely temporal information the \code{xts} objects can be used, and for purely spatial information the \code{sp} objects can be used. These will be recycled appropriately when coercing to a long format \code{data.frame}. The spatial features supported by package \pkg{sp} are two-dimensional for lines and polygons, but may be higher (three-) dimensional for spatial points, pixels and grids. \subsection{Methods} The main methods for spatio-temporal data implemented in packages \pkg{spacetime} are listed in table \ref{methods}. Their usage is illustrated in examples that follow. \begin{table} \begin{center} \begin{tabular}{|lp{10cm}|} \hline method & what it does \\ \hline \code{stConstruct} & Creates STFDF or STIDF objects from single or multiple tables \\ \code{[[}, \code{\$}, \code{\$<-} & Select or replace data values \\ \code{[} & Select spatial and/or temporal subsets, and/or data variables \\ \code{as} & coerce to other spatio-temporal objects, \code{xts}, \code{Spatial}, \code{matrix}, or \code{data.frame} \\ \code{stplot} & create spatio-temporal plots, see Section~\ref{plots} \\ \code{over} & overlay: retrieve index or data values of one object at the locations and times of another \\ \code{aggregate} & aggregate data values over particular spatial, temporal, or spatio-temporal domains \\ \hline \end{tabular} \end{center} \caption{Methods for spatio-temporal data in package \pkg{spacetime}.} \label{methods} \end{table} \subsection{Creation} Construction of spatio-temporal objects essentially needs specification of the spatial, the temporal, and the data values. The documentation of \code{stConstruct} contains examples of how this can be done from long, space-wide, and time-wide tables, or from shapefiles. A simple toy example for a full grid layout with three spatial points and four time instances is given below. First, the spatial object is created: <<>>= sp = cbind(x = c(0,0,1), y = c(0,1,1)) row.names(sp) = paste("point", 1:nrow(sp), sep="") library(sp) sp = SpatialPoints(sp) @ Then, the time points are defined as four time stamps, one hour apart, starting Aug 5 2010, 10:00 GMT. <<>>= time = as.POSIXct("2010-08-05", tz = "GMT")+3600*(10:13) @ Next, a data.frame with the data values is created: <<>>= m = c(10,20,30) # means for each of the 3 point locations values = rnorm(length(sp)*length(time), mean = rep(m, 4)) IDs = paste("ID",1:length(values), sep = "_") mydata = data.frame(values = signif(values, 3), ID=IDs) @ And finally, the \code{STFDF} object can be created by: <>= library(spacetime) @ <>= library(spacetime) stfdf = STFDF(sp, time, data = mydata) @ In this case, as no \code{endTime} is specified, it is assumed that time reflects time instances. Altnatively, one could specify \code{endTime}, as in <<>>= stfdf = STFDF(sp, time, mydata, time+60) @ where the time intervals (Section \ref{intervals}) are one minute. When given a long table, \code{stConstruct} creates an \code{STFDF} object if all space and time combinations occur only once, or else an object of class \code{STIDF}, which might be coerced into other representations. \subsection{Overlay and aggregation} Aggregation of data values to a coarser spatial or temporal form (e.g., to a coarser grid, aggregating points over administrative regions, aggregating daily data to monthly data, or aggregation along an irregular set of space-time points) can be done using the method \code{aggregate}. To obtain the required aggregation predicate, i.e., the grouping of observations in space-time, the method \code{over} is implemented for objects deriving from \code{ST}. Grouping can be done based on spatial, temporal, or spatio-temporal predicates. It takes care of the case whether time reflects time instances or time intervals (see section \ref{intervals}). These methods effectively provide a spatio-temporal equivalent to what is known in geographic information science as the {\em spatial overlay}. \subsection[Space and time selection]{Space and time selection with \code{[}} The idea behind the \code{[} method for classes in \code{sp} was that objects would behave as much as possible similar to {\em matrix} or \code{data.frame} objects. For a \code{data.frame}, the expression \code{a[i,j]} selects row(s) \code{i} and column(s) \code{j}. For objects deriving from \code{Spatial}, rows were taken as the spatial features (points, lines, polygons, pixels) and columns as the data variables\footnote{a convention that was partially broken for class \code{SpatialGridDataFrame}, where \code{a[i,j,k]} could select the $k$-th data variable of the spatial grid selection with spatial grid row(s) \code{i} and column(s) \code{j}, {\em unless} the length of \code{i} equals the number of grid cells.}. For the spatio-temporal data classes described here, \code{a[i,j,k]} selects spatial features \code{i}, temporal instances \code{j}, and data variable(s) \code{k}. Unless \code{drop=FALSE} is added to such a call, selecting a single time or single feature results in an object that is no longer spatio-temporal, but either {\em snapshot} of a particular moment, or {\em history} at a particular feature \citep{galton}. Similar to selection on spatial objects in \pkg{sp} and time series objects in \pkg{xts}, space and time indices can be defined by index or boolean vectors, but by specifying spatial areas and time periods. For instance, the selection <>= air_quality[2:3, 1:10, "PM10"] @ yields air quality data for the second and third spatial features, and the first 10 time instances. The expressions <>= air_quality[Germany, "2008::2009", "PM10"] @ with \code{Germany} a \code{Spatial} object (e.g., a \code{SpatialPolygons}) defining Germany, selects the \code{PM10} measurements for the years 2008-9, lying in \code{Germany}. For {\bf trajectory} objects of class \code{STT} or \code{STTDF}, selection is slightly different: it is assumed that trajectories are being as complete. An expression \code{obj[1:3]} will select the first three full trajectories, \code{obj[Germany, "2008::2009", "Temp"]} selects the temperature attribute for all trajectories that intersect with Germany and fall at least partly in 2008-9. \subsection[Coercion to long and wide tables]{Coercion to long and wide tables} Spatio-temporal data objects can be coerced to the corresponding purely spatial objects. Objects of class \code{STFDF} will be represented in time-wide form, where only the first (selected) data variable is retained: <<>>= xs1 = as(stfdf, "Spatial") class(xs1) xs1 @ as time values are difficult to retrieve from these column names, this object gets the proper time values as an attribute: <<>>= attr(xs1, "time") @ Objects of class \code{STSDF} or \code{STIDF} will be represented in long form, where time is added as additional column: <<>>= x = as(stfdf, "STIDF") xs2 = as(x, "Spatial") class(xs2) xs2[1:4,] @ \section{Graphs of spatio-temporal data} \label{plots} \subsection[stplot: panels, space-time plots, animation]{\code{stplot}: panels, space-time plots, animation} The \code{stplot} method can create a few specialized plot types for the classes in the \code{spacetime} package. They are: \begin{description} \item[multi-panel plots] In this form, for each time step (selected) a map is plotted in a separate panel, and the strip above the panel indicates what the panel is about. The panels share $x$- and $y$-axis, no space needs to be lost by separating white space, and a common legend is used. Three types are implemented for STFDF data: \begin{itemize} \item The $x$ and $y$ axis denote space, an example for gridded data is shown in Figure~\ref{fig:wind}, for polygon data in Figure~\ref{fig:produc}. The \code{stplot} is a wrapper around \code{spplot} in package \code{sp}, and inherits most of its options, \item The $x$ and $y$ axis denote time and value; one panel for each spatial feature, colors may indicate different variables (\code{mode="tp"}); see Figure~\ref{fig:tpts} (left), \item The $x$ and $y$ axis denote time and value; one panel for each variable, colors may denote different features (\code{mode="ts"}); see Figure~\ref{fig:tpts} (right). \end{itemize} For both cases with time is on the $y$-axis (Figure~\ref{fig:tpts}), values over time for different variables or features are connected with lines, as is usual with time series plots. This can be changed to symbols by specifying \code{type='p'}. \item[space-time plots] Space-time plots show data in a space-time cross-section, with e.g., space on the $x$-axis and time on the $y$-axis. (See also Figure~\ref{fig:st}.) %R> pdf("ts.pdf") %R> stplot(Produc.st[1:4,,c("pcap", "hwy", "water", "util")],mode="ts") %R> pdf("tp.pdf") %R> stplot(Produc.st[1:4,,c("pcap", "hwy", "water", "util")],mode="tp") Hovm\"{o}ller diagrams \citep{hov} are an example of these for full space-time lattices, i.e., objects of class \code{STFDF}. To obtain such a plot, the arguments \code{mode} and \code{scaleX} should be considered; some special care is needed when only the x- or y-axis needs to be plotted instead of the spatial index (1...n); details are found in the \code{stplot} documentation. An example of a Hovm\"{o}ller-style plot with station index along the x-axis and time along the y-axis is obtained by <>= scales=list(x=list(rot = 45)) stplot(wind.data, mode = "xt", scales = scales, xlab = NULL) @ and shown in Figure~\ref{fig:hov}. Note that the $y$-axis direction is opposite to that of regular Hovm\"{o}ller plots. \item[animated plots] Animation is another way of displaying change over time; a sequence of \code{spplot}s, one for each time step, is looped over when the parameter \code{animate} is set to a positive value (indicating the time in seconds to pause between subsequent plots). \begin{figure}[htb] \begin{center} \includegraphics[scale=.6]{wind} \end{center} \caption{ Space-time interpolations of wind (square root transformed, detrended) over Ireland using a separable product covariance model, for 10 time points regularly distributed over the month for which daily data was considered (April, 1961).} \label{fig:wind} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=.49\columnwidth]{tp} \includegraphics[width=.49\columnwidth]{ts} \end{center} \caption{ Time series for four variables and four features plotted with \code{stplot}, with \code{mode="tp"} (left) and \code{mode="ts"} (right); see also Section~\ref{panel}. } \label{fig:tpts} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[scale=.4]{hov} \end{center} \caption{Space-time (Hovm\"{o}ller) plot of wind station data.} \label{fig:hov} \end{figure} \item[Time series plots] Time series plots are a fairly common type of plot in \proglang{R}. Package \code{xts} has a plot method that allows univariate time series to be plotted. Many (if not most) plot routines in \proglang{R} support time to be along the x- or y-axis. The plot in Figure~\ref{fig:windts} was generated by using package \pkg{lattice} \citep{sarkar}, and uses a colour palette from package \pkg{RColorBrewer} \citep{neuwirth}. <>= # code to create figure 5. library(lattice) if (require(RColorBrewer, quietly = TRUE)) { b = brewer.pal(12, "Set3") par.settings = list(superpose.symbol = list(col = b, fill = b), superpose.line = list(col = b), fontsize = list(text=9)) stplot(wind.data, mode = "ts", auto.key=list(space="right"), xlab = "1961", ylab = expression(sqrt(speed)), par.settings = par.settings) } @ \begin{figure}[htb] \begin{center} \includegraphics[scale=.8]{windts} \end{center} \caption{Time series plot of daily wind speed at 12 stations, used for interpolation in Figure~\ref{fig:wind}.} \label{fig:windts} \end{figure} \end{description} \section{Spatial footprint or support, time intervals, moving objects} \label{support} \subsection{Time periods or time instances} \label{intervals} Most data structures for time series data in \proglang{R} have, explicitly or implicitly, for each record a time stamp, not a time interval. The implicit assumption seems to be (i) the time stamp is a moment, (ii) this indicates either the real moment of measurement / registration, or the start of the interval over which something is aggregated (summed, averaged, maximized). For financial ``Open, high, low, close'' data, the ``Open'' and ``Close'' refer to the values at the {\em moment} the stock exchange opens and closes, meaning time instances, whereas ``high'' and ``low'' are {\em aggregated} values -- the minimum and maximum price over the time interval between opening and closing times. Package \code{lubridate} \citep{grolemund} allows one to explicitly define and compute with time intervals (e.g., \cite{allen}). It does not provide structures to attach these intervals to time series data. As \pkg{xts} does not support these times as index, \pkg{spacetime} does also not support it. According to \href{http://en.wikipedia.org/wiki/ISO_8601}{ISO 8601:2004}, a time stamp like \code{2010-05} refers to {\em the full} month of May, 2010, and so reflects a time {\em interval}. As a selection criterion, \code{xts} will include everything inside the following interval: <<>>= library(xts) .parseISO8601('2010-05') @ and this syntax lets one define, unambiguously, yearly, monthly, daily, hourly or minute intervals, but not e.g.\~10- or 30-minute intervals. For a particular interval, the full specification is needed: <<>>= .parseISO8601('2010-05-01T13:30/2010-05-01T13:39') @ When matching two sequences of time (Figure~\ref{fig:ti}) in order to overlay or aggregate, it matters whether each of the sequences reflect instances, one of them reflects time intervals and the other instances, or both reflect time intervals. All of these cases are accommodated for. Objects in \pkg{spacetime} register both (start) time and end time. By default, objects with gridded space-time layout (Figure~\ref{fig:st}) of class or deriving from \code{STF} or \code{STS} assume interval time, and \code{STI} and \code{STT} objects assume instance time. When no end times are supplied by creation and time intervals are assumed, the assumption is that time intervals are consecutive (Figure~\ref{fig:move}), and the last interval (for which no end time is present) has a length identical to the second last interval (Figures~\ref{fig:move} and \ref{fig:ti}). \begin{figure}[htb] \begin{center} \includegraphics[scale=1.1]{ti} \end{center} \caption{Matching two time sequences, assuming \code{x} reflects time intervals, and \code{y} reflects time instances. Note that the last interval extends the last time instance of \code{x}. } \label{fig:ti} \end{figure} \subsection{Spatial support} All examples above work with spatial points, i.e., data having a point support. The assumption of data having points support is implicit for \code{SpatialPoints} features. For polygons, the assumption will be that values reflect aggregates (e.g., sums, or averages) over the polygon. For gridded data, it is ambiguous whether the value at the grid cell centre is meant (e.g., for DEM data) or an aggregate over the grid cell (typical for remote sensing imagery). The \code{Spatial*} objects of package \pkg{sp} have no {\em explicit} information about the spatial support. % What if time intervals are irregular, and values denoting aggregate? % What is the time interval for the last measurement? % Time intervals, and their relationship; possibility of overlapping % (duplicate) information. % Anyway, how do we find/identify/deal with possibility of duplicates? % Resampling? % Spatial/temporal/spatio-temporal aggregation? \section{Worked examples} \label{cases} This section shows how existing data in various formats can be converted into ST classes, and how they can be analysed and/or visualised. \subsection{North Carolina SIDS} \label{sec:nc} As an example, the North Carolina Sudden Infant Death Syndrome (sids) data will be used. These data were first analysed by \cite{sids}, and first published and analysed in a spatial setting by \cite{cressiechan}. The data are sparse in time (aggregated to 2 periods of unequal length, according to the documentation in package \code{spdep}), but have polygons in space. First, we will prepare the spatial data: <<>>= if (require(sf, quietly = TRUE)) { fname = system.file("gpkg/nc.gpkg", package = "sf")[1] nc = as(sf::st_read(fname), "Spatial") } @ then, we construct the time sequence: <<>>= time = as.POSIXct(c("1974-07-01", "1979-07-01"), tz = "GMT") endTime = as.POSIXct(c("1978-06-30", "1984-06-30"), tz = "GMT") @ and we construct the data values table, in long form, by <<>>= data = data.frame( BIR = c(nc$BIR74, nc$BIR79), NWBIR = c(nc$NWBIR74, nc$NWBIR79), SID = c(nc$SID74, nc$SID79)) @ These three components are put together by function \code{STFDF}: <<>>= nct = STFDF(sp = as(nc, "SpatialPolygons"), time, data, endTime) @ \subsection{Panel data} \label{panel} The panel data discussed in Section~\ref{sec:longwide} are imported as a full spatio-temporal \code{data.frame} (\code{STFDF}), and linked to the proper state polygons of maps. We can obtain the states polygons from package \pkg{map} \citep{maps} by: <<>>= if (require(maps, quietly = TRUE)) { states.m <- map('state', plot=FALSE, fill=TRUE) IDs <- sapply(strsplit(states.m$names, ":"), function(x) x[1]) } if (require(sf, quietly = TRUE)) states <- as(st_geometry(st_as_sf(states.m, IDs=IDs)), "Spatial") @ we obtain the time points by: <<>>= yrs = 1970:1986 time = as.POSIXct(paste(yrs, "-01-01", sep=""), tz = "GMT") @ We obtain the data table (already in long format) by <<>>= if (require(plm, quietly = TRUE)) { data("Produc") } @ When combining all this information, we do not need to reorder states because \code{states} and \code{Produc} order states alphabetically. We need to de-select District of Columbia, which is not present in \code{Produc} table (record 8): <<>>= if (require(plm, quietly = TRUE)) # deselect District of Columbia, polygon 8, which is not present in Produc: Produc.st <- STFDF(states[-8], time, Produc[order(Produc[,2], Produc[,1]),]) @ <<>>= if (require(plm, quietly = TRUE) && require(RColorBrewer, quietly = TRUE)) stplot(Produc.st[,,"unemp"], yrs, col.regions = brewer.pal(9, "YlOrRd"),cuts=9) @ produces the plot shown in Figure~\ref{fig:produc}. \begin{figure}[htb] \begin{center} \includegraphics[scale=.8]{produc} \end{center} \caption{Unemployment rate per state, over the years 1970-1986.} \label{fig:produc} \end{figure} Time and state were not removed from the data table on construction; printing these data after coercion to \code{data.frame} can then be used to verify that time and state were matched correctly. The routines in package \pkg{plm} can be used on the data, when back transformed to a \code{data.frame}, when \code{index} is used to specify which variables represent space and time (the first two columns from the \code{data.frame} no longer contain state and year). For instance, to fit a panel linear model for gross state products (gsp) to private capital stock (pcap), public capital (pc), labor input (emp) and unemployment rate (unemp), we get <<>>= if (require(plm, quietly = TRUE)) zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = as.data.frame(Produc.st), index = c("state", "year")) @ where the output of \code{summary(zz)} is left out for brevity. More details are found in \cite{croissant} and \cite{millo}. \subsection{Interpolating Irish wind} \label{sec:wind} This worked example is a modified version of the analysis presented in \code{demo(wind)} of package \pkg{gstat} \citep{pebesma04}. This demo is rather lengthy and reproduces much of the original analysis in \cite{haslett}. Here, we will reduce the material and focus on the use of spatio-temporal classes. First, we will load the wind data from package \code{gstat}. It has two tables, station locations in a \code{data.frame}, called \code{wind.loc}, and daily mean wind speed in \code{data.frame} \code{wind}. We now convert character representation (such as \code{51d56'N}) to proper numerical coordinates, and convert the station locations to a \code{SpatialPointsDataFrame} object. A plot of these data is shown in Figure~\ref{fig:windloc}. <<>>= if (require(gstat, quietly = TRUE)) { data("wind") wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]]))) wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]]))) coordinates(wind.loc) = ~x+y proj4string(wind.loc) = "+proj=longlat +datum=WGS84" } @ \begin{figure}[htb] \begin{center} <>= if (require(mapdata, quietly = TRUE)) { plot(wind.loc, xlim = c(-11,-5.4), ylim = c(51,55.5), axes=T, col="red", cex.axis =.7) map("worldHires", add=TRUE, col = grey(.5)) text(coordinates(wind.loc), pos=1, label=wind.loc$Station, cex=.7) } else plot(1) @ \end{center} \caption{Station locations for Irish wind data.} \label{fig:windloc} \end{figure} The first thing to do with the wind speed values is to reshape these data. Unlike the North Carolina SIDS data of section \ref{sec:nc}, for we have few spatial and many time points, and so the data in \code{data.frame} \code{wind} come in space-wide form with stations time series in columns: <<>>= if (require(gstat, quietly = TRUE)) { wind[1:3,] } @ We will recode the time columns to an appropriate time data structure, <<>>= if (require(gstat, quietly = TRUE)) { wind$time = ISOdate(wind$year+1900, wind$month, wind$day) wind$jday = as.numeric(format(wind$time, '%j')) } @ and then subtract a smooth time trend of daily means (not exactly equal, but similar to the trend removal in the original paper): <<>>= if (require(gstat, quietly = TRUE)) { stations = 4:15 windsqrt = sqrt(0.5148 * as.matrix(wind[stations])) # knots -> m/s Jday = 1:366 windsqrt = windsqrt - mean(windsqrt) daymeans = sapply(split(windsqrt, wind$jday), mean) meanwind = lowess(daymeans ~ Jday, f = 0.1)$y[wind$jday] velocities = apply(windsqrt, 2, function(x) { x - meanwind }) } @ Next, we will match the wind data to its location, by connecting station names to location coordinates, and create a spatial points object: <<>>= if (require(gstat, quietly = TRUE)) { wind.loc = wind.loc[match(names(wind[4:15]), wind.loc$Code),] pts = coordinates(wind.loc[match(names(wind[4:15]), wind.loc$Code),]) rownames(pts) = wind.loc$Station pts = SpatialPoints(pts, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84")) } @ Then, we project the longitude/latitude coordinates and country boundary to UTM zone 29, using \code{st_transform} in package \pkg{sf} for coordinate transformation: <<>>= utm29 = "+proj=utm +zone=29 +datum=WGS84 +ellps=WGS84" pts.sfc = st_transform(st_as_sfc(pts), utm29) pts = as(pts.sfc, "Spatial") # back to sp @ And now we can construct the spatio-temporal object from the space-wide table with velocities: <<>>= if (require(gstat, quietly = TRUE)) { wind.data = stConstruct(velocities, space = list(values = 1:ncol(velocities)), time = wind$time, SpatialObj = pts, interval = TRUE) class(wind.data) } @ For plotting purposes, we can obtain country boundaries from package \pkg{maps}: <<>>= if (require(sf, quietly = TRUE) && require(mapdata, quietly = TRUE)) { m.sf = st_as_sf(map("worldHires", xlim = c(-11.5,-6.0), ylim = c(51.3,55.0), plot=FALSE), fill = FALSE) m.sf = st_transform(m.sf, utm29) m = as(m.sf, "Spatial") } @ For interpolation, we can define a grid over the area: <<>>= if (require(gstat, quietly = TRUE)) grd = SpatialPixels(SpatialPoints(makegrid(m, n = 300)), proj4string = proj4string(m)) @ Next, we (arbitrarily) restrict observations to those of April 1961: <<>>= if (require(gstat, quietly = TRUE)) wind.data = wind.data[, "1961-04"] @ and choose 10 time points from that period to form the spatio-temporal prediction grid: <<>>= if (require(gstat, quietly = TRUE)) { n = 10 library(xts) tgrd = seq(min(index(wind.data)), max(index(wind.data)), length=n) pred.grd = STF(grd, tgrd) } @ We will interpolate with a separable exponential covariance model, with ranges 750 km and 1.5 days: <<>>= if (require(gstat, quietly = TRUE)) { v = vgmST("separable", space = vgm(1, "Exp", 750000), time = vgm(1, "Exp", 1.5 * 3600 * 24), sill=0.6) wind.ST = krigeST(values ~ 1, wind.data, pred.grd, v) colnames(wind.ST@data) <- "sqrt_speed" } @ % wind.ST = STFDF(grd, tgrd, data.frame(sqrt_speed = pred)) then creates the \code{STFDF} object with interpolated values, the results of which are shown in Figure~\ref{fig:wind}, created by <>= if (require(gstat, quietly = TRUE)) { layout = list(list("sp.lines", m, col='grey'), list("sp.points", pts, first=F, cex=.5)) stplot(wind.ST, col.regions=brewer.pal(11, "RdBu")[-c(10,11)], at=seq(-1.375,1,by=.25), par.strip.text = list(cex=.7), sp.layout = layout) } @ <>= pdf("wind.pdf", height=4.5) if (require(gstat, quietly = TRUE)) { layout = list(list("sp.lines", m, col='grey'), list("sp.points", pts, first=F, cex=.5)) print(stplot(wind.ST, col.regions=brewer.pal(11, "RdBu")[-c(10,11)], at=seq(-1.375,1,by=.25), par.strip.text = list(cex=.7), sp.layout = layout)) } else plot(1) dev.off() @ <>= pdf("windts.pdf", height = 4) if (require(gstat, quietly = TRUE)) { library(lattice) library(RColorBrewer) b = brewer.pal(12,"Set3") par.settings = list(superpose.symbol = list(col = b, fill = b), superpose.line = list(col = b), fontsize = list(text=9)) print(stplot(wind.data, mode = "ts", auto.key=list(space="right"), xlab = "1961", ylab = expression(sqrt(speed)), par.settings = par.settings)) } else plot(1) dev.off() @ <>= if (require(gstat, quietly = TRUE)) { pdf("hov.pdf") scales=list(x=list(rot=45)) stplot(wind.data, mode = "xt", scales = scales, xlab = NULL, col.regions=brewer.pal(11, "RdBu"),at = seq(-1.625,1.125,by=.25)) dev.off() } @ \subsection{Calculation of EOFs} Empirical orthogonal functions from \code{STFDF} objects can be computed in spatial form (default), e.g., from data values <>= if (require(gstat, quietly = TRUE)) eof.data = eof(wind.data) @ or alternatively from modelled, or interpolated values: <>= if (require(gstat, quietly = TRUE)) eof.int = eof(wind.ST) @ By default, spatial EOFs are competed; alternatively they can be obtained in temporal form: <>= if (require(gstat, quietly = TRUE)) eof.xts = eof(wind.ST, "temporal") @ the resulting object is of the appropriate subclass of \code{Spatial} in the spatial form, or of class \code{xts} in the temporal form. Figure~\ref{fig:eof} shows the 10 spatial EOFs obtained from the interpolated wind data of Figure~\ref{fig:wind}. \begin{figure}[htb] \begin{center} <>= if (require(gstat, quietly = TRUE)) { print(spplot(eof.int[1:4], col.regions=bpy.colors(), par.strip.text = list(cex=.5), as.table = TRUE, sp.layout = layout)) } else { plot(1) } @ \end{center} \caption{EOFs of space-time interpolations of wind over Ireland (for spatial reference, see Figure~\ref{fig:wind}), for the 10 time points at which daily data was chosen above (April, 1961).} \label{fig:eof} \end{figure} \subsection[Trajectory data: ltraj in adehabitatLT]{Trajectory data: \code{ltraj} in package \pkg{adehabitatLT}} Trajectory objects of class \code{ltraj} in package \pkg{adehabitatLT} \citep{calenge} are lists of bursts, sets of sequential, connected space-time points at which an object is registered. An example \code{ltraj} data set is obtained by\footnote{taken from \pkg{adehabitatLT}, \code{demo(mangltraj)}}: <<>>= if (require(adehabitatLT, quietly = TRUE)) { data("puechabonsp") locs = puechabonsp$relocs xy = coordinates(locs) da = as.character(locs$Date) da = as.POSIXct(strptime(as.character(locs$Date),"%y%m%d", tz = "GMT")) ltr = as.ltraj(xy, da, id = locs$Name) foo = function(dt) dt > 100*3600*24 l2 = cutltraj(ltr, "foo(dt)", nextr = TRUE) l2 } @ and these data, converted to \code{STTDF} can be plotted, as panels by \code{time} and \code{id} by <>= sttdf = as(l2, "STTDF") stplot(sttdf, by="time*id") @ which is shown in Figure~\ref{fig:stptr}. \begin{figure}[htb] \begin{center} <>= if (require(adehabitatLT, quietly = TRUE)) { sttdf = as(l2, "STTDF") print(stplot(sttdf, by="time*id")) } else plot(1) @ \end{center} \caption{Trajectories, split by id (rows) and by time (columns).} \label{fig:stptr} \end{figure} \subsection{Country shapes in cshapes} The \pkg{cshapes} \citep{weidmann} package contains a GIS dataset of country boundaries (1946-2008), and includes functions for data extraction and the computation of distance matrices. The data set consist of a \code{SpatialPolygonsDataFrame}, with the following data variables: <<>>= if (require(cshapes, quietly = TRUE)) { library(sf) cs = cshp() print(names(cs)) } @ where two data bases are used, ``COW'' (correlates of war project\footnote{Correlates of War Project. 2008. State System Membership List, v2008.1. Online, \url{http://correlatesofwar.org/}}) and ``GW'' \cite{gleditsch}. The variables COWSMONTH and COWEMONTH denote the start month and end month, respectively, according to the COW data base. In the following fragment, we create the spatio-temporal object using begin- and end-times: <<>>= if (require(cshapes, quietly = TRUE)) st = STIDF(geometry(as(cs, "Spatial")), as.POSIXct(cs$start), as.data.frame(cs), as.POSIXct(cs$end)) @ A possible query would be which countries are found at 7\textdegree East and 52\textdegree North, <<>>= if (require(cshapes, quietly = TRUE)) { pt = SpatialPoints(cbind(7, 52), CRS(proj4string(st))) as.data.frame(st[pt,,1:5]) } @ \section{Further material} \label{further} Searching past email discussion threads on the \href{https://stat.ethz.ch/mailman/listinfo/r-sig-geo}{\code{r-sig-geo}} (R Special Interest Group on using GEOgraphical data and Mapping) email list may be a good way to look for further material, before one considers posting questions. Search strings, e.g., on the google search engine may look like: \code{spacetime site:stat.ethz.ch} where the search keywords should be made more precise. The excellent book {\em Statistics for spatio-temporal data} \citep{cressie} provides a large number of methods for the analysis of mainly geostatistical data. A demo script, which can be run by <>= library(spacetime) demo(CressieWikle) @ downloads the data from the book web site, and reproduces a number of graphs shown in the book. It should be noted that the the book examples only deal with \code{STFDF} objects. Section~\ref{sec:wind} contains an example of a spatial interpolation with a spatio-temporal separable or product-sum covariance model. The functions for this are found in package \pkg{gstat}, and more information is found through <>= if (require(gstat, quietly = TRUE)) { vignette("st") } @ An example where (potentially large) data sets are proxied through R objects is given in a vignette in the \pkg{spacetime} package, obtained by <>= library(spacetime) vignette("stpg") @ A proxy object is an object that contains no data, but only references to tables in a data base. Selections on this object are translated into SQL statements that return the actually selected data. This way, the complete data set does not have to be loaded in memory (R), but can be processed part by part. Selection in the data base uses indexes on the spatial and temporal references. \label{overlay} Examples of overlay and aggregation methods for spatio-temporal data are further detailed in a separate vignette, obtained by <>= library(spacetime) vignette("sto") @ It illustrates the methods with daily air quality data taken from the European air quality data base, for 1998-2009. Aggregations are temporal, spatial, or both. %<<>>= %http://www.nws.noaa.gov/geodata/catalog/national/data/s_01au07.zip %@ % \section{TODO} % write aggregate for all 3? aggregate over time -> defer to xts? % split method for STIDF % write tests for points, lines, polygons, pixels, grid % How to do spatial selection for grids - use index, or c(rows,cols) % as spatial selector? % xts: aggregate, shift to nearest grid point, % s/t: overlay? (where,when)->index of values available? xts accepts % time index, sp will need overlay. % generic: plot, % DONE: % summary, show (automatic), spTransform (NOT: requires rgdal dependency), % stplot (as function) % Plotting methods: stplot -> interface to spplot with panel for % each time step (single attribute), or double with [attribute, time] % as panel entries. % plot -> plot.xts/plot.zoo interface? % % plot time series with little maps (micro maps) indicating region/point? ggplot? % plot maps with little time series plotted at particular points? % \subsection{coercion to external classes} % spatial panel model. \section{Discussion} \label{discussion} Handling and analyzing spatio-temporal data is often complicated by the size and complexity of these data. Also, data may come in many different forms, they may be time-rich, space-rich, and come as sets of space-time points or as trajectories. Building on existing infrastructure for spatial and temporal data, we have successfully implemented a coherent set of classes for spatio-temporal data that covers regular space-time layouts, partially regular (sparse) space-time layouts, irregular space-time layouts and trajectory data. The set is flexible in the sense that several representations of space (points, lines, polygons, grid) and time (POSIXt, Date, timeDate, yearmon, yearqtr) can be combined. We have given examples for constructing objects of these classes from various data sources, coercing them from one to another, exporting them to spatial or temporal representations, as well as visualising them in various forms. We have also shown how one can go from one form into another by ways of prediction based on a statistical model, using an example on spatio-temporal geostatistical interpolation. In addition to spatio-temporally varying information, objects of the classes can contain data values that are purely spatial or purely temporal. Selection can be done based on spatial features, time (intervals), or data variables, and follows a logic similar to that for selection on data tables (\code{data.frame}s). Challenges that remain include \begin{itemize} \item The representation of spatio-temporal polygons in a consistent way, i.e., such that each point in space-time refers to one and only one space-time feature, \item Dealing with complex developments, such as merging, splitting, and death and birth of objects (further examples are found in \cite{galton}), \item Explicitly registering the support, or footprint of spatio-temporal data, \item Annotating objects such that incorrect operations (such as the interpolation of a point process, or the weighted density estimates on a geostatistical process) can lead to warning or error messages, \item Making handling of massive data sets easier, and implementing efficient spatio-temporal indexes for them, \item Integrating package \pkg{spacetime} with other packages dealing with specific spatio-temporal classes such as \pkg{raster} and \pkg{surveillance}. \end{itemize} The classes and methods presented in this paper are a first attempt to cover a number of useful cases for spatio-temporal data. In a set of case studies it is demonstrated how they can be used, and can be useful. As software development is often opportunistic, we admittedly picked a lot of low hanging fruits, and a number of large challenges remain. We hope that these first steps will help discovering and identifying these more complex use cases. \section*{Acknowledgements} %Michael Sumner provided helpful comments on the trip example. Members from the spatio-temporal modelling lab of the Institute for Geoinformatics of the University of M\"{u}nster (Ben Gr\"{a}ler, Katharina Henneb\"{o}hl, Daniel N\"{u}st, and S\"{o}ren Gebbert) contributed in several useful discussions. Participants to the workshop {\em Handling and analyzing spatio-temporal data in \proglang{R}}, held in M\"{u}nster on Mar 21-22, 2011, are gratefully acknowledged. Several anonymous reviewers provided valuable comments. \bibliography{jss816} \end{document} http://dna.fernuni-hagen.de/Lehre-offen/Kurse/1675/KE1.pdf : We now characterize application data with respect to their “shape” in this 3D space, obtaining the following categories: 1. Events in space and time – (point, instant). Examples are archeological discoveries, plane crashes, volcano eruptions, earthquakes (at a large scale where the duration is not relevant). (STI; STIDF) 2. Locations valid for a certain period of time – (point, period). Examples are: cities built at some time, still existing or destroyed; construction sites (e.g., of buildings, highways); branches, offices, plants, or stores of a company; coal mines, oil wells, being used for some time; or “immovables”, anything that is built at some place and later destroyed. (STL; STLDF) 3. Set of location events – sequence of (point, instant). Entities of class (1) when viewed collectively. For example, the volcano eruptions of the last year. (STI; STIDF) 4. Stepwise constant locations – sequence of (point, period). Examples are: the capital of a country; the headquarter of a company; the accomodations of a traveler during a trip; the trip of an email message (assuming transfer times between nodes are zero). (STL; STLDF? ) 5. Moving entities – moving point. Examples are people, planes, cars, etc., see Table 1.7. (STT; STTDF) 6. Region events in space and time – (region, instant). E.g., a forest fire at large scale. Figure 1.11: (a) Discretely changing point and region (b) Continuously changing point and region (STT; STTDF, or STI/STIDF) 7. Regions valid for some period of time – (region, period). For example, the area closed for a certain time after a traffic accident. (STL; STLDF) 8. Set of region events – sequence of (region, instant). For example, the Olympic games viewed collectively, at a large scale. (STI; STIDF) 9. Stepwise constant regions – sequence of (region, period). For example, countries, real estate (changes of shape only through legal acts), agricultural land use, etc. (STL? STI? STT?) 10. Moving entities with extent – moving region. For example, forests (growth); forest fires at small scale (i.e., we describe the development); people in history; see Table 1.8. (STT*)