%  aRT : API R - TerraLib                                                
%  Copyright (C) 2003-2008  LEG                                          
%                                                                        
%  This program is free software; you can redistribute it and/or modify  
%  it under the terms of the GNU General Public License as published by  
%  the Free Software Foundation; either version 2 of the License, or     
%  (at your option) any later version.                                   
%                                                                        
%  This program is distributed in the hope that it will be useful,       
%  but WITHOUT ANY WARRANTY; without even the implied warranty of        
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         
%  GNU General Public License for more details.                          
%                                                                        
%  You should have received a copy of the GNU Lesser General Public      
%  License along with this library.

% \VignetteIndexEntry{Working with sp and aRT}

\documentclass{article}

\usepackage{Sweave}
\usepackage{aRT}

\title{Working with \spp\ and \aRT}
\author{Pedro Ribeiro de Andrade\\
Paulo Justiniano Ribeiro J\'{u}nior}
\date{July 20, 2009}

\begin{document}
\maketitle
\tableofcontents

\section{Introduction}
\spp\ (Pebesma \& Bivand, 2005) 
is an important package that defines standards and allows 
for exchanging information between spatial
packages in \R. 
As \aRT\ manipulates all spatial data formats, it was designed to 
follow \spp\ standards to represent spatial data,
storing and retrieving data using the \spp\ format.
Therefore, when working with \aRT\ and spatial data in \R{}
it is required that objects containing spatial data 
are converted, whenever necessary, to the \spp\ format.
For futher details on \spp\ package
see \texttt{http://r-spatial.sourceforge.net}.


On the other hand,  
\tl\ databases can contain data that cannot be directly converted to
\spp\ data. For example:
\begin{enumerate}
\item \tl\ (and therefore \aRT) requires \code{ID} in \emph{all} spatial data,
different from \spp, that requires \code{ID} only for lines and polygons.
\item \tl\ layers have support to multigeometry, meaning that each spatial element
can have more than one geometry associated. For example, a layer of cities
can store both contours and centroids.
\item geometries and attributes are stored in different objects in a \tl\
database. Geometries are stored directly inside layers, whereas attributes are
stored in tables inside layers. The reason why tables  
cannot be in the same object as geometries is because
\tl\ supports different table formats, for example \textit{static},  \textit{event}, 
and  \textit{dynamic}.
\end{enumerate}

This document illustrates how to manipulate spatial data in \aRT, 
showing how to import to and read from \tl\ databases. 
The data (and also some sentences!)
of this document are extracted from Pebesma and Bivand (2005).

%<<echo=TRUE,print=FALSE,label="library">>= 
%require(aRT)
%@

We start by loading the package and, for convenience here, 
setting \aRT\ to the silent mode which hides some information messages issued by \aRT\ functions; followed by establishing a connection to a DBMS. 
The database to be used in the examples is called ``sp'', which, if exists, 
is removed from the DBMS and recreated.
Notice this is a \tl database, i.e., a database created following the \tl\ structure.

\begin{Schunk}
\begin{Sinput}
> require(aRT)
> aRTsilent(TRUE)
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\begin{Sinput}
> con = openConn(name = "root")
> if (any(showDbs(con) == "sp")) deleteDb(con, "sp", force = T)
> db = createDb(con, "sp")
\end{Sinput}
\end{Schunk}

Further details on connections to the DBMS are provided by the 
vignette \textsl{aRTconn}.


\section{Spatial points}

\subsection{Points without attributes}
Our first example illustrates how to handle data with a \textit{points geometry}. 
For instance, consider a set of 10 points randomly generated 
on the unit square $[0,1] \times [0,1]$ and stored in a matrix
\code{xy}. The first step is to use \spp\ 
to convert into a \code{SpatialPoints} object.
\begin{Schunk}
\begin{Sinput}
> xy = matrix(runif(20), nc = 2)
> xy.sp = SpatialPoints(xy)
> xy.sp
\end{Sinput}
\begin{Soutput}
SpatialPoints:
      coords.x1 coords.x2
 [1,] 0.6739453 0.3021973
 [2,] 0.9619853 0.6170330
 [3,] 0.9228386 0.9050982
 [4,] 0.7747456 0.8476483
 [5,] 0.7248062 0.4579902
 [6,] 0.7413480 0.3863993
 [7,] 0.6302933 0.6443511
 [8,] 0.1903071 0.7170724
 [9,] 0.6977443 0.1980982
[10,] 0.3703263 0.2772779
Coordinate Reference System (CRS) arguments: NA 
\end{Soutput}
\end{Schunk}
However, an \code{SpatialPoints} object cannot be used by \aRT\ functions 
because it does not have \code{ID} and it is necessary to
convert to a \code{SpatialPointsDataFrame}.

\begin{Schunk}
\begin{Sinput}
> xy.spdf = SpatialPointsDataFrame(xy, data.frame(ID = paste(1:10)))
\end{Sinput}
\end{Schunk}

Note that here you can use both \code{xy} or \code{xy.spdf} as argument to \code{Spatial PointsDataFrame}.

The first step to store this data in a \tl\ database 
is to create a \emph{layer}. A layer is a container
that can store any geometric type and also, optionally, other types of objects. 
A layer can be created in a database using \code{createLayer()} and the function
call receives as arguments the database object and a
string defining the name to be set for the layer.
\begin{Schunk}
\begin{Sinput}
> lpoints = createLayer(db, "points")
> lpoints
\end{Sinput}
\begin{Soutput}
Object of class aRTlayer

Layer: "points"
Database: "sp"
Layer is empty
Projection Name: "NoProjection"
Projection Datum: "Spherical"
Tables: (none) 
\end{Soutput}
\end{Schunk}
Notice that we have two names, \code{"points"}
which is the name of the layer in the database, 
and \code{lpoints}, an \R{} object which can access \code{"points"}.
Next the function \code{addPoints()} is used to store the points into the layer.
Notice that after that the layer object acknowledges the points are now
stored into the database.
\begin{Schunk}
\begin{Sinput}
> addPoints(lpoints, xy.spdf)
> lpoints
\end{Sinput}
\begin{Soutput}
Object of class aRTlayer

Layer: "points"
Database: "sp"
Number of points: 10
Projection Name: "NoProjection"
Projection Datum: "Spherical"
Tables: (none) 
\end{Soutput}
\end{Schunk}
To conclude storing the points a further step creating adding a table
to the layer is still necessary
in order to  ba able to read any data from the layer.
This is a \tl\ requirement needed even when the spatial data does not have
any attributes. 
Geometries with no entry in any table cannot be retrieved
from the database.
\begin{Schunk}
\begin{Sinput}
> tpoints = createTable(lpoints, "tpoints")
> tpoints
\end{Sinput}
\begin{Soutput}
Object of class aRTtable

Table: "tpoints"
Type: static
Layer: "points"
Rows: 10
Attributes: 
    ID: character[16] (key)
\end{Soutput}
\begin{Sinput}
> lpoints
\end{Sinput}
\begin{Soutput}
Object of class aRTlayer

Layer: "points"
Database: "sp"
Number of points: 10
Projection Name: "NoProjection"
Projection Datum: "Spherical"
Tables:
    "tpoints": static
\end{Soutput}
\end{Schunk}

Now the layer has 10 points and one table, and we can retrieve the 
point coordinates using \code{getPoints()}, which returns an 
object of class \code{SpatialPointsDataFrame}:
\begin{Schunk}
\begin{Sinput}
> points = getPoints(lpoints)
> points
\end{Sinput}
\begin{Soutput}
            coordinates ID
1  (0.673945, 0.302197)  1
2  (0.370326, 0.277278) 10
3  (0.961985, 0.617033)  2
4  (0.922839, 0.905098)  3
5  (0.774746, 0.847648)  4
6   (0.724806, 0.45799)  5
7  (0.741348, 0.386399)  6
8  (0.630293, 0.644351)  7
9  (0.190307, 0.717072)  8
10 (0.697744, 0.198098)  9
\end{Soutput}
\end{Schunk}

Note that the points have a different order from the original data. That is because
the database stores the IDs as characters, therefore 10 comes before 2.

Figure~\ref{fig:points} 
illustrates two different ways for
visualising the point coordinates, on the left by
plotting the points from an \R{} object with \code{plot(points)}, 
or  directly from the layer, with \code{plot(lpoints)},
which does not requires the existence of an \R{} object containing
the points.
\begin{figure}[htb]
\begin{center}
\includegraphics{aRTsp-010}
\end{center}
\caption{Plot of a layer with points}
\label{fig:points}
\end{figure}

\subsection{Points with attributes}
Another possible way of creating a \code{SpatialPointsDataFrame} 
object is by building it by combining a \code{SpatialPoints} 
object and a data frame containing associated attributes.
In the following example we combine the object \code{xy.sp}
with a data frame \code{df} containg attributes.
Notice the \code{ID} column is required by any \code{sp} object.
\begin{Schunk}
\begin{Sinput}
> df = data.frame(z1 = round(5 + rnorm(10), 2), z2 = 0:9, 
+     ID = paste(1:10))
> xy.spdf = SpatialPointsDataFrame(xy.sp, df)
> xy.spdf
\end{Sinput}
\begin{Soutput}
            coordinates   z1 z2 ID
1  (0.673945, 0.302197) 3.10  0  1
2  (0.961985, 0.617033) 4.15  1  2
3  (0.922839, 0.905098) 3.68  2  3
4  (0.774746, 0.847648) 4.45  3  4
5   (0.724806, 0.45799) 6.62  4  5
6  (0.741348, 0.386399) 5.57  5  6
7  (0.630293, 0.644351) 3.66  6  7
8  (0.190307, 0.717072) 3.75  7  8
9  (0.697744, 0.198098) 5.19  8  9
10 (0.370326, 0.277278) 5.02  9 10
\end{Soutput}
\end{Schunk}
As before we create a layer and aad the points to it. Next, 
as our object now has attributes, we can import the table data
using \code{importTable()}.
\begin{Schunk}
\begin{Sinput}
> lpointsdf = createLayer(db, "lpointsdf")
> addPoints(lpointsdf, xy.spdf)
> tpointsdf = importTable(lpointsdf, "tpointsdf", ID = "ID", 
+     xy.spdf)
> tpointsdf
\end{Sinput}
\begin{Soutput}
Object of class aRTtable

Table: "tpointsdf"
Type: static
Layer: "lpointsdf"
Rows: 10
Attributes: 
    ID: character[16] (key)
    z1: numeric
    z2: integer
\end{Soutput}
\begin{Sinput}
> lpointsdf
\end{Sinput}
\begin{Soutput}
Object of class aRTlayer

Layer: "lpointsdf"
Database: "sp"
Number of points: 10
Projection Name: "NoProjection"
Projection Datum: "Spherical"
Tables:
    "tpointsdf": static
\end{Soutput}
\end{Schunk}

When retrieving data from the database to \R{} 
getting point coordinates and the table at once from the layer 
we can use a second argument of \code{getPoints()} with the table to be read.
\begin{Schunk}
\begin{Sinput}
> getPoints(lpointsdf, tpointsdf)
\end{Sinput}
\begin{Soutput}
            coordinates ID   z1 z2
1  (0.673945, 0.302197)  1 3.10  0
2  (0.370326, 0.277278) 10 5.02  9
3  (0.961985, 0.617033)  2 4.15  1
4  (0.922839, 0.905098)  3 3.68  2
5  (0.774746, 0.847648)  4 4.45  3
6   (0.724806, 0.45799)  5 6.62  4
7  (0.741348, 0.386399)  6 5.57  5
8  (0.630293, 0.644351)  7 3.66  6
9  (0.190307, 0.717072)  8 3.75  7
10 (0.697744, 0.198098)  9 5.19  8
\end{Soutput}
\end{Schunk}

\subsection{Doing all at once}

All the steps above can be encapsulated using
\code{importSpData()}.

\section{Grids}
(not supported yet)

\section{Lines}

\subsection{Building line objects from scratch}
In many instances, line coordinates will be retrieved from external
sources. The following example shows how to build an object of class
{\tt SpatialLines} from scratch. As objects from this class already stores
\code{ID}, they are pushed in the layer directly using \code{addLines()}.

\begin{Schunk}
\begin{Sinput}
> l1 = cbind(c(1, 2, 3), c(3, 2, 2))
> l1a = cbind(l1[, 1] + 0.05, l1[, 2] + 0.05)
> l2 = cbind(c(1, 2, 3), c(1, 1.5, 1))
> Sl1 = Line(l1)
> Sl1a = Line(l1a)
> Sl2 = Line(l2)
> S1 = Lines(list(Sl1), ID = "a")
> S2 = Lines(list(Sl2), ID = "b")
> S3 = Lines(list(Sl1a), ID = "c")
> Sl = SpatialLines(list(S1, S2, S3))
> llines = createLayer(db, "llines")
> addLines(llines, Sl)
> createTable(llines, "llines")
\end{Sinput}
\begin{Soutput}
Object of class aRTtable

Table: "llines"
Type: static
Layer: "llines"
Rows: 3
Attributes: 
    ID: character[16] (key)
\end{Soutput}
\end{Schunk}

\begin{figure}[htb]
\begin{center}
\includegraphics{aRTsp-015}
\end{center}
\caption{Plot of a layer with lines}
\label{fig:points}
\end{figure}


\subsection{Building line objects with attributes}

The same as polygons

\section{Polygons}
\subsection{Building from scratch}
The following example shows how a set of polygons are built from scratch.
Note that {\tt Sr4} has the opposite direction (right) as the other three;
it is meant to represent a hole in the {\tt Sr3} polygon.
\begin{Schunk}
\begin{Sinput}
> Sr1 = Polygon(cbind(c(2, 4, 4, 1, 2), c(2, 3, 5, 4, 
+     2)))
> Sr2 = Polygon(cbind(c(5, 4, 2, 5), c(2, 3, 2, 2)))
> Sr3 = Polygon(cbind(c(4, 4, 5, 10, 4), c(5, 3, 2, 5, 
+     5)))
> Sr4 = Polygon(cbind(c(5, 6, 6, 5, 5), c(4, 4, 3, 3, 
+     4)), hole = TRUE)
> Srs1 = Polygons(list(Sr1), "s1")
> Srs2 = Polygons(list(Sr2), "s2")
> Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
> SR = SpatialPolygons(list(Srs1, Srs2, Srs3), 1:3)
> lrings = createLayer(db, "lrings")
> addPolygons(lrings, SR)
> trings = createTable(lrings, "trings")
> lrings
\end{Sinput}
\begin{Soutput}
Object of class aRTlayer

Layer: "lrings"
Database: "sp"
Number of polygons: 4
Projection Name: "NoProjection"
Projection Datum: "Spherical"
Tables:
    "trings": static
\end{Soutput}
\begin{Sinput}
> pols = getPolygons(lrings)
\end{Sinput}
\end{Schunk}

\begin{figure}
\begin{center}
\includegraphics{aRTsp-018}
\end{center}
\caption{plot of a layer with polygons}
\label{fig:polygons}
\end{figure}

\subsection{Polygons with attributes}
Polygons with attributes, objects of class {\tt SpatialPolygonsDataFrame},
are built from the {\tt SpatialPolygons} object (topology) and the attributes
(data.frame):

To import the attributes, we need to create a table, but, due to the internal
differences of \spp\ data storage
we need to insert \code{SrDf} manually, creating both table and 
the two integer columns before inserting the data:

\section*{References}
\begin{description}
\item Chambers, J.M., 1998, Programming with data, a guide to the S language.
Springer, New York.

\item Pebesma, E.J.  and Bivand, R.S., 2005, Classes and methods for spatial data in R, 
R-News 5 (2), pp. 9-13.

\end{description}

\end{document}

% vim:syntax=tex
