====== Exemplo: chuva no Paraná e Krigagem ======

Exemplo: Banco de Dados de Chuva no PR\\
Banco existente no SGBD\\
(banco geográfico no formato TerraLib)

//
//
** 1. Conectando banco, visualizando no R e criando vistasno TV **

**//Fazer no TV://**
  - abrir o terraView e conectar ao banco "chuvaPR"
  - se existir alguma vista já criada, visualize!

<code R>
require(aRT)
</code>

Estabelecendo conexão com o DBMS (SGBD)
<code R>
con <- openConn()
con
</code>

Listando os bancos existentes no DBMS
<code R>
showDbs(con)
</code>

Conectando a um banco da dados
<code R>
pr <- openDb(con, "chuvaPR", update=T)
</code>

Visualizando informações deste banco
<code R>
pr
</code>

Abrindo o layer com contorno do estado
<code R>
pr.l_contorno <- openLayer(pr, "contorno")
pr.l_contorno
visualizando dado do layer
plot(pr.l_contorno)
</code>

Abrindo outro layer, das estações meteorológicas e dados de precipitação
<code R>
pr.l_dados <- openLayer(pr, "dados")
pr.l_dados
</code>

Visualizando, sobrepondo ao gráfico anterior
<code R>
plot(pr.l_dados, add=T)
</code>

Apagando uma "vista" pré-existente do banco
<code R>
deleteView(pr, "Parana")  type Y at the R prompt
</code>

**//Fazer no TV://**
  - reconectar o banco o TV e ver que a vista "sumiu"!!!

Criando uma vista/tema com o contorno ...
<code R>
tema_pol <- createTheme(pr.l_contorno, "contorno", view="Parana")
setVisual(tema_pol, visualPolygons())
</code>

... e adicionando os pontos à vista já criada
<code R>
tema_dados <- createTheme(pr.l_dados, "estacoes", view = "Parana")
setVisual(tema_dados, visualPoints(pch=21, color="black", size=5))
</code>

**// Fazer no TV://**
  - reconectar o banco o TV e ver que "apareceu" a vista agora chamada "parana" !!!

//
//
** 2. Operações da TL via aRT **

Ilustrando como (algumas) operações da terraLib podem ser usadas
<code R>
getProj(pr.l_contorno)
getMetric(pr.l_contorno, "area")
getMetric(pr.l_contorno, "length")

getDistance(pr.l_contorno, id=as.character(c(1,2)), pr.l_dados)
getDistance(pr.l_contorno, id=as.character(c(1,3)), pr.l_dados)
getDistance(pr.l_contorno, id=as.character(c(2,3)), pr.l_dados)
</code>

//
//
** 3. Trazendo e manipulando dados e fazendo análises no R **

Abrindo o tema e trazendo dados para um objeto do R\\
Duas possibilidades:
  - via "theme"
<code R>
pr
th_dados <- openTheme(pr, "estacoes")
th_dados
dados <- getData(th_dados)
head(dados)
</code>
  - via tabela do banco
<code R>
pr.l_dados
tb_dados <- openTable(pr.l_dados, "t_dados")
tb_dados
dat <- getData(tb_dados)
dat[1:5,]
pts <- getPoints(pr.l_dados)
dat[1:5,]
</code>

Trazendo o polígono com a borda do espaço para um objeto do R
<code R>
pol <- getPolygons(pr.l_contorno)
pol
plot(pol)
</code>

Fazendo agora uma interpolação dos dados de chuva via geoestatística/krigagem

  * Carregando o pacote geoR
<code R>
require(geoR)
</code>
  * Convertendo dados os dados para o formato "geodata" da geoR (conveniente) 
<code R>
geo <- as.geodata(dados, data.col=2)
geo$borders <- pol@polygons[[1]]@Polygons[[1]]@coords
plot(geo)
</code>
  * Estimando parâmetros 
<code R>
ml <- likfit(geo, trend="1st", ini=c(1000, 100)) um pouco demorado
definindo e visualizando grid de predição
loc0 <- pred_grid(geo$borders, by=10)
points(geo, bor=borders)
points(loc0, pch=".", col=2) veja gráfico!
</code>
  * e fazendo krigagem
<code R>
kc <- krige.conv(geo, loc=loc0, krige=krige.control(obj=ml), bor=geo$borders)
</code>
  * visualizando no R
<code R>
image(kc, col=terrain.colors(15), coords=parana$coords) ver gráfico
</code>

Preparing the predictions for aRT to transfer as a raster to the DBMS
<code R>
georpred <- .prepare.graph.kriging(locations=loc0, borders=parana$borders, values=kc$pred)
names(georpred)[3] <- "z"
</code>

Criando um novo layer no banco para armazenar o grid de predição
<code R>
pr.l_pred <- createLayer(pr, "predicao")
addRaster(pr.l_pred, georpred)
</code>

Checking the current status of the DB
<code R>
pr
</code>

**// Fazer no TV://**
  - reconectar o banco o TV e ver que um novo layer foi criado

Optional : setting to a TV view
<code R>
th <- createTheme(pr.l_pred, "raster", view="parana")
setVisual(th, visualRaster(color = terrain.colors(15)), mode="r")
pr
</code>

**// Fazer no TV://**
  - reconectar o banco o TV e ver que um novo tema foi adicionado a vista
  - visualize!!!!


//
//
** 4. Adicionando uma nova coluna (com predições) à tabela de dados **\\
Para exemplificar isto vamos fazer a krigagem nos pontos onde há dados
<code R>
kc0 <- krige.conv(geo, loc=geo$coords, krige=krige.control(obj=ml), borders=NULL)
prs <- data.frame(id=rownames(geo$coords), pred=kc0$pred)
createColumn(tb_dados, "preds", type="numeric")
updateColumns(tb_dados, prs)
</code>

veja o novo status da tabela
<code R>
tb_dados
</code>


**//Fazer no TV://**
  - reconectar o banco e ver que há agora uma nova tabela no layer de dados

//
//
** 5. Some plots **

Plotting directly from the data-base (not using R objects)
<code R>
plot(pr.l_pred, col=terrain.colors(15))
plot(pr.l_dados, add=T)
plot(pr.l_contorno, add=T) veja a figura!
</code>


Deletando um layer
<code R>
deleteLayer(con, "pr.l_pred")
</code>