%  aRT : API R - TerraLib                                                
%  Copyright (C) 2003-2005  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.

\documentclass{article}

\usepackage{Sweave}

\newcommand{\aRT}{{\tt aRT}}
\newcommand{\aRTgeo}{{\tt aRTgeometry}}
\newcommand{\spp}{{\tt sp}}
\newcommand{\code}[1]{{\tt #1}}
\newcommand{\tl}{\textsf{TerraLib}}

\title{Working with \spp\ and \aRT}
\author{Pedro Ribeiro de Andrade Neto\\
Paulo Justiniano Ribeiro Júnior}
\date{August 16, 2005}

\begin{document}
\maketitle

\section{Introduction}
\spp\ is an important package for exchanging information between spatial
packages. As \aRT\ manipulates all spatial data formats, it must be
fully connected to \spp. All data in \aRT\ is in \spp\ format, but 
\tl\ databases can contain data that cannot be directly converted to
\spp\ data. For example:
\begin{enumerate}
\item \aRT\ requires \code{ID} in \emph{all} spatial data, different from
\spp, that requires \code{ID} only with lines and polygons.
\item \tl\ layers have support to multigeometry, meaning that each element
can have more than one geometry associated. For example, a layer of cities
can store their contours and centroids.
\item Geometries and attributes are stored in different places in a \tl\
database. Geometries are stored directly in layers, while attributes are
stored in tables inside layers. They cannot be in the same object because
\tl\ supports different types of table (static, event, dynamics, ...).
\end{enumerate}

This document shows how to manipulate data in \aRT, showing how to import to
and read from \tl\ databases. The data (and also some sentences) used in this
document is based on \emph{S Classes and Methods for Spatial Data: the \spp\
Package}, by Pebesma and Bivand.

<<echo=TRUE,print=FALSE>>= 
library(aRT)
@

First we establish a connection to a DBMS. The database
to be used in the exemples is called ``sp''. We remove it if exists
and then we create a new one.

<<echo=TRUE,print=FALSE>>= 
aRTsilent(TRUE) # hiding info messages
con = openConn()
if(any(showDbs(con)=="sp")) deleteDb(con, "sp", force=T)
db = createDb(con, "sp")
@

<<echo=FALSE,print=FALSE>>= 
set.seed(13331)
# library(lattice)
@

\section{Spatial points}
We can generate a set of 10 points on the unit square $[0,1] \times
[0,1]$, and convert into a
{\tt SpatialPoints} object by
<<echo=TRUE,print=FALSE>>= 
xc = round(runif(10), 2)
yc = round(runif(10), 2)
xy = cbind(xc, yc)
xy.sp = SpatialPoints(xy)
xy.sp
@

To store this data we need to create a \emph{layer}, a container that can store
any geometric type and create other objects. A layer can be created using
\code{createLayer()}:
<<echo=TRUE>>=
lpoints = createLayer(db,"points")
lpoints
@
Note that we have two names, \code{"points"}, the name in the database, 
and \code{lpoints}, the R object which can access \code{"points"}.
The function \code{addPoints()} is used to store the points in a layer.
<<echo=TRUE>>=
addPoints(lpoints,xy.sp)
@

To read any data from the layer, we need to generate a table,
even if the spatial data does not have attributes (it is a \tl\
requirement).
Geometries with no entry in any table cannot be retrieved
from the database. The argument \code{genid=TRUE} at \code{createTable()} 
forces \aRT\ to generate a column and fill it with the geometries unique
ids.
<<echo=TRUE>>=
tpoints = createTable(lpoints, "tpoints", genid=TRUE)
tpoints
lpoints
@

Now the layer has 10 points and one table, and we can read the data
using \code{getGeometry}:
<<echo=TRUE>>=
getGeometry(lpoints)
@
\begin{figure}[htb]
\begin{center}
<<echo=FALSE,print=FALSE,fig=TRUE>>=
plot(lpoints, pch="x",xlab="",ylab="")
@
\end{center}
\caption{plot of a layer with points}
\label{fig:points}
\end{figure}
As the points do not have an associated geometry, we do not care about 
their order, and the order is indeed different from the original.
If we convert it to a {\tt SpatialPointsDataFrame} we can check what is
different:

<<echo=TRUE>>=
getGeometry(lpoints)
@

One way of creating a {\tt SpatialPointsDataFrame} object is by building
it from a a {\tt SpatialPoints} object and a data.frame containing
the attributes:
<<echo=TRUE>>=
df = data.frame(z1 = round(5 + rnorm(10), 2), z2 = 0:9, ID = 1:10)
xy.spdf = SpatialPointsDataFrame(xy.sp, df)
xy.spdf

lpointsdf = createLayer(db, "lpointsdf")

addPoints(lpointsdf, xy.spdf)

tpointsdf = importTable(lpointsdf, "tpointsdf", id="ID", xy.spdf)
tpointsdf
getGeometry(lpointsdf)
getData(tpointsdf)
@

As \code{SpatialPointsDataFrame} does not have \code{ID}, a
simple merge of the layer with the table will not work\footnote{E-mail
para eles, as outras classes tem id e fazem o merge, 
entao acho que nao é tarefa do aRT.}.
<<echo=TRUE>>=
#SpatialPointsDataFrame(getGeometry(lpointsdf), getData(tpointsdf))

#tpointsdf = importTable(lpoints, "tpointsdf", xy.spdf)

#SpatialPointsDataFrame(getGeometry(lpointsdf), getData(tpointsdf))
@

\section{Grids}


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

<<fig=FALSE>>=
l1 = cbind(c(1,2,3),c(3,2,2))
l1a = cbind(l1[,1]+.05,l1[,2]+.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, Sl1a), ID="a")
S2 = Lines(list(Sl2), ID="b")
Sl = SpatialLines(list(S1,S2))

llines = createLayer(db,"llines")

addLines(llines, Sl)

createTable(llines, "llines", gen=T)
@

\subsection{Building line objects with attributes}

The same as polygons

\section{Polygons}
<<echo=FALSE,print=FALSE>>= 
set.seed(25)
@
\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.
<<fig=FALSE>>=

#genPol = function(quant)
#{
#    res = list()
#    res = lapply(1:quant, function(x) {
#        from = round(runif(2,1,20),2)
#        to   = from + round(runif(2,1,5),2)
#
#        Sr1 = Polygon(cbind(c(from[1],from[1],to[1],to[1],  from[1]),
#                          c(from[2],to[2],  to[2],from[2],from[2])))
#        Srs1 = Polygons(list(Sr1),paste(x))
#    })
#    SR = SpatialPolygons(res, 1:quant)
#}

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", genid=T)#FALSE)

#th=createTheme(lrings, "trings")
lrings

#sr = as.SpatialPolygon(getGeometry(lrings))

#getGeometry(lrings)
@


\begin{figure}
\begin{center}
<<echo=FALSE,print=FALSE,fig=TRUE>>=
plot(lrings, col="gray", xlab="",ylab="")
@
\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):
<<>>=
attr = data.frame(a=1:3, b=3:1, row.names=c("s1","s2","s3/4"))
SrDf = SpatialPolygonsDataFrame(SR, attr)
lringsdf = createLayer(db, "lringsdf")
addPolygons(lringsdf, SrDf)
@

To import the attributes, we need to create a table, but, due to the internal
differences of \spp\ data storage
<<fig=FALSE>>=
xy.spdf@data
SrDf@data
class(SrDf@data)
@
we need to insert \code{SrDf} manually, creating both table and 
the two integer columns before inserting the data:
<<fig=FALSE>>=
tringsdf = createTable(lringsdf, "tringsdf", id="ID", gen=F)
createColumn(tringsdf, "a", "i")
createColumn(tringsdf, "b", "i")
addRows(tringsdf, SrDf@data)
tringsdf
#lringsdf

getData(tringsdf)

#sr = as.SpatialPolygons(getGeometry(lringsdf))
#summary(sr)
#sr = as.SpatialPolygonsDataFrame(getGeometry(lringsdf))
#sr
@

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

\end{document}

% vim:syntax=tex
