% dimensions(x) returns number of spatial dimensions
% y = transform(x, "proj4string")
% bbox(x)
% coordinates(x) ; <-
% rings(x) ; <-
% method to retrieve lines? --> Lines()?
% gridded(x)  ; <-
% 
\documentclass{article}

\usepackage{graphicx}
\usepackage[colorlinks=true,urlcolor=blue]{hyperref}

% \VignetteIndexEntry{ sp: overlay and aggregation }

\usepackage{color}

\usepackage{Sweave}
\newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\newcommand{\code}[1]{{\tt #1}}
\let\pkg=\strong

\title{\bf Map overlay and \\ spatial aggregation in {\tt sp}}
\author{Edzer Pebesma\footnote{Institute for Geoinformatics,
University of Muenster, Weseler Strasse 253, 48151 M\"{u}nster, Germany.
{\tt edzer.pebesma@uni-muenster.de}}}
\date{\today}

\begin{document}

\maketitle

\begin{abstract}
Numerical ``map overlay'' combines spatial features from one map
layer with the attribute (numerical) properties of another. This
vignette explains the R method ``over'', which provides a consistent
way to retrieve indices or attributes from a given spatial object
(map layer) at the locations of another. Using this, the R generic
``aggregate'' is extended for spatial data, so that any spatial
properties can be used to define an aggregation predicate, and
any R function can be used as aggregation function.
\end{abstract}

\tableofcontents

\section{Introduction}
According to the free e-book by Davidson (2008),
\begin{quotation} 
{\em An overlay is a clear sheet of plastic or semi-transparent paper. It
is used to display supplemental map and tactical information related
to military operations. It is often used as a supplement to orders
given in the field. Information is plotted on the overlay at the
same scale as on the map, aerial photograph, or other graphic being
used. When the overlay is placed over the graphic, the details
plotted on the overlay are shown in their true position. }
\end{quotation}
This suggests that {\em map overlay} is concerned with combining
two, or possibly more, map layers by putting them on top of each
other. This kind of overlay can be obtained in R e.g. by plotting
one map layer, and plotting a second map layer on top of it. If the
second one contains polygons, transparent colours can be used to
avoid hiding of the first layer. When using the {\tt spplot} command,
the {\tt sp.layout} argument can be used to combine multiple layers.

O'Sullivan and Unwin (2003) argue in chapter 10 (Putting maps
together: map overlay) that map overlay has to do with the
combination of two (or more) maps. They mainly focus on the
combination of the selection criteria stemming from several map
layers, e.g. finding the deciduous forest area that is less than
5 km from the nearest road. They call this {\em boolean overlays}.

One could look at this problem as a polygon-polygon overlay, where we
are looking for the intersection of the polygons with the deciduous
forest with the polygons delineating the area less than 5 km from
a road. Other possibilities are to represent one or both coverages
as grid maps, and find the grid cells for which both criteria are
valid (grid-grid overlay). A third possibility would be that one
of the criteria is represented by a polygon, and the other by a
grid (polygon-grid overlay, or grid-polygon overlay). In the end, as
O'Sullivan and Unwin argue, we can overlay any spatial type (points,
lines, polygons, pixels/grids) with any other. In addition, we can
address spatial attributes (as the case of grid data), or only the
geometry (as in the case of the polygon-polygon intersection).

This vignette will explain how the {\tt over} method in package
{\tt sp} can be used to compute map overlays, meaning that instead
of overlaying maps visually, the digital information that comes
from combining two digital map layers is retrieved. From there,
methods to {\em aggregate} (compute summary statistics; Heuvelink
and Pebesma, 1999) over a spatial domain will be developed and
demonstrated.  Pebesma (2012) describes overlay and aggregation
for spatio-temporal data.

\section{Geometry overlays}
We will use the word {\em geometry} to denote the purely spatial
characteristics, meaning that attributes (qualities, properties of
something at a particular location) are ignored. With {\em location}
we denote a point, line, polygon or grid cell. Section \ref{attr}
will discuss how to retrieve and possibly aggregate or summarize
attributes found there.

Given two geometries, {\tt A} and {\tt B}, the following equivalent
commands
<<eval=FALSE>>=
A %over% B
over(A, B)
@
retrieve the geometry (location) indices of \code{B} at the locations
of \code{A}. More in particular, an integer vector of length
{\tt length(A)} is returned, with {\tt NA} values for locations in {\tt A}
not matching with locations in {\tt B} (e.g. those points outside
a set of polygons). 

Selecting points of \code{A} {\em inside} or {\em on} some geometry
\code{B} (e.g.  a set of polygons) {\tt B} is done by
<<eval=FALSE>>=
A[B,]
@
which is short for
<<eval=FALSE>>=
A[!is.na(over(A,B)),]
@
We will now illustrate this with toy data created by
<<keep.source=TRUE>>=
library(sp)
x = c(0.5, 0.5, 1.2, 1.5)
y = c(1.5, 0.5, 0.5, 0.5)
xy = cbind(x,y)
dimnames(xy)[[1]] = c("a", "b", "c", "d")
pts = SpatialPoints(xy)

