Exemplo: chuva no Paraná e Krigagem

1. exemplo: Banco de Dados de Chuva no PR Banco existente no SGBD (banco geográfico no formato TerraLib) ## 2. Conectando e inspecionando o banco ##

para ilustrar o uso vamos fazer as operações no R via aRT e acompanhar visualizando no terraView

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

fazer no TV: a. abrir o terraView e conectar ao banco "chuvaPR" ## b. se existir alguma vista já criada, visualize!

require(aRT)

estabelecendo conexão com o DBMS (SGBD) con ← openConn() con

listando os bancos existentes no DBMS showDbs(con)

conectando a um banco da dados pr ← openDb(con, "chuvaPR", update=T)

visualizando informações deste banco pr

abrindo o layer com contorno do estado pr.l_contorno ← openLayer(pr, "contorno") pr.l_contorno visualizando dado do layer plot(pr.l_contorno)

abrindo outro layer, das estações meteorológicas e dados de precipitação pr.l_dados ← openLayer(pr, "dados") pr.l_dados

visualizando, sobrepondo ao gráfico anterior plot(pr.l_dados, add=T)

apagando uma "vista" pré-existente do banco deleteView(pr, "Parana") type Y at the R prompt

fazer no TV: reconectar o banco o TV e ver que a vista "sumiu"!!!

criando uma vista/tema com o contorno … tema_pol ← createTheme(pr.l_contorno, "contorno", view="Parana") setVisual(tema_pol, visualPolygons())

… e adicionando os pontos à vista já criada tema_dados ← createTheme(pr.l_dados, "estacoes", view = "Parana") setVisual(tema_dados, visualPoints(pch=21, color="black", size=5))

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 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)

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

abrindo o tema e trazendo dados para um objeto do R

duas possibilidades: a. via "theme" pr th_dados ← openTheme(pr, "estacoes") th_dados dados ← getData(th_dados) head(dados)

b. via tabela do banco 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,]

trazendo o polígono com a borda do espaço para um objeto do R pol ← getPolygons(pr.l_contorno) pol plot(pol)

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

carregando o pacote geoR require(geoR)

convertendo dados os dados para o formato "geodata" da geoR (conveniente) geo ← as.geodata(dados, data.col=2) geo$borders ← pol@polygons1@Polygons1@coords plot(geo)

estimando parâmetros 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! e fazendo krigagem kc ← krige.conv(geo, loc=loc0, krige=krige.control(obj=ml), bor=geo$borders)

visualizando no R image(kc, col=terrain.colors(15), coords=parana$coords) ver gráfico

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

criando um novo layer no banco para armazenar o grid de predição pr.l_pred ← createLayer(pr, "predicao") addRaster(pr.l_pred, georpred)

checking the current status of the DB pr

fazer no TV: reconectar o banco o TV e ver que um novo layer foi criado

Optional : setting to a TV view

th ← createTheme(pr.l_pred, "raster", view="parana") setVisual(th, visualRaster(color = terrain.colors(15)), mode="r") pr

fazer no TV: reconectar o banco o TV e ver que um novo tema foi adicionado a vista visualize!!!!

## Adicionando uma nova coluna (com predições) à tabela de dados ##

para exemplicifar isto vamos fazer a krigagem nos pontos onde há dados 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)

veja o novo status da tabela tb_dados

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

## 4. Some plots ##

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

## Deletando um layer

# deleteLayer(con, "pr.l_pred")