xpol = c(0,1,1,0,0)
ypol = c(0,0,1,1,0)
pol = SpatialPolygons(list(
	Polygons(list(Polygon(cbind(xpol-1.05,ypol))), ID="x1"),
	Polygons(list(Polygon(cbind(xpol,ypol))), ID="x2"),
	Polygons(list(Polygon(cbind(xpol,ypol-1.05))), ID="x3"),
	Polygons(list(Polygon(cbind(xpol+1.05,ypol))), ID="x4"),
	Polygons(list(Polygon(cbind(xpol+.4, ypol+.1))), ID="x5")
	))
@
and shown in figure \ref{fig:toy}.

\begin{figure}[htb]
<<fig=TRUE,echo=FALSE>>=
plot(pol, xlim = c(-1.1, 2.1), ylim = c(-1.1, 1.6), border=2:6, axes=TRUE)
points(pts, col='red')
text(c(-1,0.1,0.1,1.1,0.45), c(0,0,-1.05,0,0.1), 
	c("x1", "x2", "x3", "x4", "x5"))
text(coordinates(pts), pos=1, row.names(pts))
@
\caption{ Toy data: points (a-d), and (overlapping) polygons (x1-x5) }
\label{fig:toy}
\end{figure}

Now, the polygons \code{pol} in which points \code{pts} lie are
<<>>=
over(pts, pol)
@
As points \code{b} and \code{c} fall in two overlapping polygons, we can
retrieve the complete information as a list:
<<>>=
over(pts, pol, returnList = TRUE)
@
and the appropriate points falling in any of the polygons are selected by
<<>>=
pts[pol]
@
The reverse, identical sequence of commands for 
selecting polygons \code{pol} that have 
(one or more) points of \code{pts} in them is done by
<<>>=
over(pol, pts)
over(pol, pts, returnList = TRUE)
row.names(pol[pts])
@

\section{Using \code{over} to extract attributes}
\label{attr}

This section shows how \code{over(x,y)} is used to extract attribute
values of argument \code{y} at locations of \code{x}.  The return
value is either an (aggregated) data frame, or a list.

We now create an example \code{SpatialPointsDataFrame} and a
\code{SpatialPolygonsDataFrame} using the toy data created earlier:
<<>>=
zdf = data.frame(z1 = 1:4, z2=4:1, f = c("a", "a", "b", "b"),
	row.names = c("a", "b", "c", "d"))
zdf
ptsdf = SpatialPointsDataFrame(pts, zdf)

zpl = data.frame(z = c(10, 15, 25, 3, 0), zz=1:5, 
	f = c("z", "q", "r", "z", "q"), row.names = c("x1", "x2", "x3", "x4", "x5"))
zpl
poldf = SpatialPolygonsDataFrame(pol, zpl)
@

In the simplest example 
<<>>=
over(pts, poldf)
@
a \code{data.frame} is created with each row corresponding to
the first element of the \code{poldf} attributes at locations
in \code{pts}.

As an alternative, we can pass a user-defined function to process
the table (selecting those columns to which the function makes sense):
<<>>=
over(pts, poldf[1:2], fn = mean)
@

To obtain the complete list of table entries at each point of \code{pts}, 
we use the \code{returnList} argument:
<<>>=
over(pts, poldf, returnList = TRUE)
@

The same actions can be done when the arguments are reversed:
<<>>=
over(pol, ptsdf)
over(pol, ptsdf[1:2], fn = mean)
@

\section{Lines, and Polygon-Polygon overlays with {\tt rgeos}}

Package \code{sp} provides many of the \code{over} methods, but
not all.  Package \code{rgeos} provides functions to compute
geometry intersections, i.e.  for any set of (points, lines,
polygons) to determine whether they have one ore more points in
common. The \code{over} methods not provided by package \code{sp}
can thus be completed; it involves {\em any} \code{over} methods
where a \code{SpatialLines} object is involved (either as \code{x}
or \code{y}), and the method where \code{x} and \code{y} are both of
class \code{SpatialPolygons} (table \ref{tab}). For this purpose,
objects of class \code{SpatialPixels} or \code{SpatialGrid} are
converted to \code{SpatialPolygons}. A toy example combines polygons
with lines, created by
<<>>=
l1 = Lines(Line(coordinates(pts)), "L1")
l2 = Lines(Line(rbind(c(1,1.5), c(1.5,1.5))), "L2")
L = SpatialLines(list(l1,l2))
@
and shown in figure \ref{fig:lines}.

\begin{table}
\centering
\begin{tabular}{|l|ccccc|} \hline
              & y: Points & y: Lines & y: Polygons & y: Pixels & y: Grid \\
              \hline
x: Points        &   s    &    r     &      s      &      s    &  s  \\
x: Lines         &   r    &    r     &      r      &     r:y   & r:y \\
x: Polygons      &   s    &    r     &      r      &     s:y   & s:y \\
x: Pixels        &  s:x   &   r:x    &     s:x     &     s:x   & s:x \\ \hline
x: Grid          &  s:x   &   r:x    &     s:x     &     s:x   & s:x \\ \hline
\end{tabular}
\caption{ \code{over} methods implemented for different \code{x}
and \code{y} arguments.  s: provided by \pkg{sp}; r: provided by
\pkg{rgeos}. s:x or s:y indicates that the x or y argument is converted to grid
cell center points; r:x or r:y indicate grids or pixels are converted to
polygons. }
\label{tab}
\end{table}

\begin{figure}[htb]
<<fig=TRUE,echo=FALSE>>=
plot(pol, xlim = c(-1.1, 2.1), ylim = c(-1.1, 1.6), border=2:6, axes=TRUE)
text(c(-1,0.1,0.1,1.1,0.45), c(0,0,-1.05,0,0.1), c("x1", "x2", "x3", "x4", "x5"))
lines(L, col = 'green')
text(c(0.52, 1.52), c(1.5, 1.5), c("L1", "L2"))
@
\caption{ Toy data: two lines and (overlapping) polygons (x1-x5) }
\label{fig:lines}
\end{figure}

The set of \code{over} operations on the polygons, lines and points 
is shown below (note that lists and vectors are named in this case):
<<>>=
library(rgeos)
over(pol, pol)
over(pol, pol,returnList = TRUE)
over(pol, L)
over(L, pol)
over(L, pol, returnList = TRUE)
over(L, L)
over(pts, L)
over(L, pts)
@

Another example overlays a line with a grid, shown in figure \ref{fig:grid}.
\begin{figure}
<<fig=TRUE>>=
data(meuse.grid)
gridded(meuse.grid) = ~x+y
Pt = list(x = c(178274.9,181639.6), y = c(329760.4,333343.7))
sl = SpatialLines(list(Lines(Line(cbind(Pt$x,Pt$y)), "L1")))
image(meuse.grid)
xo = over(sl, geometry(meuse.grid), returnList = TRUE)
image(meuse.grid[xo[[1]], ],col=grey(0.5),add=T)
lines(sl)
@
\caption{ Overlay of line with grid, identifying cells crossed (or touched)
by the line }
\label{fig:grid}
\end{figure}


\section{Aggregation}
In the following example, the values of a fine grid with 40 m x 40
m cells are aggregated to a course grid with 400 m x 400 m cells.
<<>>=
data(meuse.grid)
gridded(meuse.grid) = ~x+y
off = gridparameters(meuse.grid)$cellcentre.offset + 20
gt = GridTopology(off, c(400,400), c(8,11))
SG = SpatialGrid(gt)
agg = aggregate(meuse.grid[3], SG)
@
Figure \ref{fig:agg} shows the result of this aggregation
(\code{agg}, in colors) and the points (+) of the original grid
(\code{meuse.grid}). Function \code{aggregate} aggregates its first
argument over the geometries of the second argument, and returns
a geometry with attributes.  The default aggregation function
(\code{mean}) can be overridden.

\begin{figure}[htb]
<<fig=TRUE,echo=FALSE>>=
image(agg)
points(meuse.grid, pch = 3, cex=.2, col = "#80808080")
@
\caption{ aggregation over meuse.grid distance values to a 400 m
x 400 m grid}
\label{fig:agg}
\end{figure}

An example of the aggregated values of \code{meuse.grid} along
(or under) the line shown in Figure \ref{fig:line} are
<<>>=
sl.agg = aggregate(meuse.grid, sl)
class(sl.agg)
as.data.frame(sl.agg)
@
Function \code{aggregate} returns a spatial object of the same
class of \code{sl} (\code{SpatialLines}), and \code{as.data.frame}
shows the attribute table as a \code{data.frame}.

\section*{References}
\begin{itemize}
\item O'Sullivan, D., Unwin, D. (2003) Geographical Information
Analysis. Wiley, NJ.
\item
Davidson, R., 2008.  Reading topographic maps. Free e-book from:
\url{http://www.map-reading.com/}
\item
Heuvelink, G.B.M., and E.J. Pebesma, 1999.  Spatial aggregation
and soil process modelling Geoderma 89, 1-2, 47-65.
\url{http://dx.doi.org/10.1016/S0016-7061(98)00077-9}
\item
Pebesma, E., 2012.  Spatio-temporal overlay and
aggregation.  Package vignette for package spacetime,
\url{http://cran.r-project.org/web/packages/spacetime/vignettes/sto.pdf}

\end{itemize}

\end{document}
